nvartrk=7 in ESE task
[u/mrichter/AliRoot.git] / TPC / AliTPCComparison2.C
CommitLineData
6e9ab565 1/****************************************************************************
2 * Very important, delicate and rather obscure macro. *
3 * *
4 * Creates list of "trackable" tracks, *
5 * sorts tracks for matching with the ITS, *
6 * calculates efficiency, resolutions etc. *
7 * *
8 * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
9 * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de *
10 ****************************************************************************/
11
12#if !defined(__CINT__) || defined(__MAKECINT__)
13#include<fstream.h>
14#include<TROOT.h>
15#include<TCanvas.h>
16#include<TClassTable.h>
17#include<TClonesArray.h>
18#include<TFile.h>
19#include<TFolder.h>
20#include<TH1.h>
21#include<TObject.h>
22#include<TObjArray.h>
23#include<TParticle.h>
24#include<TTree.h>
25#include <AliMagF.h>
26#include <AliRun.h>
6e9ab565 27#include <alles.h>
28#endif
29
30/*
31 multievent comparison of reconstructed tracks with "exact tracks"
32 from silulation
33
34*/
35
36
37struct GoodTrackTPC {
38 Int_t fEventN; //event number
39 Int_t lab;
40 Int_t code;
41 Float_t px,py,pz;
42 Float_t x,y,z;
43};
44
45enum tagprimary {kPrimaryCharged = 0x4000};
46
47Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn);
48
49Int_t FindPrimaries(Int_t nparticles);
50
51Int_t AliTPCComparison2(Int_t firstev=0, Int_t eventn=1) {
52
53 const Int_t MAX = 20000;
54 const Bool_t kOLD = kFALSE; // if tracking was done with a previous version
55 /***********************************************************************/
56
57 TFile *inkin = TFile::Open("rfio:galice.root");
58 // if(gAlice)delete gAlice; COMMENTED BECAUSE OF A BUG (IN COMPILED MODE)
59 gAlice = (AliRun*)inkin->Get("gAlice");
60 cout<<"AliRun object found on file "<<gAlice<<endl;
c84a5e9e 61 AliKalmanTrack::SetFieldMap(gAlice->Field());
6e9ab565 62 inkin->Close();
63 /*
64 delete gAlice; COMMENTED BECAUSE OF A BUG IN COMPILED MODE
65 gAlice = 0;
66 */
67 /***********************************************************************/
68 cerr<<"Doing comparison...\n";
69 Int_t i;
70 gBenchmark->Start("AliTPCComparison2");
71
72 TFile *cf=0;
73 AliTPCParam *digp=0;
74 if(kOLD){
75 cf=TFile::Open("AliTPCclusters.root");
76 if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
7a09f434 77 digp= (AliTPCParam*)cf->Get("75x40_100x60_150x60");
6e9ab565 78 if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
79 }
80 ///////////
6e9ab565 81 TObjArray tarray(MAX);
82 AliTPCtrack *iotrack=0;
83 Int_t nentr= 0;
84 Int_t eventptr[1000];
85 TFile *tf=TFile::Open("AliTPCtracks.root");
86 TTree *tracktree=0;
87
88 eventptr[firstev]=0;
89 eventptr[firstev+1]=0;
90 for (Int_t event=firstev;event<eventn; event++){
91 cout<<"================================================"<<endl;
92 cout<<"Processing event "<<event<<endl;
93 cout<<"================================================"<<endl;
6e9ab565 94
95 char tname[100];
96 if (eventn==-1) {
97 sprintf(tname,"TreeT_TPC");
98 }
99 else {
100 sprintf(tname,"TreeT_TPC_%d",event);
101 }
102
103 // Load tracks
104 if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
105
106 tracktree=(TTree*)tf->Get(tname);
107 if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
108 TBranch *tbranch=tracktree->GetBranch("tracks");
109 Int_t nentr0=(Int_t)tracktree->GetEntries();
110 nentr+=nentr0;
111 for (i=0; i<nentr0; i++) {
112 iotrack=new AliTPCtrack;
113 tbranch->SetAddress(&iotrack);
114 tracktree->GetEvent(i);
6e9ab565 115 tarray.AddLast(iotrack);
116 }
117 eventptr[event+1] = nentr; //store end of the event
6e9ab565 118 delete tracktree;
119 }
120 tf->Close();
121 if(kOLD)cf->Close();
122
123 GoodTrackTPC gt[MAX];
124
125 Int_t ngood=0;
126 ifstream in("good_tracks_tpc");
127 if (in) {
128 cerr<<"Reading good tracks...\n";
129 while (in>>gt[ngood].fEventN>>gt[ngood].lab>>gt[ngood].code>>
130 gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
131 gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
132 ngood++;
133 cerr<<ngood<<"\r";
7a09f434 134 //cout<<ngood<<"\r";
6e9ab565 135 if (ngood==MAX) {
136 cerr<<"Too many good tracks !\n";
137 break;
138 }
139 }
140 cout<<endl;
141 if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
142 } else {
143 cerr<<"Marking good tracks (this will take a while)...\n";
144 ngood=good_tracks(gt,45000,firstev,eventn);
7a09f434 145 printf("Goood %d\n", ngood);
6e9ab565 146 ofstream out("good_tracks_tpc");
6e9ab565 147 if (out) {
148 cout<<"File good_tracks_tpc opened\n";
149 for (Int_t ngd=0; ngd<ngood; ngd++) {
6e9ab565 150 out<<gt[ngd].fEventN<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
151 gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '<<
152 gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
6e9ab565 153 }
154 } else cerr<<"Can not open file (good_tracks_tpc) !\n";
7a09f434 155 out<<flush;
6e9ab565 156 out.close();
7a09f434 157
158 ofstream out2("good_tracks_tpc_par");
159
160 if (out2) {
161 //cout<<"File good_tracks_tpc opened\n";
162 for (Int_t ngd=0; ngd<ngood; ngd++) {
163 Float_t pt = TMath::Sqrt(gt[ngd].px*gt[ngd].px+gt[ngd].py*gt[ngd].py);
164 Float_t angle = 0;
165 if (TMath::Abs(pt)>0.01) angle = gt[ngd].pz/pt;
166 out2<<gt[ngd].fEventN<<"\t"<<gt[ngd].lab<<"\t"<<gt[ngd].code<<"\t"<<
167 pt<<"\t"<<angle<<"\t"<<endl;
168 }
169 } else cerr<<"Can not open file (good_tracks_tpc) !\n";
170 out2<<flush;
6e9ab565 171 out2.close();
7a09f434 172
6e9ab565 173 }
174 cerr<<"Number of good tracks : "<<ngood<<endl;
7a09f434 175 cout<<"Number of good tracks : "<<ngood<<endl;
6e9ab565 176 if(ngood==0)return 5;
7a09f434 177 TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
178 TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
6e9ab565 179 TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
180 hpt->SetFillColor(2);
181 TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60);
182 hmpt->SetFillColor(6);
183
184 AliTPCtrack *trk=(AliTPCtrack*)tarray.UncheckedAt(0);
185
186 Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
187 Double_t pmax=6.0+pmin;
188
189 TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
190 TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
191 TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
192 TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
193 hg->SetLineColor(4); hg->SetLineWidth(2);
194 TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax);
195 hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
196
197 TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
198 TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,250.);
199 hep->SetMarkerStyle(8);
200 hep->SetMarkerSize(0.4);
201
202 Int_t mingood=ngood;
203 Int_t track_notfound[MAX], itrack_notfound=0;
204 Int_t track_fake[MAX], itrack_fake=0;
205 Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound=0;
206 while (ngood--) {
207 Int_t lab=gt[ngood].lab,tlab=-1;
208 Float_t ptg=
209 TMath::Sqrt(gt[ngood].px*gt[ngood].px + gt[ngood].py*gt[ngood].py);
210
211 if (ptg<pmin) continue; // Should be 1e-33 instead pmin?
212
213 hgood->Fill(ptg);
214 Int_t ievent = gt[ngood].fEventN;
215
216 AliTPCtrack *track=0;
217 for (i=eventptr[ievent]; i<eventptr[ievent+1]; i++) {
218
219 track=(AliTPCtrack*)tarray.UncheckedAt(i);
220 tlab=track->GetLabel();
221 if (lab==TMath::Abs(tlab)) break;
222 }
223 if (i==eventptr[ievent+1]) {
224 track_notfound[itrack_notfound++]=lab;
225 continue;
226 }
227
228 Int_t micount=0;
229 Int_t mi;
230 AliTPCtrack * mitrack;
231 for (mi=eventptr[ievent]; mi<eventptr[ievent+1]; mi++) {
232
233 mitrack=(AliTPCtrack*)tarray.UncheckedAt(mi);
234 if (lab==TMath::Abs(mitrack->GetLabel())) micount++;
235 }
236 if (micount>1) {
237 track_multifound[itrack_multifound]=lab;
238 track_multifound_n[itrack_multifound]=micount;
239 itrack_multifound++;
240 }
241
242 //
243 Double_t xk=gt[ngood].x;
7a09f434 244 if (!track) continue;
245 // printf("Track =%p\n",track);
6e9ab565 246 track->PropagateTo(xk);
247
248 if (lab==tlab) hfound->Fill(ptg);
249 else {
250 track_fake[itrack_fake++]=lab;
251 hfake->Fill(ptg);
252 }
253 Double_t par[5]; track->GetExternalParameters(xk,par);
254 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
255 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
256 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
257 Float_t lam=TMath::ATan(par[3]);
258 Float_t pt_1=TMath::Abs(par[4]);
259
260 if (TMath::Abs(gt[ngood].code)==11 && ptg>4.) {
261 hmpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
7a09f434 262 }
263 //else
264 {
6e9ab565 265 Float_t phig=TMath::ATan2(gt[ngood].py,gt[ngood].px);
266 hp->Fill((phi - phig)*1000.);
267
268 Float_t lamg=TMath::ATan2(gt[ngood].pz,ptg);
269 hl->Fill((lam - lamg)*1000.);
270
271 hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
272 }
273
274 Float_t mom=1./(pt_1*TMath::Cos(lam));
275 Float_t dedx=track->GetdEdx();
276 hep->Fill(mom,dedx,1.);
277 if (TMath::Abs(gt[ngood].code)==211)
278 if (mom>0.4 && mom<0.5) {
279 he->Fill(dedx,1.);
280 }
281
282 }
283
284 cout<<"\nList of Not found tracks :\n";
285 for ( i = 0; i< itrack_notfound; i++){
286 cout<<track_notfound[i]<<"\t";
287 if ((i%5)==4) cout<<"\n";
288 }
289 cout<<"\nList of fake tracks :\n";
290 for ( i = 0; i< itrack_fake; i++){
291 cout<<track_fake[i]<<"\t";
292 if ((i%5)==4) cout<<"\n";
293 }
294 cout<<"\nList of multiple found tracks :\n";
295 for ( i=0; i<itrack_multifound; i++) {
296 cout<<"id. "<<track_multifound[i]
297 <<" found - "<<track_multifound_n[i]<<"times\n";
298 }
299
300 Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
301 if (ng!=0) cerr<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
7a09f434 302 if (ng!=0) cout<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
6e9ab565 303 cout<<"Total number of found tracks ="<<nentr<<endl;
304 cout<<"Total number of \"good\" tracks ="
305 <<mingood<<" (selected for comparison: "<<ng<<')'<<endl<<endl;
306
307
308 gStyle->SetOptStat(111110);
309 gStyle->SetOptFit(1);
310
311 TCanvas *c1=new TCanvas("c1","",0,0,700,850);
312
313 TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
314 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
315 hp->SetFillColor(4); hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd();
316
317 TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
318 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
319 hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd();
320
321 TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
322 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
323 hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd();
324
325 TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
326 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
327 hmpt->SetXTitle("(%)"); hmpt->Fit("gaus"); c1->cd();
328
329 TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
330 p5->SetFillColor(41); p5->SetFrameFillColor(10);
331 hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
332 hg->Divide(hfound,hgood,1,1.,"b");
333 hf->Divide(hfake,hgood,1,1.,"b");
334 hg->SetMaximum(1.4);
335 hg->SetYTitle("Tracking efficiency");
336 hg->SetXTitle("Pt (GeV/c)");
337 hg->Draw();
338
339 TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
340 line1->Draw("same");
341 TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
342 line2->Draw("same");
343
344 hf->SetFillColor(1);
345 hf->SetFillStyle(3013);
346 hf->SetLineColor(2);
347 hf->SetLineWidth(2);
348 hf->Draw("histsame");
349 TText *text = new TText(0.461176,0.248448,"Fake tracks");
350 text->SetTextSize(0.05);
351 text->Draw();
352 text = new TText(0.453919,1.11408,"Good tracks");
353 text->SetTextSize(0.05);
354 text->Draw();
355
356
357
358 TCanvas *c2=new TCanvas("c2","",320,32,530,590);
359
360 TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw();
361 p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
362 he->SetFillColor(2); he->SetFillStyle(3005);
363 he->SetXTitle("Arbitrary Units");
364 he->Fit("gaus"); c2->cd();
365
366 TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw();
367 p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
368 hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
369 hep->Draw(); c1->cd();
370
371 gBenchmark->Stop("AliTPCComparison2");
372 gBenchmark->Show("AliTPCComparison2");
373
374
375 return 0;
376}
377
378
379Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn) {
380 //eventn - number of events in file
381
7a09f434 382 TFile *file=TFile::Open("galice.root");
6e9ab565 383 if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
384 // delete gAlice; gAlice = 0;
385 if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
386 cerr<<"gAlice have not been found on galice.root !\n";
387 exit(5);
388 }
1c53abe2 389
390 TFile *fdigit = TFile::Open("digits.root");
391 file->cd();
6e9ab565 392
393 AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
394 Int_t ver = TPC->IsVersion();
395 cerr<<"TPC version "<<ver<<" has been found !\n";
396
397 AliTPCParam *digp=(AliTPCParam*)file->Get("75x40_100x60");
7a09f434 398 if(digp){
399 cerr<<"2 pad-lenght geom hits with 3 pad-length geom digits...\n";
400 delete digp;
401 digp = new AliTPCParamSR();
402 }
403 else
404 {
405 digp =(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60");
406 }
407
6e9ab565 408 if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
409 TPC->SetParam(digp);
410
411 Int_t nrow_up=digp->GetNRowUp();
412 Int_t nrows=digp->GetNRowLow()+nrow_up;
413 Int_t zero=digp->GetZeroSup();
414 Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
415 Int_t good_number=Int_t(0.4*nrows);
416
417 Int_t nt=0; //reset counter
418 char treeName[100]; //declare event identifier
419
420 for (Int_t event=firstev;event<eventn;event++){
421 Int_t np=gAlice->GetEvent(event);
422 Int_t nopr = FindPrimaries(np);
423 cout<<"function good_tracks -- event "<<event<<", Charged primaries:"<<nopr<<"\n";
424
425 Int_t *good=new Int_t[np];
426 for (Int_t ii=0; ii<np; ii++) good[ii]=0;
427
428
429 Int_t sectors_by_rows=0;
430 TTree *TD=0;
431 AliSimDigits da, *digits=&da;
432 Int_t *count=0;
433 switch (ver) {
434 case 1:
435 {
436 TFile *cf=TFile::Open("AliTPCclusters.root");
437 if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
438 AliTPCClustersArray *ca=new AliTPCClustersArray;
439 ca->Setup(digp);
440 ca->SetClusterType("AliTPCcluster");
441 ca->ConnectTree("TreeC_TPC_0");
442 Int_t nrows=Int_t(ca->GetTree()->GetEntries());
443 for (Int_t n=0; n<nrows; n++) {
444 AliSegmentID *s=ca->LoadEntry(n);
445 Int_t sec,row;
446 digp->AdjustSectorRow(s->GetID(),sec,row);
447 AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
448 Int_t ncl=clrow.GetArray()->GetEntriesFast();
449 while (ncl--) {
450 AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
451 Int_t lab=c->GetLabel(0);
452 if (lab<0) continue; //noise cluster
453 lab=TMath::Abs(lab);
454 if (sec>=digp->GetNInnerSector())
455 if (row==nrow_up-1 ) good[lab]|=0x1000;
456 if (sec>=digp->GetNInnerSector())
457 if (row==nrow_up-1-gap) good[lab]|=0x800;
458 good[lab]++;
459 }
460 ca->ClearRow(sec,row);
461 }
462 delete ca;
463 cf->Close();
464 }
465 break;
466 case 2:
467 {
7a09f434 468 sprintf(treeName,"TreeD_75x40_100x60_150x60_%d",event);
88cb7938 469 TD=(TTree*)fdigit->Get(treeName); // To be revised according to
470 // NewIO schema M.Kowalski
6e9ab565 471 TD->GetBranch("Segment")->SetAddress(&digits);
472 count = new Int_t[np];
473 Int_t i;
474 for (i=0; i<np; i++) count[i]=0;
475 sectors_by_rows=(Int_t)TD->GetEntries();
476 for (i=0; i<sectors_by_rows; i++) {
477 if (!TD->GetEvent(i)) continue;
478 Int_t sec,row;
479 digp->AdjustSectorRow(digits->GetID(),sec,row);
480 digits->First();
1c53abe2 481 digits->ExpandTrackBuffer();
6e9ab565 482 do { //Many thanks to J.Chudoba who noticed this
483 Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
1c53abe2 484 // Short_t dig = digits->GetDigit(it,ip);
485 Short_t dig = digits->CurrentDigit();
486
487 Int_t idx0=digits->GetTrackIDFast(it,ip,0)-2;
488 Int_t idx1=digits->GetTrackIDFast(it,ip,1)-2;
489 Int_t idx2=digits->GetTrackIDFast(it,ip,2)-2;
6e9ab565 490 if (idx0>=0 && dig>=zero && idx0<np ) count[idx0]+=1; //background events - higher track ID's
491 if (idx1>=0 && dig>=zero && idx1<np ) count[idx1]+=1;
492 if (idx2>=0 && dig>=zero && idx2<np ) count[idx2]+=1;
493 } while (digits->Next());
494 for (Int_t j=0; j<np; j++) {
495 if (count[j]>1) {
496 if (sec>=digp->GetNInnerSector())
497 if (row==nrow_up-1 ) good[j]|=0x4000;
498 if (sec>=digp->GetNInnerSector())
499 if (row==nrow_up-1-gap) good[j]|=0x1000;
500
501 if (sec>=digp->GetNInnerSector())
502 if (row==nrow_up-1-shift) good[j]|=0x2000;
503 if (sec>=digp->GetNInnerSector())
504 if (row==nrow_up-1-gap-shift) good[j]|=0x800;
505 good[j]++;
506 }
507 count[j]=0;
508 }
509
510 }
511 delete[] count;
512 delete TD; //Thanks to Mariana Bondila
513 }
514 break;
515 default:
516 cerr<<"Invalid TPC version !\n";
517 file->Close();
518 exit(7);
519 }
520
521 /** select tracks which are "good" enough **/
7a09f434 522 //printf("\t %d \n",np);
6e9ab565 523 for (Int_t i=0; i<np; i++) {
524 if ((good[i]&0x5000) != 0x5000)
525 if ((good[i]&0x2800) != 0x2800) continue;
526 if ((good[i]&0x7FF ) < good_number) continue;
527
528 TParticle *p = (TParticle*)gAlice->Particle(i);
529
530 if (p->Pt()<0.100) continue;
531 if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
532 // THIS IS GOOD IF YOU WANT SECONDARIES
533 Int_t j=p->GetFirstMother();
534 if (j>=0) {
535 TParticle *pp = (TParticle*)gAlice->Particle(j);
536 //if (!(pp->TestBit(kPrimaryCharged))) continue; //only one decay is allowed
537 }
538
7a09f434 539 if(!(p->TestBit(kPrimaryCharged)))continue; // only primaries
540 // printf("1");
6e9ab565 541
542 gt[nt].fEventN=event;
543 gt[nt].lab=i;
544 gt[nt].code=p->GetPdgCode();
545 gt[nt].px=0.; gt[nt].py=0.; gt[nt].pz=0.;
527abbe9 546 gt[nt].x=0.; gt[nt].y=0.; gt[nt].z=0.;
6e9ab565 547 nt++;
548 if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
549 cerr<<np-i<<" \r";
550 }
551
552 /** check if there is also information at the entrance of the TPC **/
553 TTree *TH=gAlice->TreeH(); np=(Int_t)TH->GetEntries();
554 cout<<endl;
555 for (Int_t i=0; i<np; i++) {
556 TPC->ResetHits();
557 TH->GetEvent(i);
558 AliTPChit *phit = (AliTPChit*)TPC->FirstHit(-1);
559 for ( ; phit; phit=(AliTPChit*)TPC->NextHit() ) {
560 if (phit->fQ !=0. ) continue;
561
562 Double_t px=phit->X(), py=phit->Y(), pz=phit->Z();
563
564 if ((phit=(AliTPChit*)TPC->NextHit())==0) break;
565 if (phit->fQ != 0.) continue;
566
567 Double_t x=phit->X(), y=phit->Y(), z=phit->Z();
568 if (TMath::Sqrt(x*x+y*y)>90.) continue;
569
570 Int_t j, lab=phit->Track();
571 for (j=0; j<nt; j++) {if (gt[j].fEventN==event && gt[j].lab==lab) break;}
572 if (j==nt) continue;
7a09f434 573 //printf("1-");
6e9ab565 574 // (px,py,pz) - in global coordinate system, (x,y,z) - in local !
575 gt[j].px=px; gt[j].py=py; gt[j].pz=pz;
576 Float_t cs,sn; digp->AdjustCosSin(phit->fSector,cs,sn);
577 gt[j].x = x*cs + y*sn;
578 gt[j].y =-x*sn + y*cs;
579 gt[j].z = z;
580 }
581 cerr<<np-i<<" \r";
582 }
7a09f434 583 //printf("\n%d\n",nt);
6e9ab565 584 cout<<endl;
585 delete[] good;
586 } // /// loop on events
587 delete gAlice; gAlice=0;
588
589 file->Close();
590 gBenchmark->Stop("AliTPCComparison2");
591 gBenchmark->Show("AliTPCComparison2");
592 return nt;
593
594}
595
596Int_t FindPrimaries(Int_t nparticles){
597
598 // cuts:
599 Double_t vertcut = 0.001;
600 Double_t decacut = 3.;
601 Double_t timecut = 0.;
602 Int_t nprch1=0;
603 TParticle * part = gAlice->Particle(0);
604 Double_t xori = part->Vx();
605 Double_t yori = part->Vy();
606 Double_t zori = part->Vz();
607 for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
608
609 part = gAlice->Particle(iprim);
1c53abe2 610 Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
611 if (ptot<0.01) continue;
6e9ab565 612 char *xxx=strstr(part->GetName(),"XXX");
613 if(xxx)continue;
614
615 TParticlePDG *ppdg = part->GetPDG();
616 if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
617
618 Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
619 if(dist>vertcut)continue; // cut on the vertex
620
621 if(part->T()>timecut)continue;
622
1c53abe2 623 //Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
6e9ab565 624 if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
625
626 Bool_t prmch = kTRUE; // candidate primary track
627 Int_t fidau=part->GetFirstDaughter(); // cut on daughters
628 Int_t lasdau=0;
629 Int_t ndau=0;
630 if(fidau>=0){
631 lasdau=part->GetLastDaughter();
632 ndau=lasdau-fidau+1;
633 }
634 if(ndau>0){
635 for(Int_t j=fidau;j<=lasdau;j++){
636 TParticle *dau=gAlice->Particle(j);
637 Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
638 if(distd<decacut)prmch=kFALSE; // eliminate if the decay is near the vertex
639 }
640 }
641
642 if(prmch){
643 nprch1++;
644 part->SetBit(kPrimaryCharged);
645 }
646 }
647
648 return nprch1;
649}
650