[u/mrichter/AliRoot.git] / ITS / TPCtracks.C

struct GoodTrack {
  Int_t lab;
  Int_t code;
  Float_t px,py,pz;
  Float_t x,y,z;
  Float_t pxg,pyg,pzg,ptg;
  Bool_t flag;
};
Int_t good_tracks(GoodTrack *gt, Int_t max);


Int_t TPCtracks() {
   cerr<<"Doing comparison...\n";

  // Connect the Root Galice file containing Geometry, Kine and Hits


   if (gClassTable->GetID("AliRun") < 0) {
      gROOT->LoadMacro("loadlibs.C");
      loadlibs();
   }


   TFile *cf=TFile::Open("AliTPCclusters.root");
   if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
   AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
   if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }

// Load clusters  
   AliTPCClustersArray *ca=new AliTPCClustersArray;
   ca->Setup(digp);
   ca->SetClusterType("AliTPCcluster");
   ca->ConnectTree("Segment Tree");
   Int_t nentr=Int_t(ca->GetTree()->GetEntries());
   for (Int_t i=0; i<nentr; i++) {
       AliSegmentID *s=ca->LoadEntry(i);
   }

// Load tracks
   TFile *tf=TFile::Open("AliTPCtracks.root");
   if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
   TObjArray tarray(2000);
   TTree *tracktree=(TTree*)tf->Get("TreeT");
   TBranch *tbranch=tracktree->GetBranch("tracks");
   Int_t nentr=tracktree->GetEntries();
   for (i=0; i<nentr; i++) {
       AliTPCtrack *iotrack=new AliTPCtrack;
       tbranch->SetAddress(&iotrack);
       tracktree->GetEvent(i);
       iotrack->CookLabel(ca);
       tarray.AddLast(iotrack);
   }   

   cerr<<"Number of found tracks "<<nentr<<endl;
   tarray.Sort();

   TFile *pfile = new TFile("tpctracks.root","RECREATE"); 

   TTree tracktree1("TreeT","Tree with TPC tracks");
   tracktree1.Branch("tracks","AliTPCtrack",&iotrack,32000,0);

   for (i=0; i<nentr; i++) {
       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
       tracktree1.Fill();
   }   

   tracktree1.Write();

   pfile->Close();

   tarray.Clear();

   tf->Close();

   delete ca;
   cf->Close();

   printf("after cf close\n");


/////////////////////////////////////////////////////////////////////////
   GoodTrack gt[7000];
   Int_t ngood=0;

   cerr<<"Marking good tracks (this will take a while)...\n";
   ngood=good_tracks(gt,7000);
   ofstream out("good_tracks");
   if (out) {
         for (Int_t ngd=0; ngd<ngood; ngd++)            
	    out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
	       <<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
	       <<gt[ngd].x  <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
	       <<gt[ngd].pxg  <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
	       <<gt[ngd].ptg <<gt[ngd].flag<<endl;
   } else cerr<<"Can not open file (good_tracks) !\n";
   out.close();

   cerr<<"Number of good tracks : "<<ngood<<endl;

   printf("before return in tpctracks\n");

   return 0;
}

//---------------------------------

Int_t good_tracks(GoodTrack *gt, Int_t max) {
   Int_t nt=0;

   //TFile *file=TFile::Open("rfio:galice.root");  // via server
   TFile *file=TFile::Open("galice.root");
   if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}

   if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
     cerr<<"gAlice have not been found on galice.root !\n";
     exit(5);
   }

   gAlice->GetEvent(0);   

   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
   Int_t ver = TPC->IsVersion(); 
   cerr<<"TPC version "<<ver<<" has been found !\n";

   AliTPCParam *digp=(AliTPCParam*)file->Get("75x40_100x60");
   if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
   TPC->SetParam(digp);

   Int_t nrow_up=digp->GetNRowUp();
   Int_t nrows=digp->GetNRowLow()+nrow_up;
   Int_t zero=digp->GetZeroSup();
   Int_t gap=Int_t(0.125*nrows);
   Int_t good_number=Int_t(0.4*nrows);

   TClonesArray *particles=gAlice->Particles(); 
   Int_t np=particles->GetEntriesFast();
   Int_t *good=new Int_t[np];
   for (Int_t ii=0; ii<np; ii++) good[ii]=0;

   switch (ver) {
   case 1:
     {
      TFile *cf=TFile::Open("AliTPCclusters.root");
      if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
      AliTPCClustersArray *ca=new AliTPCClustersArray;
      ca->Setup(digp);
      ca->SetClusterType("AliTPCcluster");
      ca->ConnectTree("Segment Tree");
      Int_t nrows=Int_t(ca->GetTree()->GetEntries());
      for (Int_t n=0; n<nrows; n++) {
          AliSegmentID *s=ca->LoadEntry(n);
          Int_t sec,row;
          digp->AdjustSectorRow(s->GetID(),sec,row);
          AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
          Int_t ncl=clrow.GetArray()->GetEntriesFast();
          while (ncl--) {
              AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
              Int_t lab=c->GetLabel(0);
              if (lab<0) continue; //noise cluster
              lab=TMath::Abs(lab);
              if (sec>=digp->GetNInnerSector())
              if (row==nrow_up-1    ) good[lab]|=0x1000;
              if (sec>=digp->GetNInnerSector())
              if (row==nrow_up-1-gap) good[lab]|=0x800;
              good[lab]++;
          }
          ca->ClearRow(sec,row);
      }
      cf->Close();
     }
      break;
   case 2:
      TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
      AliSimDigits da, *digits=&da;
      TD->GetBranch("Segment")->SetAddress(&digits);
      Int_t *count = new Int_t[np];
      Int_t i;
      for (i=0; i<np; i++) count[i]=0;
      Int_t sectors_by_rows=(Int_t)TD->GetEntries();
      for (i=0; i<sectors_by_rows; i++) {
          if (!TD->GetEvent(i)) continue;
          Int_t sec,row;
          digp->AdjustSectorRow(digits->GetID(),sec,row);
          cerr<<sec<<' '<<row<<"                                     \r";
          digits->First();
          while (digits->Next()) {
              Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
              Short_t dig = digits->GetDigit(it,ip);
              Int_t idx0=digits->GetTrackID(it,ip,0); 
              Int_t idx1=digits->GetTrackID(it,ip,1);
              Int_t idx2=digits->GetTrackID(it,ip,2);
              if (idx0>=0 && dig>=zero) count[idx0]+=1;
              if (idx1>=0 && dig>=zero) count[idx1]+=1;
              if (idx2>=0 && dig>=zero) count[idx2]+=1;
          }
          for (Int_t j=0; j<np; j++) {
              if (count[j]>1) {
                 if (sec>=digp->GetNInnerSector())
                 if (row==nrow_up-1    ) good[j]|=0x1000;
                 if (sec>=digp->GetNInnerSector())
                 if (row==nrow_up-1-gap) good[j]|=0x800;
                 good[j]++;
              }
              count[j]=0;
          }
      }
      delete[] count;
      break;
   default:
      cerr<<"Invalid TPC version !\n";
      file->Close();
      exit(7);
   }

   TTree *TH=gAlice->TreeH();
   TClonesArray *hits=TPC->Hits();
   Int_t npart=TH->GetEntries();

   while (npart--) {
      TPC->ResetHits();
      TH->GetEvent(npart);
      Int_t nhits=hits->GetEntriesFast();
      if (nhits==0) continue;
      AliTPChit *hit;
      Int_t j;
      for (j=0; j<nhits-1; j++) {
         hit=(AliTPChit*)hits->UncheckedAt(j);
         if (hit->fQ==0.) break;
      }
      if (j==nhits-1) continue;
      AliTPChit *hit1=(AliTPChit*)hits->UncheckedAt(j+1);
      if (hit1->fQ != 0.) continue;
     // Int_t i=hit->fTrack;
	   Int_t i=hit->Track();  //modificata in accordo a nuovo AliTPCComparison
      TParticle *p = (TParticle*)particles->UncheckedAt(i);
      if (p->GetFirstMother()>=0) continue;  //secondary particle
      if (good[i] < 0x1000+0x800+2+good_number) continue;
      if (p->Pt()<0.100) continue;
      if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;

      gt[nt].lab=i;
      gt[nt].code=p->GetPdgCode();
//**** px py pz - in global coordinate system, x y z - in local !
     // gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ;  //modificato tenendo conto di AliTPCComparison
      gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();	  
      Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
      //gt[nt].x = hit1->fX*cs + hit1->fY*sn;
     // gt[nt].y =-hit1->fX*sn + hit1->fY*cs;  //modificato tenedo conto di AliTPCComaprison
      //gt[nt].z = hit1->fZ;
      gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
      gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
      gt[nt].z = hit1->Z();		
      gt[nt].pxg = p->Px();
      gt[nt].pyg = p->Py();
      gt[nt].pzg = p->Pz();				
      gt[nt].ptg = p->Pt();		
      gt[nt].flag = 0;		
      nt++;

      cerr<<i<<"                \r";
      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
   }
   delete[] good;

   printf("before delete gAlice\n");

   delete gAlice; gAlice=0;

   printf("after delete gAlice\n");
   file->Close();
   printf("after file close\n");
   return nt;
}

//--------------------------------------
Commit	Line	Data
d1c69928	1	struct GoodTrack {
	2	Int_t lab;
	3	Int_t code;
	4	Float_t px,py,pz;
	5	Float_t x,y,z;
	6	Float_t pxg,pyg,pzg,ptg;
	7	Bool_t flag;
	8	};
	9	Int_t good_tracks(GoodTrack *gt, Int_t max);
	10
	11
	12	Int_t TPCtracks() {
	13	cerr<<"Doing comparison...\n";
	14
	15	// Connect the Root Galice file containing Geometry, Kine and Hits
	16
	17
	18
	19	if (gClassTable->GetID("AliRun") < 0) {
	20	gROOT->LoadMacro("loadlibs.C");
	21	loadlibs();
	22	}
	23
	24
	25	TFile *cf=TFile::Open("AliTPCclusters.root");
	26	if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
	27	AliTPCParam digp= (AliTPCParam)cf->Get("75x40_100x60");
	28	if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
	29
	30	// Load clusters
	31	AliTPCClustersArray *ca=new AliTPCClustersArray;
	32	ca->Setup(digp);
	33	ca->SetClusterType("AliTPCcluster");
	34	ca->ConnectTree("Segment Tree");
	35	Int_t nentr=Int_t(ca->GetTree()->GetEntries());
	36	for (Int_t i=0; i<nentr; i++) {
	37	AliSegmentID *s=ca->LoadEntry(i);
	38	}
	39
	40	// Load tracks
	41	TFile *tf=TFile::Open("AliTPCtracks.root");
	42	if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
	43	TObjArray tarray(2000);
	44	TTree tracktree=(TTree)tf->Get("TreeT");
	45	TBranch *tbranch=tracktree->GetBranch("tracks");
	46	Int_t nentr=tracktree->GetEntries();
	47	for (i=0; i<nentr; i++) {
	48	AliTPCtrack *iotrack=new AliTPCtrack;
	49	tbranch->SetAddress(&iotrack);
	50	tracktree->GetEvent(i);
	51	iotrack->CookLabel(ca);
	52	tarray.AddLast(iotrack);
	53	}
	54
	55	cerr<<"Number of found tracks "<<nentr<<endl;
	56	tarray.Sort();
	57
	58	TFile *pfile = new TFile("tpctracks.root","RECREATE");
	59
	60	TTree tracktree1("TreeT","Tree with TPC tracks");
	61	tracktree1.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
	62
	63	for (i=0; i<nentr; i++) {
	64	iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
65	tracktree1.Fill();
66	}
67
68	tracktree1.Write();
69
70	pfile->Close();
71
72	tarray.Clear();
73
74	tf->Close();
75
76	delete ca;
77	cf->Close();
78
79	printf("after cf close\n");
80
81
82	/////////////////////////////////////////////////////////////////////////
83	GoodTrack gt[7000];
84	Int_t ngood=0;
85
86	cerr<<"Marking good tracks (this will take a while)...\n";
87	ngood=good_tracks(gt,7000);
88	ofstream out("good_tracks");
89	if (out) {
90	for (Int_t ngd=0; ngd<ngood; ngd++)
91	out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
92	<<gt[ngd].px <<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
93	<<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<' '
94	<<gt[ngd].pxg <<' '<<gt[ngd].pyg <<' '<<gt[ngd].pzg <<' '
95	<<gt[ngd].ptg <<gt[ngd].flag<<endl;
96	} else cerr<<"Can not open file (good_tracks) !\n";
97	out.close();
98
99	cerr<<"Number of good tracks : "<<ngood<<endl;
100
101	printf("before return in tpctracks\n");
102
103	return 0;
104	}
105
106	//---------------------------------
107
108	Int_t good_tracks(GoodTrack *gt, Int_t max) {
109	Int_t nt=0;
110
111	//TFile *file=TFile::Open("rfio:galice.root"); // via server
112	TFile *file=TFile::Open("galice.root");
113	if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
114
115	if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
116	cerr<<"gAlice have not been found on galice.root !\n";
117	exit(5);
118	}
119
120	gAlice->GetEvent(0);
121
122	AliTPC TPC=(AliTPC)gAlice->GetDetector("TPC");
123	Int_t ver = TPC->IsVersion();
124	cerr<<"TPC version "<<ver<<" has been found !\n";
125
126	AliTPCParam digp=(AliTPCParam)file->Get("75x40_100x60");
127	if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
128	TPC->SetParam(digp);
129
130	Int_t nrow_up=digp->GetNRowUp();
131	Int_t nrows=digp->GetNRowLow()+nrow_up;
132	Int_t zero=digp->GetZeroSup();
133	Int_t gap=Int_t(0.125*nrows);
134	Int_t good_number=Int_t(0.4*nrows);
135
136	TClonesArray *particles=gAlice->Particles();
137	Int_t np=particles->GetEntriesFast();
138	Int_t *good=new Int_t[np];
139	for (Int_t ii=0; ii<np; ii++) good[ii]=0;
140
141	switch (ver) {
142	case 1:
143	{
144	TFile *cf=TFile::Open("AliTPCclusters.root");
145	if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
146	AliTPCClustersArray *ca=new AliTPCClustersArray;
147	ca->Setup(digp);
148	ca->SetClusterType("AliTPCcluster");
149	ca->ConnectTree("Segment Tree");
150	Int_t nrows=Int_t(ca->GetTree()->GetEntries());
151	for (Int_t n=0; n<nrows; n++) {
152	AliSegmentID *s=ca->LoadEntry(n);
153	Int_t sec,row;
154	digp->AdjustSectorRow(s->GetID(),sec,row);
155	AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
156	Int_t ncl=clrow.GetArray()->GetEntriesFast();
157	while (ncl--) {
158	AliTPCcluster c=(AliTPCcluster)clrow[ncl];
159	Int_t lab=c->GetLabel(0);
160	if (lab<0) continue; //noise cluster
161	lab=TMath::Abs(lab);
162	if (sec>=digp->GetNInnerSector())
163	if (row==nrow_up-1 ) good[lab]\|=0x1000;
164	if (sec>=digp->GetNInnerSector())
165	if (row==nrow_up-1-gap) good[lab]\|=0x800;
166	good[lab]++;
167	}
168	ca->ClearRow(sec,row);
169	}
170	cf->Close();
171	}
172	break;
173	case 2:
174	TTree TD=(TTree)gDirectory->Get("TreeD_75x40_100x60");
175	AliSimDigits da, *digits=&da;
176	TD->GetBranch("Segment")->SetAddress(&digits);
177	Int_t *count = new Int_t[np];
178	Int_t i;
179	for (i=0; i<np; i++) count[i]=0;
180	Int_t sectors_by_rows=(Int_t)TD->GetEntries();
181	for (i=0; i<sectors_by_rows; i++) {
182	if (!TD->GetEvent(i)) continue;
183	Int_t sec,row;
184	digp->AdjustSectorRow(digits->GetID(),sec,row);
185	cerr<<sec<<' '<<row<<" \r";
186	digits->First();
187	while (digits->Next()) {
188	Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
189	Short_t dig = digits->GetDigit(it,ip);
190	Int_t idx0=digits->GetTrackID(it,ip,0);
191	Int_t idx1=digits->GetTrackID(it,ip,1);
192	Int_t idx2=digits->GetTrackID(it,ip,2);
193	if (idx0>=0 && dig>=zero) count[idx0]+=1;
194	if (idx1>=0 && dig>=zero) count[idx1]+=1;
195	if (idx2>=0 && dig>=zero) count[idx2]+=1;
196	}
197	for (Int_t j=0; j<np; j++) {
198	if (count[j]>1) {
199	if (sec>=digp->GetNInnerSector())
200	if (row==nrow_up-1 ) good[j]\|=0x1000;
201	if (sec>=digp->GetNInnerSector())
202	if (row==nrow_up-1-gap) good[j]\|=0x800;
203	good[j]++;
204	}
205	count[j]=0;
206	}
207	}
208	delete[] count;
209	break;
210	default:
211	cerr<<"Invalid TPC version !\n";
212	file->Close();
213	exit(7);
214	}
215
216	TTree *TH=gAlice->TreeH();
217	TClonesArray *hits=TPC->Hits();
218	Int_t npart=TH->GetEntries();
219
220	while (npart--) {
221	TPC->ResetHits();
222	TH->GetEvent(npart);
223	Int_t nhits=hits->GetEntriesFast();
224	if (nhits==0) continue;
225	AliTPChit *hit;
226	Int_t j;
227	for (j=0; j<nhits-1; j++) {
228	hit=(AliTPChit*)hits->UncheckedAt(j);
229	if (hit->fQ==0.) break;
230	}
231	if (j==nhits-1) continue;
232	AliTPChit hit1=(AliTPChit)hits->UncheckedAt(j+1);
233	if (hit1->fQ != 0.) continue;
234	// Int_t i=hit->fTrack;
235	Int_t i=hit->Track(); //modificata in accordo a nuovo AliTPCComparison
236	TParticle p = (TParticle)particles->UncheckedAt(i);
237	if (p->GetFirstMother()>=0) continue; //secondary particle
238	if (good[i] < 0x1000+0x800+2+good_number) continue;
239	if (p->Pt()<0.100) continue;
240	if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
241
242	gt[nt].lab=i;
243	gt[nt].code=p->GetPdgCode();
244	//**** px py pz - in global coordinate system, x y z - in local !
245	// gt[nt].px=hit->fX; gt[nt].py=hit->fY; gt[nt].pz=hit->fZ; //modificato tenendo conto di AliTPCComparison
246	gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
247	Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
248	//gt[nt].x = hit1->fXcs + hit1->fYsn;
249	// gt[nt].y =-hit1->fXsn + hit1->fYcs; //modificato tenedo conto di AliTPCComaprison
250	//gt[nt].z = hit1->fZ;
251	gt[nt].x = hit1->X()cs + hit1->Y()sn;
252	gt[nt].y =-hit1->X()sn + hit1->Y()cs;
253	gt[nt].z = hit1->Z();
254	gt[nt].pxg = p->Px();
255	gt[nt].pyg = p->Py();
256	gt[nt].pzg = p->Pz();
257	gt[nt].ptg = p->Pt();
258	gt[nt].flag = 0;
259	nt++;
260
261	cerr<<i<<" \r";
262	if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
263	}
264	delete[] good;
265
266	printf("before delete gAlice\n");
267
268	delete gAlice; gAlice=0;
269
270	printf("after delete gAlice\n");
271	file->Close();
272	printf("after file close\n");
273	return nt;
274	}
275
276	//--------------------------------------
277