EdbDataProc Class Reference

emulsion data processing More...

#include <EdbDataSet.h>

Inheritance diagram for EdbDataProc:

Collaboration diagram for EdbDataProc:

Public Member Functions
void	AjustZ (int doZ)
void	Align (int alignFlag)
void	AlignLinkTracks (int alg=0, int alignFlag=0)
int	CheckCCD ()
int	CheckShrinkage (EdbPVRec *ali, int couple, float &shr1, float &shr2)
void	CloseCouplesTree (TTree *tree)
void	CorrectAngles ()
	EdbDataProc ()
	============================================================================== More...
	EdbDataProc (const char *file)
	============================================================================== More...
	EdbDataProc (int npl, TArrayI &ids, TArrayF &zs)
EdbPVRec *	ExtractDataVolume (EdbSegP &v, int plmin, int plmax, float accept[4], int datatype=0)
EdbPVRec *	ExtractDataVolume (EdbTrackP &tr, float binx=20, float bint=10, int datatype=0)
EdbPVRec *	ExtractDataVolumeF (EdbTrackP &tr, float binx=20, float bint=10, int datatype=0)
void	FillCouplesTree (TTree *tree, EdbPVRec *al, int fillraw=0)
void	FineAlignment (int doFine)
void	FineAlignmentTracks ()
EdbDataSet *	GetDataSet ()
EdbPVRec *	GetPVR () const
EdbVertex *	GetVertexFromTree (EdbVertexRec &vertexrec, const char *fname, const int vertexID)
int	InitVolume (EdbPVRec *ali, int datatype=0, TIndex2 *trseg=0)
int	InitVolume (int datatype=0, const char *rcut="1")
int	InitVolumeRaw (EdbPVRec *ali)
int	InitVolumeTracks (EdbPVRec *ali, const char *rcut)
int	Link ()
int	Link (EdbDataPiece &piece)
void	LinkRawTracks (int alg=0)
void	LinkTracks (int alg=0, float p=-1.)
void	LinkTracksC (int alg=0, float p=-1.)
TIndex2 *	MakeTracksSegmentsList (EdbPVRec &ali)
int	NoUpdate () const
int	Process ()
EdbPVRec *	PVR () const
void	SetNoUpdate (int nu)
void	SetPVR (EdbPVRec *pvr)
int	ShrinkCorr ()
virtual	~EdbDataProc ()

Static Public Member Functions
static int	LinkTracksWithFlag (EdbPVRec *ali, float p, float probmin, int nsegmin, int maxgap, int flag, float mass=0.1396)
static int	MakeTracksTree (EdbPVRec *ali=0, const char *file="linked_tracks.root")
static int	MakeTracksTree (TObjArray &tracks, float xv=0, float yv=0, const char *file="linked_tracks.root")
static int	MakeVertexTree (TObjArray &vtxarr, const char *file)
static int	ReadTracksTree (EdbPVRec &ali, const char *fname="linked_tracks.root", const char *rcut="t.eFlag>-1&&nseg>2&&t.eProb>.01")
static int	ReadVertexTree (EdbVertexRec &vertexrec, const char *fname, const char *rcut, map< int, EdbTrackP * > trackID_map)
static int	ReadVertexTree (EdbVertexRec &vertexrec, const char *fname, const char *rcut, TObjArray *builttracks=0)

Private Attributes
EdbDataSet *	eDataSet
int	eNoUpdate
EdbPVRec *	ePVR

Detailed Description

emulsion data processing

Constructor & Destructor Documentation

EdbDataProc() [1/3]

EdbDataProc::EdbDataProc

(

)

==============================================================================

{
  eDataSet=0;
  eNoUpdate=1;
  ePVR     = new EdbPVRec();
}

EdbDataProc() [2/3]

EdbDataProc::EdbDataProc	(	int	npl,
		TArrayI &	ids,
		TArrayF &	zs
	)

---

Special constructor (requested by Nicolai) set everything as default as possible

{
 
  eDataSet  = new EdbDataSet();
  ePVR      = 0;
  eNoUpdate = 0;
 
  if( npl>ids.GetSize() )  return;
  if( npl>zs.GetSize() )   return;
 
  ePVR = new EdbPVRec();
  EdbPattern *pat =0;
  for(int i=0; i<npl; i++) {
    pat = new EdbPattern( 0,0, zs.At(i) );
    pat->SetPID(ids.At(i));
    ePVR->AddPattern(pat);
  }
  ePVR->SetPatternsID();
 
  ePVR->SetScanCond( new EdbScanCond() );
  ePVR->SetCouplesAll();
}

EdbDataProc() [3/3]

EdbDataProc::EdbDataProc

(

const char *

file

)

==============================================================================

{
  eDataSet = new EdbDataSet(file);
  ePVR     = 0;
  eNoUpdate=0;
}

~EdbDataProc()

EdbDataProc::~EdbDataProc ( )

virtual

---

{
  if(eDataSet) delete eDataSet;
  if(ePVR)     delete ePVR;
}

Member Function Documentation

AjustZ()

void EdbDataProc::AjustZ

(

int

doZ

)

---

{
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolume(ali);
  ali->Link();
  ali->FillTracksCell();
  ali->SelectLongTracks(doZ);
  ali->MakeSummaryTracks();
  ali->FineCorrZnew();
  for(int i=0; i<ali->Npatterns(); i++) 
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateZPar(0,ali->GetPattern(i)->Z());
}

Align()

void EdbDataProc::Align

(

int

doAlign

)

---

{
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolume(ali);
 
  ali->Align(doAlign);
 
  if(gEDBDEBUGLEVEL>1) ali->PrintAff();
  EdbAffine2D  aff;
  for(int i=0; i<ali->Npatterns(); i++) {
    ali->GetPattern(i)->GetKeep(aff);
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffPar(0,aff);
  }
}

AlignLinkTracks()

void EdbDataProc::AlignLinkTracks	(	int	alg = 0,
		int	doAlign = 0
	)

---

{
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolume(ali);
 
  ali->Align(doAlign);
  if(gEDBDEBUGLEVEL>1) ali->PrintAff();
  EdbAffine2D  aff;
  for(int i=0; i<ali->Npatterns(); i++) {
    ali->GetPattern(i)->GetKeep(aff);
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffPar(0,aff);
  }
 
  ali->Link();
  Log(2,"EdbDataProc::AlignLinkTracks","link ok");
 
//    if(alg==1) {
//      ali->MakeHoles();
//      ali->Link();
//      printf("link ok\n");
//    } 
 
  ali->FillTracksCell();
 
  TTree *cptree=EdbDataPiece::InitCouplesTree("linked_couples.root","RECREATE");
  FillCouplesTree(cptree, ali,0);
  CloseCouplesTree(cptree);
 
  MakeTracksTree(ali);
}

CheckCCD()

int EdbDataProc::CheckCCD

(

)

---

{
  if(!eDataSet) return 0;
  EdbDataPiece *piece;
  int ndef=0;
  int np=eDataSet->N();
  for(int i=0; i<np; i++) {
    piece = eDataSet->GetPiece(i);
    ndef = piece->CheckCCD();
    if(ndef<0) { Log(2,"EdbDataProc::CheckCCD","skip piece"); continue; }
    Log(1,"EdbDataProc::CheckCCD","piece %s: eliminated defects: %d",piece->GetName(),ndef); 
    piece->WriteCuts();
    piece->CloseRun();
  }
  return np;
}

CheckShrinkage()

int EdbDataProc::CheckShrinkage	(	EdbPVRec *	ali,
		int	couple,
		float &	shr1,
		float &	shr2
	)

---

{
  EdbPatCouple  *patc = ali->GetCouple(couple);
  if(!patc) return 0;
 
  EdbSegCouple *sc=0;
  EdbSegP *s1=0, *s2=0;
  double dz,tx,ty,t,t1,t2,sumt1=0,sumt2=0;
  int    nsum=0;
 
  int  nc=patc->Ncouples();
  for( int ic=0; ic<nc; ic++ ) {
    sc = patc->GetSegCouple(ic);
 
    if(sc->CHI2()>1.5)  continue;
 
    s1 = patc->Pat1()->GetSegment(sc->ID1());
    s2 = patc->Pat2()->GetSegment(sc->ID2());
 
    dz = s2->Z() - s1->Z();
    tx = (s2->X() - s1->X())/dz;
    ty = (s2->Y() - s1->Y())/dz;
 
    t  = TMath::Sqrt( tx*tx + ty*ty );
    if(t<.1) continue;
    if(t>.45) continue;
    t1 = TMath::Sqrt( s1->TX()*s1->TX() + s1->TY()*s1->TY() );
    t2 = TMath::Sqrt( s2->TX()*s2->TX() + s2->TY()*s2->TY() );
 
    nsum++;
    sumt1 += t1/t;
    sumt2 += t2/t;
  }
 
  if(nsum<1)  return 0;
 
  shr1 = sumt1/nsum;
  shr2 = sumt2/nsum;
 
  return nsum;
}

CloseCouplesTree()

void EdbDataProc::CloseCouplesTree

(

TTree *

tree

)

---

{
  tree->AutoSave();
  TFile *f=0;
  f = tree->GetCurrentFile();
  if (f) SafeDelete(f);
  tree=0;
}

CorrectAngles()

void EdbDataProc::CorrectAngles

(

)

---

{
  EdbDataPiece *piece;
  int npieces = eDataSet->N();
  Log(2,"EdbDataProc::CorrectAngles","npieces = %d",npieces);
  if(!npieces) return;
 
  for(int i=0; i<npieces; i++ ) {
    piece = eDataSet->GetPiece(i);
    if(piece) piece->CorrectAngles();
  }
}

ExtractDataVolume() [1/2]

EdbPVRec * EdbDataProc::ExtractDataVolume	(	EdbSegP &	seg,
		int	plmin,
		int	plmax,
		float	acc[4],
		int	datatype = 0
	)

---

datatype: 0 - base track
1 - up
2 - down
100 - raw?

{
 
  if(!ePVR)               return 0;
  int npat = ePVR->Npatterns();
 
  EdbPVRec *ali = new EdbPVRec();
 
  Log(2,"EdbDataProc::ExtractDataVolume","Select Data %d in EdbPVRec with %d patterns:",
     datatype, npat);
 
  float dz, min[5],  max[5];
  EdbPattern *pat=0;
 
  for(int i=0; i<npat; i++) {
 
    pat = ePVR->GetPattern(i);
    if(!pat)                   continue;
 
    dz = pat->Z() - seg.Z();
 
    min[0] = seg.X() + dz*seg.TX() - acc[0];
    min[1] = seg.Y() + dz*seg.TY() - acc[1];
    max[0] = seg.X() + dz*seg.TX() + acc[0];
    max[1] = seg.Y() + dz*seg.TY() + acc[1];
 
    min[2] = -.6;    max[2] = .6;
    min[3] = -.6;    max[3] = .6;
    min[4] =  0.;    max[4] = 100.;
 
    ali->AddPattern( pat->ExtractSubPattern(min,max) );
  }
  
  return ali;
}

ExtractDataVolume() [2/2]

EdbPVRec * EdbDataProc::ExtractDataVolume	(	EdbTrackP &	tr,
		float	binx = 20,
		float	bint = 10,
		int	datatype = 0
	)

---

datatype: 0 - base track
1 - up
2 - down
100 - raw?

{
 
  if(!ePVR)               return 0;
  int npat = ePVR->Npatterns();
 
  EdbPVRec *ali = new EdbPVRec();
 
  Log(2,"EdbDataProc::ExtractDataVolume","Select Data %d in EdbPVRec with %d patterns",
     datatype, npat);
 
  float min[5],  max[5];
 
  float dx  = binx*TMath::Sqrt(tr.SX()*tr.N());
  float dy  = binx*TMath::Sqrt(tr.SY()*tr.N());
  float dtx = bint*TMath::Sqrt(tr.STX()*tr.N());
  float dty = bint*TMath::Sqrt(tr.STY()*tr.N());
 
  float dz  = (tr.GetSegment(tr.N()-1)->Z() - tr.GetSegment(0)->Z());
  float tx  = (tr.GetSegment(tr.N()-1)->X() - tr.GetSegment(0)->X())/dz;
  float ty  = (tr.GetSegment(tr.N()-1)->Y() - tr.GetSegment(0)->Y())/dz;
 
  Log(2,"EdbDataProc::ExtractDataVolume","select segments with acceptance: %f %f [microns]   %f %f [mrad]",
     dx,dy,dtx*1000,dty*1000);
 
  EdbPattern *pat  = 0;
  EdbPattern *spat = 0;
  int nseg =0;
 
  for(int i=0; i<npat; i++) {
 
    pat = ePVR->GetPattern(i);
    if(!pat)                   continue;
 
    dz = pat->Z() - tr.Z() + 107.;        //TODO
 
    min[0] = tr.X() + dz*tx - dx;
    max[0] = tr.X() + dz*tx + dx;
    min[1] = tr.Y() + dz*ty - dy;
    max[1] = tr.Y() + dz*ty + dy;
    min[2] = tx - dtx;
    max[2] = tx + dtx;
    min[3] = ty - dty;
    max[3] = ty + dty;
    min[4] =  0.;    
    max[4] = 100.;
 
    spat = pat->ExtractSubPattern(min,max);
    nseg += spat->N();
    ali->AddPattern( spat );
  }
  Log(2,"EdbDataProc::ExtractDataVolume","%d segments are selected",nseg);
  return ali;
}

ExtractDataVolumeF()

EdbPVRec * EdbDataProc::ExtractDataVolumeF	(	EdbTrackP &	tr,
		float	binx = 20,
		float	bint = 10,
		int	datatype = 0
	)

---

{
  if(!ePVR)               return 0;
  int npat = ePVR->Npatterns();
  EdbPVRec *ali = new EdbPVRec();
  Log(2,"EdbDataProc::ExtractDataVolumeF","Select Data %d in EdbPVRec with %d patterns:",
     datatype, npat);
 
  EdbSegP ss; // the "selector" segment 
  ss.SetCOV( tr.GetSegment(0)->COV() );
 
  float dz  = (tr.GetSegment(tr.N()-1)->Z() - tr.GetSegment(0)->Z());
  float tx  = (tr.GetSegment(tr.N()-1)->X() - tr.GetSegment(0)->X())/dz;
  float ty  = (tr.GetSegment(tr.N()-1)->Y() - tr.GetSegment(0)->Y())/dz;
 
  ss.SetTX(tx);
  ss.SetTY(ty);
  ss.SetX(tr.X());
  ss.SetY(tr.Y());
  ss.SetZ(tr.Z());
 
  if(gEDBDEBUGLEVEL>2) ss.Print();
 
  EdbPattern *pat  = 0;
  EdbPattern *spat = 0;
  int nseg =0;
 
  TObjArray arr;
 
  for(int i=0; i<npat; i++) {
    pat = ePVR->GetPattern(i);
    if(!pat)                   continue;
    ss.PropagateTo(pat->Z());
    nseg += pat->FindCompliments(ss,arr,binx,bint);
  }
 
  spat = new EdbPattern(0,0,0);
  for(int i=0; i<arr.GetEntriesFast(); i++)
    spat->AddSegment( *((EdbSegP*)(arr.At(i))) );
 
  ali->AddPattern( spat );
 
  Log(2,"EdbDataProc::ExtractDataVolumeF","%d segments are selected",nseg);
  return ali;
}

FillCouplesTree()

void EdbDataProc::FillCouplesTree	(	TTree *	tree,
		EdbPVRec *	al,
		int	fillraw = 0
	)

---

{
  tree->GetDirectory()->cd();
 
  Log(2,"EdbDataProc::FillCouplesTree","fill couples tree...");
 
  EdbPatCouple *patc=0;
  float xv = al->X();
  float yv = al->Y();
 
  int pid1,pid2;
  EdbSegCouple *cp=0;
  EdbSegP      *s1=0;
  EdbSegP      *s2=0;
  EdbSegP      *s=0;
 
  tree->SetBranchAddress("pid1",&pid1);
  tree->SetBranchAddress("pid2",&pid2);
  tree->SetBranchAddress("xv"  ,&xv);
  tree->SetBranchAddress("yv"  ,&yv);
  tree->SetBranchAddress("cp"  ,&cp);
  tree->SetBranchAddress("s1." ,&s1);
  tree->SetBranchAddress("s2." ,&s2);
  tree->SetBranchAddress("s."  ,&s );
 
  if(fillraw) {
    // **** fill tree with raw segments ****
    EdbPattern *pat=0;
    int nic;
    int nip=al->Npatterns();
    for( int ip=0; ip<nip; ip++ ) {
      pat  = al->GetPattern(ip);
      pid1 = pat->ID();
      pid2 = -1;
      nic=pat->N();
      for( int ic=0; ic<nic; ic++ ) {
        s1 = pat->GetSegment(ic);
        tree->Fill();
      }
    }
  }
 
  // **** fill tree with found couples ****
 
  s = new EdbSegP();
 
  int nip=al->Ncouples();
  for( int ip=0; ip<nip; ip++ ) {
    patc = al->GetCouple(ip);
    pid1 = patc->Pat1()->ID();
    pid2 = patc->Pat2()->ID();
 
    int nic=patc->Ncouples();
    for( int ic=0; ic<nic; ic++ ) {
      cp = patc->GetSegCouple(ic);
      s1 = patc->Pat1()->GetSegment(cp->ID1());
      s2 = patc->Pat2()->GetSegment(cp->ID2());
      s = cp->eS;
      tree->SetBranchAddress("s."  ,&s );
 
      s->SetID(tree->GetEntries());             // basetrack id will be the tree entry number
      EdbTraceBack::SetBaseTrackVid( *s, 0, 0, tree->GetEntries() );   //TODO: plate, piece if available
      tree->Fill();
    }
  }
 
}

FineAlignment()

void EdbDataProc::FineAlignment

(

int

doFine

)

---

{
  EdbPVRec    ali;
  InitVolume(&ali);
  ali.Link();
  Log(2,"EdbDataProc::FineAlignment","link ok");
  ali.FillTracksCell();
 
  EdbAffine2D aff;
 
  int fctr=0;
  int fcMin = 49;
 
  int i;
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  for( i=0; i<ali.Npatterns(); i++ ) {
    fctr = ali.FineCorrXY(i,aff,doFine);
  }
 
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  for( i=0; i<ali.Npatterns(); i++ ) {
    fctr = ali.FineCorrXY(i,aff,doFine);
  }
 
  if(fctr>fcMin) {
    for( i=0; i<ali.Npatterns(); i++) {
      ali.GetPattern(i)->GetKeep(aff);
      if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffPar(0,aff);
    }
  }
 
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  float dz=0;
  float z = ali.GetPattern(ali.Npatterns()-1)->Z();
  for( i=ali.Npatterns()-2; i>=0; i-- ) {
    fctr = ali.FineCorrZ(i,dz);
    if(fctr<=fcMin) break;
    z -= dz;
    Log(2,"EdbDataProc::FineAlignment","dz = %f  z = %f",dz,z);
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateZPar(0,z);
  }
 
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  for( i=0; i<ali.Npatterns(); i++ ) {
    fctr = ali.FineCorrTXTY(i,aff);
    if(fctr<=fcMin) break;
    if(gEDBDEBUGLEVEL>2) aff.Print();
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffTPar(0,aff);
  }
 
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  dz=0;
  z = ali.GetPattern(ali.Npatterns()-1)->Z();
  for( i=ali.Npatterns()-2; i>=0; i-- ) {
    fctr = ali.FineCorrZ(i,dz);
    if(fctr<=fcMin) break;
    z -= dz;
    Log(2,"EdbDataProc::FineAlignment","dz = %f  z = %f",dz,z);
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateZPar(0,z);
  }
 
  ali.SelectLongTracks(ali.Npatterns());
  ali.MakeSummaryTracks();
  float shr=0;
  for( i=0; i<ali.Npatterns(); i++ ) {
    fctr = ali.FineCorrShr(i,shr);
    if(fctr<=fcMin) break;
    eDataSet->GetPiece(i)->CorrectShrinkage(0,shr);
    if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateShrPar(0);
  }
 
}

FineAlignmentTracks()

void EdbDataProc::FineAlignmentTracks

(

)

---

{
  InitVolume(100);
  EdbPVRec *ali = PVR();
 
  EdbAffine2D aff,afft;
 
  int fctr=0,fctr0=10000, fcMin=50;
 
  for( int i=0; i<ali->Npatterns(); i++ ) {
    fctr = ali->FineCorrF(i,aff,afft);
    if(fctr<fctr0) fctr0=fctr;
    Log(2,"EdbDataProc::FineAlignmentTracks","fctr = %d",fctr);
  }
 
  if(fctr0>fcMin) {
    for( int i=0; i<ali->Npatterns(); i++) {
      ali->GetPattern(i)->GetKeep(aff);
      if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffPar(0,aff);
 
      //if(!NoUpdate())   eDataSet->GetPiece(i)->UpdateAffTPar(0,afft);
    }
  }
}

GetDataSet()

EdbDataSet * EdbDataProc::GetDataSet ( )

inline

{return eDataSet;}

GetPVR()

EdbPVRec * EdbDataProc::GetPVR ( ) const

inline

{return ePVR;}

GetVertexFromTree()

EdbVertex * EdbDataProc::GetVertexFromTree	(	EdbVertexRec &	vertexrec,
		const char *	fname,
		const int	vertexID
	)

{
  //reading vertex tree, getting only the vertex I need
  map<int,EdbTrackP*>emptymap;
  ReadVertexTree(vertexrec, fname, Form("vID==%i",vertexID),emptymap); //vertex is added to EdbPVRec
 
  EdbVertex *myvertex;
  EdbPVRec *ali = vertexrec.ePVR;
  //looking for myvertex (when called more than once, ali keeps adding vertices)
  for (int ivtx = 0; ivtx<ali->eVTX->GetEntries();ivtx++){
    EdbVertex *v1 = (EdbVertex*) ali->eVTX->At(ivtx);
    if (v1->ID()==vertexID) myvertex = v1;
  }
  
  return myvertex;
 
 }

InitVolume() [1/2]

int EdbDataProc::InitVolume	(	EdbPVRec *	ali,
		int	datatype = 0,
		TIndex2 *	trseg = 0
	)

---

TO CHECK!!!

{
  if (!eDataSet) Log(2,"EdbDataProc::InitVolume","Warning. eDataSet is NULL"); 
  
  EdbScanCond *cond = eDataSet->GetPiece(0)->GetCond(0);
  ali->SetScanCond( cond );
  
  EdbDataPiece *piece;
  int npieces = eDataSet->N();
  Log(2,"EdbDataProc::InitVolume","npieces = %d",npieces);
  if(!npieces) return 0;
 
  //TTree *cptree=0;
  int i;
 
  EdbPattern *pat=0;
  EdbPattern *p1=0;
  EdbPattern *p2=0;
 
  for(i=0; i<npieces; i++ ) {
    piece = eDataSet->GetPiece(i);
    if(!piece->InitCouplesTree("READ"))  Log(1," EdbDataProc::InitVolume","no tree %d",i);
 
    pat = new EdbPattern( 0.,0., piece->GetLayer(0)->Z(),100 );
    pat->SetPID(i);
    //pat->SetSegmentsPID();  //TO CHECK!!
    if(datatype==10) {
      p1 = new EdbPattern( 0.,0., piece->GetLayer(0)->Z() + piece->GetLayer(1)->Z() );
      p2 = new EdbPattern( 0.,0., piece->GetLayer(0)->Z() + piece->GetLayer(2)->Z() );
      p1->SetPID(i);
      p2->SetPID(i);
    }
 
    piece->GetCPData_new( pat,p1,p2,trseg );
 
    ali->AddPattern( pat );
    if(datatype==10) {
      ali->AddPattern( p1 );
      ali->AddPattern( p2 );
    }
    piece->CloseCPData();
  }
 
  float x0 = eDataSet->GetPiece(npieces-1)->GetLayer(0)->X();
  float y0 = eDataSet->GetPiece(npieces-1)->GetLayer(0)->Y();
  //ali->Centralize();
  ali->Centralize(x0,y0);
 
  int npat = ali->Npatterns();
  for(i=0; i<npat; i++ ) {
    pat = ali->GetPattern(i);
    pat->SetSegmentsZ();   
    pat->Transform(    eDataSet->GetPiece(pat->PID())->GetLayer(0)->GetAffineXY()   );
    pat->TransformA(   eDataSet->GetPiece(pat->PID())->GetLayer(0)->GetAffineTXTY() );
    pat->TransformShr( eDataSet->GetPiece(pat->PID())->GetLayer(0)->Shr()  );
  }
  ali->SetPatternsID();
 
  ali->SetSegmentsErrors();
 
  ali->SetCouplesAll();
  ali->SetChi2Max(cond->Chi2PMax());
  ali->SetOffsetsMax(cond->OffX(),cond->OffY());
  return npieces;
}

InitVolume() [2/2]

int EdbDataProc::InitVolume	(	int	datatype = 0,
		const char *	rcut = "1"
	)

---

datatype: 0 - couples basetrack data
10 - couples full data
100 - tracks only
1000 - tracks and basetracks

{
 
  if(!ePVR) ePVR = new EdbPVRec();
  if(datatype==100)  return InitVolumeTracks(ePVR, rcut);
 
  if(datatype==1000) {
    InitVolumeTracks(ePVR, rcut);
    TIndex2 *trseg=MakeTracksSegmentsList( *ePVR );
    delete ePVR;
    ePVR = new EdbPVRec();
    int n = InitVolume(ePVR, datatype, trseg);
    n += InitVolumeTracks(ePVR, rcut);
    return n;
  }
 
  return InitVolume(ePVR, datatype);
}

InitVolumeRaw()

int EdbDataProc::InitVolumeRaw

(

EdbPVRec *

ali

)

---

{
  EdbScanCond *cond = eDataSet->GetPiece(0)->GetCond(1);
  ali->SetScanCond( cond );
  
  EdbDataPiece *piece;
  int npieces = eDataSet->N();
  Log(2,"EdbDataProc::InitVolumeRaw","npieces = %d",npieces);
  if(!npieces) return 0;
 
  for(int i=0; i<npieces; i++ ) {
    piece = eDataSet->GetPiece(i);
    piece->GetRawData(ali);
  }
 
  float x0 = eDataSet->GetPiece(npieces-1)->GetLayer(0)->X();
  float y0 = eDataSet->GetPiece(npieces-1)->GetLayer(0)->Y();
  ali->Centralize(x0,y0);
 
  for(int j=0; j<npieces; j++ ) {
    for(int ip=0; ip<2; ip++ ) {
      int i = 2*j+ip;
      ali->GetPattern(i)->SetSegmentsZ();
      ali->GetPattern(i)->Transform(    eDataSet->GetPiece(j)->GetLayer(0)->GetAffineXY()   );
      ali->GetPattern(i)->TransformA(   eDataSet->GetPiece(j)->GetLayer(0)->GetAffineTXTY() );
      ali->GetPattern(i)->TransformShr( eDataSet->GetPiece(j)->GetLayer(0)->Shr()  );
    }
  }
  ali->SetPatternsID();
  ali->SetSegmentsErrors();
 
  ali->SetCouplesAll();
  ali->SetChi2Max(cond->Chi2PMax());
  ali->SetOffsetsMax(cond->OffX(),cond->OffY());
  return npieces;
}

InitVolumeTracks()

int EdbDataProc::InitVolumeTracks	(	EdbPVRec *	ali,
		const char *	rcut
	)

---

{
  EdbScanCond *cond=0;
  if(eDataSet) {
    cond= eDataSet->GetPiece(0)->GetCond(0);
    ali->SetScanCond( cond );
  
    EdbDataPiece *piece;
    int npieces = eDataSet->N();
    Log(2,"EdbDataProc::InitVolumeTracks","npieces = %d",npieces);
    if(!npieces) return 0;
 
    if( ali->Npatterns() == 0) {                    // case of segments already readed
      EdbAffine2D *a=0;
      EdbPattern  *pat=0;
      for(int i=0; i<npieces; i++ ) {
    piece = eDataSet->GetPiece(i);
    pat = new EdbPattern( 0.,0., piece->GetLayer(0)->Z(),100 );
    a = piece->GetLayer(0)->GetAffineXY();
    pat->SetKeep( a->A11(),a->A12(),a->A21(),a->A22(), a->B1(),a->B2() );
    pat->SetPID(i);
    ali->AddPattern( pat );
      }
      ali->SetPatternsID();
    }
  }
  else {
    cond = new EdbScanCond();
    ali->SetScanCond( cond );
  }
 
  ReadTracksTree( *ali, "linked_tracks.root", rcut );
 
  ali->SetSegmentsTracks();  // reset the id of tracks????? to check that it is really necessary!!
  ali->SetSegmentsErrors();
  ali->SetCouplesAll();
  ali->SetChi2Max(cond->Chi2PMax());
  ali->SetOffsetsMax(cond->OffX(),cond->OffY());
  return ali->Npatterns();
}

Link() [1/2]

int EdbDataProc::Link

(

)

{
  if (!eDataSet) return 0;
  EdbDataPiece *piece;
  Int_t np = eDataSet->N();
  for (Int_t i = 0; i < np; i++) {
    piece = eDataSet->GetPiece(i);
    if (piece->TakePiecePar() >= 0) {
      if (gEDBDEBUGLEVEL > 2) piece->Print();
      piece->WriteCuts();
      if (piece->Flag() == 1) Link(*piece);
    }
    piece->CloseRun();
  }
  return np;
}

Link() [2/2]

int EdbDataProc::Link

(

EdbDataPiece &

piece

)

!! 0

{
  EdbPVRec  *ali;
 
  //const char *file_name=piece.GetNameCP();
  //TTree *cptree=EdbDataPiece::InitCouplesTree(file_name,"RECREATE");
  //if(!cptree) return 0;
 
  if (!piece.InitCouplesTree("RECREATE")) return 0;
 
  EdbScanCond *cond = piece.GetCond(1);
  int    ntot=0, nareas=0;
  float  shr1=1,shr2=1;
  double shrtot1=0., shrtot2=0.;
  int    nshr=0,nshrtot=0;
 
  for( int irun=0; irun<piece.Nruns(); irun++ ) {
    nareas = piece.MakeLinkListCoord(irun);
    //nareas = piece.MakeLinkListArea(irun);
    Log(2,"EdbDataProc::Link","%d areas mapped",nareas);
    if (nareas <= 0) continue;
 
    for (Int_t i = 0; i < nareas; i++ ) {
 
      ali      = new EdbPVRec();
      ali->SetScanCond( cond );
 
      piece.GetAreaData(ali,i,1);
      piece.GetAreaData(ali,i,2);
 
      //ali->SetSegmentsErrors();
 
      ali->SetCouplesAll();
      ali->SetChi2Max(cond->Chi2PMax());
 
      for(int ic=0; ic<ali->Ncouples(); ic++) 
    ali->GetCouple(ic)->SetCHI2mode(cond->Chi2Mode());
 
      ali->Link();
 
      if( ShrinkCorr() ) {
    shr1 = 1;
    shr2 = 1;
    nshr = CheckShrinkage( ali,0, shr1, shr2 );
    if(nshr) {
      nshrtot += nshr;
      shrtot1 += nshr*shr1;
      shrtot2 += nshr*shr2;
    }
      }
 
      FillCouplesTree( piece.eCouplesTree, ali,piece.GetOUTPUT());  
 
      delete ali;
    }
    ntot+=nareas;
    if(nshrtot>3) {
      shrtot1 = shrtot1/nshrtot;
      shrtot2 = shrtot2/nshrtot;
      if(nshrtot<20) Log(1,"EdbDataProc::Link","WARNING: unreliable shrinkage correction - low statistics");
      Log(2,"EdbDataProc::Link","Shrinkage correction(%d): %f %f", nshrtot, (float)shrtot1,(float)shrtot2);
      piece.CorrectShrinkage( 1, (float)shrtot1 );
      piece.CorrectShrinkage( 2, (float)shrtot2 );
      if(!NoUpdate())   piece.UpdateShrPar(1);
      if(!NoUpdate())   piece.UpdateShrPar(2);
    }
  }
  int ncp = piece.eCouplesTree->GetEntries();
  piece.CloseCPData();
  return ncp;
}

LinkRawTracks()

void EdbDataProc::LinkRawTracks

(

int

alg = 0

)

---

{
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolumeRaw(ali);
  ali->Link();
  Log(2,"EdbDataProc::LinkRawTracks","link ok");
 
//    if(alg==1) {
//      ali->MakeHoles();
//      ali->Link();
//      printf("link ok\n");
//    }
 
  ali->FillTracksCell();
 
  TTree *cptree=EdbDataPiece::InitCouplesTree("linked_couples.root","RECREATE");
  FillCouplesTree(cptree, ali,0);
  CloseCouplesTree(cptree);
 
  MakeTracksTree(ali);
}

LinkTracks()

void EdbDataProc::LinkTracks	(	int	alg = 0,
		float	p = -1.
	)

---

{
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolume(ali);
 
  Log(2,"EdbDataProc::LinkTracks","tracking:  alg = %d \t p= %f", alg, p);
 
  ali->SetCouplesPeriodic(0,1);
  int ntr=0;
 
  ntr = LinkTracksWithFlag( ali, p, 0.05, 2, 3, 0 );
  //ntr = LinkTracksWithFlag( ali, p/4., 0.01, 2, 3, 1 );
  //if(merge>0) ali->MergeTracks(merge);
 
  if(alg==-1) { // this is a test option produce couples tree from the first 2 linked patterns
    TTree *cptree=EdbDataPiece::InitCouplesTree("linked_couples.root","RECREATE");
    FillCouplesTree(cptree, ali,0);
    CloseCouplesTree(cptree);
  }
 
  float mass=0.139;                    //TODO!
  ali->FitTracks( p, mass );         // is important to call it before MakeTracksTree!
 
  MakeTracksTree(ali);
  
  // After the EdbPVRec is created and filled, set it also as the EdbDataProc own object.
  SetPVR(ali);
  if (gEDBDEBUGLEVEL>2) {
    ali->Print();
  }
  Log(2,"EdbDataProc::LinkTracks","EdbDataProc::LinkTracks...Done.");
  return;
}

LinkTracksC()

void EdbDataProc::LinkTracksC	(	int	alg = 0,
		float	p = -1.
	)

---

special tracking procedure for carbonium dataset

{
 
  EdbPVRec    *ali  = new EdbPVRec();
  InitVolume(ali);
 
  Log(2,"EdbDataProc::LinkTracksC","carbonium tracking:  alg = %d \t p= %f", alg, p);
 
  ali->SetCouplesPeriodic(0,1);
  LinkTracksWithFlag( ali, p, 0.01, 2, 3, 0 );
 
  ali->SetCouplesPeriodic(0,3);
  LinkTracksWithFlag( ali, p, 0.01, 2, 6, 1 );
 
  ali->SetCouplesPeriodic(1,3);
  LinkTracksWithFlag( ali, p, 0.01, 2, 6, 2 );
 
  ali->SetCouplesPeriodic(2,3);
  LinkTracksWithFlag( ali, p, 0.01, 2, 6, 3 );
 
  float mass=0.139;                    //TODO!
  ali->FitTracks( p, mass );         // is important to call it before MakeTracksTree!
 
  MakeTracksTree(ali);
}

LinkTracksWithFlag()

int EdbDataProc::LinkTracksWithFlag ( EdbPVRec *  ali, float  p, float  probmin, int  nsegmin, int  maxgap, int  flag, float  mass = 0.1396  )

static

---

p - momentum expected for tracks
probmin - minimum probability to accept tracks
nsegmin - min number of segments/track
maxgap - max gap permitted for propagation
flag - to be assigned for the tracks found in this pass

{
 
  int ntr0 = ali->Ntracks();
  ali->Link();
  ali->FillTracksCell();
  ali->MakeTracks(nsegmin,flag);
 
  int       noProp=0;
  if(p<0)   noProp=1;
  //if(p<0.01) p=4.;  //TODO: this protection should be out of this function?
 
  int ntr = ali->Ntracks();
  float X0 =  ali->GetScanCond()->RadX0();
  EdbTrackP *tr=0;
  for(int itr=ntr0; itr<ntr; itr++) {
    tr = ali->GetTrack(itr);
    tr->ClearF();
    if(p<0) tr->SetP(1.);
    else    tr->SetP(p);
    tr->SetM(mass);
    tr->FitTrackKFS(false,X0);
  }
 
  if(!noProp)
    if( ali->Npatterns()>2 ) {
      ali->FillCell(50,50,0.015,0.015);
      for(int i=0; i<10; i++) 
    if( ali->PropagateTracks(ali->Npatterns()-1,2, probmin, maxgap ) <1) break;
    }
 
  /*
  ntr = ali->Ntracks();
  for(int i=0; i<ntr; i++) {
    tr = ali->GetTrack(i);
    if(tr->Flag()<0) continue;
    if(tr->N()<5)  tr->SetP(0.5); 
    else           tr->SetP(tr->P_MS());
  }
 
  for(int i=0; i<10; i++) 
    if( ali->PropagateTracks(ali->Npatterns()-1,2, probmin, maxgap ) <1) break;
  */
 
  return ntr-ntr0;
}

MakeTracksSegmentsList()

TIndex2 * EdbDataProc::MakeTracksSegmentsList

(

EdbPVRec &

ali

)

{
  TIndex2 *itracks = 0;
  int nsegtot=0;
  EdbTrackP *tr=0;
  int ntr = ali.eTracks->GetEntries();
 
  if(!ntr) return 0;
  Double_t *w    = new Double_t[500000];
  itracks = new TIndex2();
 
  for(int i=0; i<ntr; i++) {
    tr = (EdbTrackP*)(ali.eTracks->At(i));
    if(!tr) continue;
    EdbSegP *s=0;
    int nseg=tr->N();
    for(int j=0; j<nseg; j++) {
      s = tr->GetSegment(j);
      if(!s) continue;
      w[nsegtot] = TIndex2::BuildValue( s->Vid(0), s->Vid(1) );
      nsegtot++;
    }
  }
  itracks->Set(0);
  itracks->BuildIndex(nsegtot,w);
  delete[] w;
  return itracks;
}

MakeTracksTree() [1/2]

int EdbDataProc::MakeTracksTree ( EdbPVRec *  ali = 0, const char *  file = "linked_tracks.root"  )

static

{
  if(!ali) return 0;
  TObjArray *trarr = ali->eTracks;
  if(!trarr) return 0;
  float xv=ali->X();
  float yv=ali->Y();
  
  return MakeTracksTree(*trarr,xv,yv, file);
}

MakeTracksTree() [2/2]

int EdbDataProc::MakeTracksTree ( TObjArray &  tracks, float  xv = 0, float  yv = 0, const char *  file = "linked_tracks.root"  )

static

{
  Log(2,"EdbDataProc::MakeTracksTree","write tracks into %s ... ",file);
  TFile fil(file,"RECREATE");
  TTree *tracks= new TTree("tracks","tracks");
 
  EdbTrackP    *track = new EdbTrackP(8);   // why is this track initialised with 8 segments by default ???
  EdbSegP      *tr = new EdbSegP();
  TClonesArray *segments  = new TClonesArray("EdbSegP");
  TClonesArray *segmentsf = new TClonesArray("EdbSegP");
 
  int   nseg,trid,npl,n0;
  float w=0.;
 
  tracks->Branch("trid",&trid,"trid/I");
  tracks->Branch("nseg",&nseg,"nseg/I");
  tracks->Branch("npl",&npl,"npl/I");
  tracks->Branch("n0",&n0,"n0/I");
  tracks->Branch("xv",&xv,"xv/F");
  tracks->Branch("yv",&yv,"yv/F");
  tracks->Branch("w",&w,"w/F");
  tracks->Branch("t.","EdbSegP",&tr,32000,99); // track data is now filled after using Copy method
  //tracks->Branch("t.","EdbSegP",&track,32000,99);
  tracks->Branch("s", &segments);
  tracks->Branch("sf",&segmentsf);
 
  int ntr = trarr.GetEntriesFast();
  
  for(int itr=0; itr<ntr; itr++) {
    
    track = (EdbTrackP*)(trarr.At(itr));
    tr->Copy(*track); //Using Copy method from EdbSegP avoids the cast
    tr->ForceCOV(track->COV());
    trid = track->ID();
    nseg = track->N();
    npl  = track->Npl();
    n0   = track->N0();
    
    segments->Clear("C");
    segmentsf->Clear("C");
    nseg = track->N();
    w    = track->Wgrains();
    EdbSegP *s=0,*sf=0;
    for(int is=0; is<nseg; is++) {
      s = track->GetSegment(is);
      if(s) new((*segments)[is])  EdbSegP( *s );
      sf = track->GetSegmentF(is);
      if(sf) new((*segmentsf)[is])  EdbSegP( *sf );
    }
 
    tr->SetVid( 0, tracks->GetEntries() );  // put track counter in t.eVid[1]
    tracks->Fill();
    // track->Clear(); // if this Clear() is written at this stage, the corresponding track in
    // the EdbPattern Volume (ali, eTracks) will also be cleared! This behavior is not
    // wanted, is it? So we do comment this out here ... (Frank, 09 26 2016)
  }
 
  tracks->Write();
  fil.Close();
  Log(2,"EdbDataProc::MakeTracksTree","%d tracks are written",ntr);
  return ntr; 
}

MakeVertexTree()

int EdbDataProc::MakeVertexTree ( TObjArray &  vtxarr, const char *  file  )

static

{
  Log(2,"EdbDataProc::MakeVertexTree","write vertices into %s ... ",file);
  TFile fil(file,"RECREATE");
  TTree *vtx= new TTree("vtx","Reconstructed vertices in emulion");
  
  EdbSegP      *tr = new EdbSegP();
  EdbTrackP *track = NULL;
 
  TClonesArray *tracks = new TClonesArray("EdbSegP"); //tracks are saved as EdbSegP, like in MakeTracksTree
  TClonesArray *segments  = new TClonesArray("EdbSegP");
  TClonesArray *segmentsf = new TClonesArray("EdbSegP");
 
  //list of variables
 
  Float_t vx, vy, vz; //true reconstructed vertex position
  Float_t meanvx, meanvy, meanvz; //track connection point
  Int_t vID = 0;
  Float_t maxaperture;
  Float_t probability;
  Int_t n;
  Int_t flag; //we need also the flag to check vertex type
  const Int_t maxdim = 1000; //maximum number of tracks associated to a vertex that can be saved in the tree
  //big arrays for containers of track variables
  Int_t TrackID[maxdim];
  Int_t nholes[maxdim];
  Int_t maxgap[maxdim];
  Int_t nseg[maxdim];
  Int_t npl[maxdim];
  Float_t TX[maxdim];
  Float_t TY[maxdim];
  Float_t impactparameter[maxdim];
  Int_t incoming[maxdim]; //coming from the vertex
  //MC information
  Int_t MCEventID[maxdim];
  Int_t MCTrackID[maxdim];
  Int_t MCTrackPdgCode[maxdim];
  Int_t MCMotherID[maxdim];
 
  //list of branches
  vtx->Branch("vID",&vID,"vID/I");
  vtx->Branch("flag",&flag,"flag/I");
  vtx->Branch("vx",&vx,"vx/F");
  vtx->Branch("vy",&vy,"vy/F");
  vtx->Branch("vz",&vz,"vz/F");
  vtx->Branch("meanvx",&meanvx,"meanvx/F");
  vtx->Branch("meanvy",&meanvy,"meanvy/F");
  vtx->Branch("meanvz",&meanvz,"meanvz/F");
  vtx->Branch("maxaperture",&maxaperture,"maxaperture/F");
  //vtx->Branch("maxrmsthetaspace",&maxrmsthetaspace,"maxrmsthetaspace/F");
  vtx->Branch("probability",&probability,"probability/F");
  vtx->Branch("n",&n,"n/I");
  vtx->Branch("t.",&tracks); // track data is now filled after using Copy method
  //tracks->Branch("t.","EdbSegP",&track,32000,99);
  vtx->Branch("s", &segments);
  vtx->Branch("sf",&segmentsf);
  //track variables (they are array with the number of tracks as size)
  vtx->Branch("TrackID",&TrackID,"TrackID[n]/I");
  vtx->Branch("nseg",&nseg,"nseg[n]/I");
  vtx->Branch("npl",&npl,"npl[n]/I");
  vtx->Branch("nholes",&nholes,"nholes[n]/I"); //even more variables in this tree
  vtx->Branch("maxgap",&maxgap,"maxgap[n]/I"); 
  vtx->Branch("incoming",&incoming,"incoming[n]/I");
  vtx->Branch("impactparameter",&impactparameter,"impactparameter[n]/F");
  //inserting MCtrue information
  vtx->Branch("MCEventID", &MCEventID, "MCEventID[n]/I");
  vtx->Branch("MCTrackID",&MCTrackID,"MCTrackID[n]/I");
  vtx->Branch("MCTrackPdgCode",&MCTrackPdgCode,"MCTrackPdgCode[n]/I");
  vtx->Branch("MCMotherID",&MCMotherID,"MCMotherID[n]/I");
 
  int nvtx = vtxarr.GetEntriesFast();
  //START TREE FILLING
  for(int ivtx=0; ivtx<nvtx; ivtx++) {
    int itotalseg = 0;
    tracks->Clear("C");
    segments->Clear("C");
    segmentsf->Clear("C");
 
    EdbVertex* vertex = (EdbVertex*) (vtxarr.At(ivtx));
    if(vertex->Flag()<0) continue; //saving only 'true' vertices in the tree file and in the object
    vx=vertex->VX();
    vy=vertex->VY();
    vz=vertex->VZ();
    meanvx=vertex->X();
    meanvy=vertex->Y();
    meanvz=vertex->Z();
    n=vertex->N();
    maxaperture = vertex->MaxAperture();
    probability = vertex->V()->prob();  
    flag = vertex->Flag();
    //adding vertex to list to be saved
    //loop on tracks //now it can be done offline (again)
    for (int itrk = 0; itrk < n; itrk++){
     //getting track and track variables to fill branches
     track = vertex->GetTrack(itrk);
     tr->Copy(*track);     
//     tr->ForceCOV(track->COV());
     if(tr) new((*tracks)[itrk])  EdbSegP( *tr ); //adding track to trackclonesarray
     ((EdbSegP*) tracks->At(itrk))->ForceCOV(track->COV()); //default copy does NOT save COV corrrectly
     TrackID[itrk] = track->Track(); //the eTrack attribute of EdbSegP now allows association to original root tree
     nseg[itrk] = track->N();
     npl[itrk]  = track->Npl();
     Int_t zpos = vertex->GetVTa(itrk)->Zpos();
     incoming[itrk] = zpos;
     nholes[itrk] = track->N0();
     maxgap[itrk] = track->CheckMaxGap();
     impactparameter[itrk] = vertex->GetVTa(itrk)->Imp();
     //Storing MCTrue information (of course for real data these values have no sense)
     MCEventID[itrk] = track->MCEvt();
     MCTrackID[itrk] = track->MCTrack();
     if(MCEventID[itrk]>-1){ //vertices from MC simulation, used Vid[0] and Aid[0] to store these information
      MCTrackPdgCode[itrk] = track->Vid(0);
      MCMotherID[itrk] = track->Aid(0); //used to store MotherID information
     }
     else { //vertices from data, set default values 
      MCTrackPdgCode[itrk] = -999;
      MCMotherID[itrk] = -999;
     }
    //w    = track->Wgrains();
     EdbSegP *s=0,*sf=0;
    //loop on segments
     for(int is=0; is<nseg[itrk]; is++) {
      s = track->GetSegment(is);     
      if(s) new((*segments)[itotalseg])  EdbSegP( *s );
      ((EdbSegP*) segments->At(itotalseg))->ForceCOV(s->COV()); //default copy does NOT save COV corrrectly
      sf = track->GetSegmentF(is);
      if(sf) new((*segmentsf)[itotalseg])  EdbSegP( *sf );
      ((EdbSegP*) segmentsf->At(itotalseg))->ForceCOV(sf->COV()); //default copy does NOT save COV corrrectly
      itotalseg++;
      }
 
   // track->SetVid( 0, tracks->GetEntries() );  // put track counter in t.eVid[1]
     }
    vtx->Fill();
    vID++; //now the number of elements in the tree and vertices is the same. vID starts from 0
    }
  vtx->Write();
  fil.Close();
  Log(2,"EdbDataProc::MakeVertexTree","%d vertices are written",nvtx);
  return nvtx; 
}

NoUpdate()

int EdbDataProc::NoUpdate ( ) const

inline

{ return eNoUpdate; }

Process()

int EdbDataProc::Process ( )

inline

{ return Link(); }  // to be removed

PVR()

EdbPVRec * EdbDataProc::PVR ( ) const

inline

{return ePVR;}

ReadTracksTree()

int EdbDataProc::ReadTracksTree ( EdbPVRec &  ali, const char *  fname = "linked_tracks.root", const char *  rcut = "t.eFlag>-1&&nseg>2&&t.eProb>.01"  )

static

{
  TFile f(fname);
  if(f.IsZombie()) { Log(1,"EdbDataProc::ReadTracksTree","Error open file %s", fname);  return 0; }
  
  TTree *tracks = (TTree*)f.Get("tracks");
  Int_t   trid=0;
  Int_t   nseg=0;
  Int_t   npl=0;
  Int_t   n0=0;
  Float_t xv=0.;
  Float_t yv=0.;
 
  int nentr = (int)(tracks->GetEntries());
  TCut cut = rcut;
  tracks->Draw(">>lst", cut );
  TEventList *lst = (TEventList*)gDirectory->GetList()->FindObject("lst");
  int nlst =lst->GetN();
  if(cut) Log(2,"EdbDataProc::ReadTracksTree","select %d of %d tracks by cut %s",nlst, nentr, cut.GetTitle() );
 
  TClonesArray *seg  = new TClonesArray("EdbSegP", 60);
  TClonesArray *segf = new TClonesArray("EdbSegP", 60);
  EdbSegP *trk=0;
  EdbSegP *s1=0;
  EdbSegP *s1f=0;
 
  tracks->SetBranchAddress("trid", &trid);
  tracks->SetBranchAddress("nseg", &nseg);
  tracks->SetBranchAddress("npl", &npl);
  tracks->SetBranchAddress("n0", &n0);
  tracks->SetBranchAddress("xv", &xv);
  tracks->SetBranchAddress("yv", &yv);
  tracks->SetBranchAddress("sf", &segf);
  tracks->SetBranchAddress("s",  &seg);
  tracks->SetBranchAddress("t.", &trk);
 
  EdbPattern *pat=0;
  int entr=0;
  for (int j=0; j<nlst; j++){
    entr = lst->GetEntry(j);
    tracks->GetEntry(entr);
 
    EdbTrackP *tr1 = new EdbTrackP();
    ((EdbSegP*)tr1)->Copy(*trk);
    tr1->ForceCOV(trk->COV());
    tr1->SetM(0.139);                 //TODO
 
    for(int i=0; i<nseg; i++) {
      s1 = (EdbSegP*)(seg->At(i));
      s1f = (EdbSegP*)(segf->At(i));
      pat = ali.GetPattern( s1->PID() );
      if(!pat) { 
    Log(1,"EdbDataProc::ReadTracksTree","WARNING: no pattern with pid %d: creating new one!",s1->PID()); 
    pat = new EdbPattern( 0., 0., s1->Z() );
    pat->SetID(s1->PID());
    pat->SetScanID(s1->ScanID());
    ali.AddPatternAt(pat,s1->PID());
      }
      //This way it should be more clear which segment instance is added to the track
      EdbSegP *segtobeadded = pat->AddSegment(*s1);
      EdbSegP *segftobeadded = new EdbSegP(*s1f);
      segtobeadded->ForceCOV(s1->COV());
      segftobeadded->ForceCOV(s1f->COV());
 
      tr1->AddSegment( segtobeadded );
      tr1->AddSegmentF( segftobeadded );
    }
    tr1->SetSegmentsTrack(tr1->ID());
    tr1->SetCounters();
    //tr1->FitTrackKFS(true);
    tr1->SetTrack(trid); //providing trackid to eTrack so it will not be lost when trackID resets
    ali.AddTrack(tr1);
  }
 
  Log(2,"EdbDataProc::ReadTracksTree","%d tracks are read",nlst);
  return nlst;
}

ReadVertexTree() [1/2]

int EdbDataProc::ReadVertexTree ( EdbVertexRec &  vertexrec, const char *  fname, const char *  rcut, map< int, EdbTrackP * >  trackID_map  )

static

{
  TFile f(fname);
  if(f.IsZombie()) { Log(1,"EdbDataProc::ReadVertexTree","Error open file %s", fname);  return 0; }
  
  TTree *vtx = (TTree*)f.Get("vtx");
 
  //vertex variables
  Float_t vx, vy, vz;
  Int_t vID = 0;
  Int_t n = 0;
  Int_t flag = 0;
  //track variables
 
  const Int_t maxdim = 1000; //maximum number of tracks associated to a vertex that can be saved in the tree
  //big arrays for containers of track variables
  Int_t TrackID[maxdim];
  Int_t nholes[maxdim];
  Int_t maxgap[maxdim];
  Int_t nseg[maxdim];
  Float_t TX[maxdim];
  Float_t TY[maxdim];
  Float_t impactparameter[maxdim];
  Int_t incoming[maxdim]; //coming from the vertex
  //MC information
  Int_t MCEventID[maxdim];
  Int_t MCTrackID[maxdim];
  Int_t MCMotherID[maxdim];
/*
  Int_t   trid=0;
  Int_t   nseg=0;
  Int_t   npl=0;
  Int_t   n0=0;
  Float_t xv=0.;
  Float_t yv=0.;*/
 
  int nentr = (int)(vtx->GetEntries());
  TCut cut = rcut;
  vtx->Draw(">>lst", cut );
  TEventList *lst = (TEventList*)gDirectory->GetList()->FindObject("lst");
  int nlst =lst->GetN();
  if(cut) Log(2,"EdbDataProc::ReadVtxTree","select %d of %d vertices by cut %s",nlst, nentr, cut.GetTitle() );
 
  TClonesArray *tracks = new TClonesArray("EdbSegP");
  TClonesArray *seg  = new TClonesArray("EdbSegP", 60);
  TClonesArray *segf = new TClonesArray("EdbSegP", 60);
  
  EdbSegP *trk=0;
  EdbSegP *s1=0; 
  EdbSegP *s1f=0;
  //vertex variables
  vtx->SetBranchAddress("vID",&vID);
  vtx->SetBranchAddress("flag",&flag);
  vtx->SetBranchAddress("vx",&vx);
  vtx->SetBranchAddress("vy",&vy);
  vtx->SetBranchAddress("vz",&vz);
  vtx->SetBranchAddress("n",&n);
  //arrays
  vtx->SetBranchAddress("nseg",&nseg);
  vtx->SetBranchAddress("TrackID",&TrackID);
  vtx->SetBranchAddress("incoming",&incoming);
  //TClones arrays
  vtx->SetBranchAddress("sf", &segf);
  vtx->SetBranchAddress("s",  &seg);
  vtx->SetBranchAddress("t.", &tracks);
 
  EdbPattern *pat=0;
  int entr=0;
  //now each tree entry is a vertex
 
  TObjArray * vertextracks = new TObjArray(maxdim);
  float expz = 0.;
 
  EdbPVRec *ali = vertexrec.ePVR; //I get EdbPVR from EdbVertexRec from now on
  for (int j=0; j<nlst; j++){
    //RESET COUNTERS!
    int itotalseg = 0; //counter for all the segments
    vertextracks->Clear("C");
 
    entr = lst->GetEntry(j);
    vtx->GetEntry(entr);
    EdbVertex *v1 = new EdbVertex();
    //loop on all tracks associated to the vertex
    for (int itrk = 0; itrk < n; itrk++){
     EdbTrackP *tr1;
     if(!trackID_map[TrackID[itrk]]){
      tr1= new EdbTrackP();
      trk = (EdbSegP*) tracks->At(itrk); //getting track itrk
      ((EdbSegP*)tr1)->Copy(*trk);
      tr1->ForceCOV(trk->COV());
      tr1->SetM(0.139);                 //TODO
 
      for(int i=0; i<nseg[itrk]; i++) {
       s1 = (EdbSegP*)(seg->At(itotalseg));
       s1f = (EdbSegP*)(segf->At(itotalseg));
       pat = ali->GetPattern( s1->PID() );
       if(!pat) { 
          Log(1,"EdbDataProc::ReadVertexTree","WARNING: no pattern with pid %d: creating new one!",s1->PID()); 
          pat = new EdbPattern( 0., 0., s1->Z() );
          pat->SetID(s1->PID());
          pat->SetScanID(s1->ScanID());
          ali->AddPatternAt(pat,s1->PID());
      }
 
        //This way it should be more clear which segment instance is added to the track
        EdbSegP *segtobeadded = pat->AddSegment(*s1);
        EdbSegP *segftobeadded = new EdbSegP(*s1f);
        segtobeadded->ForceCOV(s1->COV());
        segftobeadded->ForceCOV(s1f->COV());
 
        tr1->AddSegment(segtobeadded);
        tr1->AddSegmentF(segftobeadded);
        itotalseg++; //increasing counter of number of segments
      }
      tr1->SetSegmentsTrack(tr1->ID());
      tr1->SetCounters();
     //tr1->FitTrackKFS(true);
      tr1->SetTrack(TrackID[itrk]); //providing trackid to eTrack so it will not be lost when trackID resets
      ali->AddTrack(tr1);
      trackID_map[TrackID[itrk]] = tr1;
     }
     else{ //track already added to ali, add it only to vertex
      tr1 = trackID_map[TrackID[itrk]];
      itotalseg+= nseg[itrk];
     }
     vertextracks->Add(tr1);     
     //v1 = vertexrec->AddTrackToVertex(v1, tr1, incoming[itrk]);
    }
    //setting vertex parameters and saving vertex
    //v1 = vertexrec.Make1Vertex(*vertextracks,vz) (copied here instead, by putting zpos directly as for in vertexing);
    v1->SetXYZ( 0,0, vz );
    int ntr = vertextracks->GetEntries();
    for(int i=0; i<ntr; i++) {
      EdbTrackP *t = (EdbTrackP*)vertextracks->At(i);
      EdbVTA *vta = new EdbVTA(t,v1);
      vta->SetFlag(2);
      v1->AddVTA(vta);
      vta->SetZpos(incoming[i]);
      t->AddVTA(vta);
    }
    if( vertexrec.MakeV(*v1) )  vertexrec.AddVertex(v1);
    //else { SafeDelete(v); return 0; }                                // vertex is not valid
    //return v;
    v1->SetID(vID);
    v1->SetFlag(flag);
    /* make vertex manually, old version
    int ntr = vertextracks->GetEntries();
    v1->SetXYZ(vx,vy,vz);  
    for(int i=0; i<ntr; i++) {
     EdbTrackP *t = (EdbTrackP*)vertextracks->At(i);
     EdbVTA *vta = new EdbVTA(t,v1);
     vta->SetFlag(2);
     v1->AddVTA(vta);
     (t->Z() >= v1->Z())? vta->SetZpos(1) : vta->SetZpos(0);
     t->AddVTA(vta);
    } */
    
    ali->AddVertex(v1);
  }
 
  Log(2,"EdbDataProc::ReadVtxTree","%d vertices are read",nlst);
  return nlst;
}

ReadVertexTree() [2/2]

int EdbDataProc::ReadVertexTree ( EdbVertexRec &  vertexrec, const char *  fname, const char *  rcut, TObjArray *  builttracks = 0  )

static

                                                                                                                        {
 
  map<int,EdbTrackP*>emptymap;
  if (builttracks){
   int ntracks = builttracks->GetEntries();
  
   for (int itrk = 0; itrk < ntracks; itrk++){
        EdbTrackP * track = (EdbTrackP*) builttracks->At(itrk);
        emptymap[track->Track()] = track; //adding this trackID to map
    }
  }
 
  int nlst = ReadVertexTree( vertexrec, fname, rcut,emptymap);
  return nlst;
 }

SetNoUpdate()

void EdbDataProc::SetNoUpdate ( int  nu )

inline

{ eNoUpdate=nu; }

SetPVR()

void EdbDataProc::SetPVR ( EdbPVRec *  pvr )

inline

{ePVR=pvr;}

ShrinkCorr()

int EdbDataProc::ShrinkCorr ( )

inline

{return 1;}

Member Data Documentation

eDataSet

EdbDataSet* EdbDataProc::eDataSet

private

eNoUpdate

int EdbDataProc::eNoUpdate

private

ePVR

EdbPVRec* EdbDataProc::ePVR

private

---

The documentation for this class was generated from the following files:

• /home/antonio/fedra_doxygen/src/libEIO/EdbDataSet.h
• /home/antonio/fedra_doxygen/src/libEIO/EdbDataSet.cxx