Fix for memory leak in fClusters
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7 #include <TFile.h>
8 #include <TH1.h>
9 #include <TParticle.h>
10 #include <TTree.h>
11 #include <TClonesArray.h>
12
13 #include <AliRun.h>
14 #include <AliSimDigits.h>
15 #include <AliTPC.h>
16 #include <AliTPCcluster.h>
17 #include <AliTPCClustersArray.h>
18 #include <AliTPCClustersRow.h>
19 #include <AliTPCParam.h>
20 #include <AliComplexCluster.h>
21 #include <AliStack.h>
22
23 #if __GNUC__ == 3
24 #include <fstream>
25 #include <iosfwd>
26 #else
27 #include <fstream.h>
28 #endif
29
30 #include "AliL3Logging.h"
31 #include "AliL3Transform.h"
32 #include "AliL3SpacePointData.h"
33 #include "AliL3Track.h"
34 #include "AliL3FileHandler.h"
35 #include "AliL3TrackArray.h"
36 #include "AliL3Evaluate.h"
37
38 #if __GNUC__ == 3
39 using namespace std;
40 #endif
41
42 /** \class AliL3Evaluate
43 <pre>
44 //_____________________________________________________________
45 // AliL3Evaluate
46 //
47 // Evaluation class for tracking. Plots efficiencies etc..
48 //
49 </pre>
50 */
51
52 ClassImp(AliL3Evaluate)
53
54 AliL3Evaluate::AliL3Evaluate()
55 {
56   fTracks = 0;
57   fNFastPoints = 0;
58   fMcIndex = 0;
59   fMcId = 0;
60   fMinSlice=0;
61   fMaxSlice=0;
62   fGoodTracks=0;
63   fNtupleRes=0;
64   fDigitsTree=0;
65   fPtRes=0;
66   fNGoodTracksPt=0;
67   fNFoundTracksPt=0;
68   fNFakeTracksPt=0;
69   fTrackEffPt=0;
70   fFakeTrackEffPt=0;
71   fNGoodTracksEta=0;
72   fNFoundTracksEta=0;
73   fNFakeTracksEta=0;
74   fTrackEffEta=0;
75   fFakeTrackEffEta=0;
76   fMcIndex=0;
77   fMcId=0;
78   fNtuppel=0;
79   fStandardComparison=kTRUE;
80 }
81
82 AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice)
83 {
84
85   if(slice)
86     {
87       fMinSlice=slice[0];
88       fMaxSlice=slice[1];
89     }
90   else
91     {
92       fMinSlice=0;
93       fMaxSlice=35;
94     }
95   sprintf(fPath,"%s",datapath);
96   fMaxFalseClusters = 0.1;
97   fGoodFound = 0;
98   fGoodGen = 0;
99   fMinPointsOnTrack = min_clusters;
100   fMinHitsFromParticle = minhits;
101   fMinGoodPt = minpt;
102   fMaxGoodPt = maxpt;
103   memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
104   fTracks=0;
105   fGoodTracks=0;
106   fNtupleRes=0;
107   fDigitsTree=0;
108   fPtRes=0;
109   fNGoodTracksPt=0;
110   fNFoundTracksPt=0;
111   fNFakeTracksPt=0;
112   fTrackEffPt=0;
113   fFakeTrackEffPt=0;
114   fNGoodTracksEta=0;
115   fNFoundTracksEta=0;
116   fNFakeTracksEta=0;
117   fTrackEffEta=0;
118   fFakeTrackEffEta=0;
119   fMcIndex=0;
120   fMcId=0;
121   fNtuppel=0;
122   fStandardComparison=kTRUE;
123 }
124
125 void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
126 {
127   Char_t fname[1024];
128   AliL3FileHandler *clusterfile[36][6];
129   for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
130     {
131       for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
132         {
133           Int_t patch;
134           if(sp==kTRUE)
135             patch=-1;
136           else
137             patch=p;
138           
139           if(fClusters[s][p])
140             delete fClusters[s][p];
141           fClusters[s][p] = 0;
142           clusterfile[s][p] = new AliL3FileHandler();
143           if(event<0)
144             sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
145           else
146             sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
147           if(!clusterfile[s][p]->SetBinaryInput(fname))
148             {
149               LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
150                 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; 
151               delete clusterfile[s][p];
152               clusterfile[s][p] = 0; 
153               continue;
154             }
155           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
156           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
157           clusterfile[s][p]->CloseBinaryInput();
158           if(sp==kTRUE)
159             break;
160         }
161     }
162    
163   sprintf(fname,"%s/tracks_%d.raw",fPath,event);
164   AliL3FileHandler *tfile = new AliL3FileHandler();
165   if(!tfile->SetBinaryInput(fname)){
166     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
167       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG; 
168     return;
169   }
170   if(fTracks)
171     delete fTracks;
172   fTracks = new AliL3TrackArray();
173   tfile->Binary2TrackArray(fTracks);
174   tfile->CloseBinaryInput();
175   delete tfile;
176   fTracks->QSort();
177 }
178
179
180 AliL3Evaluate::~AliL3Evaluate()
181 {
182   if(fGoodTracks) delete fGoodTracks;
183   if(fDigitsTree) fDigitsTree->Delete();
184   if(fTracks) delete fTracks;
185   if(fPtRes) delete fPtRes;
186   if(fNGoodTracksPt) delete fNGoodTracksPt;
187   if(fNFoundTracksPt) delete fNFoundTracksPt;
188   if(fNFakeTracksPt) delete fNFakeTracksPt;
189   if(fTrackEffPt) delete fTrackEffPt;
190   if(fFakeTrackEffPt) delete fFakeTrackEffPt;
191   if(fNGoodTracksEta) delete fNGoodTracksEta;
192   if(fNFoundTracksEta) delete fNFoundTracksEta;
193   if(fNFakeTracksEta) delete fNFakeTracksEta;
194   if(fTrackEffEta) delete fTrackEffEta;
195   if(fFakeTrackEffEta) delete fFakeTrackEffEta;
196   if(fMcIndex) delete [] fMcIndex;
197   if(fMcId)    delete [] fMcId;
198   if(fNtuppel) delete fNtuppel;
199   if(fNtupleRes) delete fNtupleRes;
200   for(Int_t s=0; s<=35; s++)
201     for(Int_t p=0; p<6; p++)
202       if(fClusters[s][p]) delete fClusters[s][p];
203 }
204
205 void AliL3Evaluate::AssignPIDs()
206 {
207   fTracks->QSort();
208   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignPIDs","Track Loop")
209     <<"Assigning pid to the found tracks...."<<ENDLOG;
210   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
211     {
212       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
213       if(!track) continue; 
214       if(track->GetNumberOfPoints() < fMinPointsOnTrack)
215         track->SetPID(0);
216       else {
217         Float_t pid = GetTrackPID(track);
218         track->SetPID(pid);
219       }
220     }
221 }
222   
223 void AliL3Evaluate::AssignIDs()
224 {
225   //Assign MC id to the tracks.
226 #ifndef do_mc
227   cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
228   return;
229 #endif
230   fGoodFound=0;
231   fTracks->QSort();
232   LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
233     <<"Assigning MC id to the found tracks...."<<ENDLOG;
234   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
235     {
236       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
237       if(!track) continue; 
238       if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
239       
240       fGoodFound++;
241       Int_t tID = GetMCTrackLabel(track);
242       track->SetMCid(tID);
243     }
244   //cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
245 }
246
247 Float_t AliL3Evaluate::GetTrackPID(AliL3Track *track)
248 {
249   track->CalculateHelix();
250   // Track dEdx
251   Int_t nc=track->GetNHits();
252   UInt_t *hits = track->GetHitNumbers();
253   Float_t sampleDEdx[159];
254   for (Int_t iHit = 0; iHit < nc; iHit++) {
255     UInt_t hitID = hits[iHit];
256     Int_t iSector = (hitID>>25) & 0x7f;
257     Int_t patch = (hitID>>22) & 0x7;
258     UInt_t position = hitID&0x3fffff;
259     AliL3SpacePointData *points = fClusters[iSector][patch];
260     if(!points) continue; 
261     if(position>=fNcl[iSector][patch]) 
262       {
263         LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
264           <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
265         continue;
266       }
267     UChar_t padrow = points[position].fPadRow;
268     Float_t pWidth = AliL3Transform::GetPadPitchWidthLow();
269     if (padrow>63)
270       pWidth = AliL3Transform::GetPadPitchWidthUp(); 
271     Float_t corr=1.; if (padrow>63) corr=0.67;
272     sampleDEdx[iHit] = points[position].fCharge/pWidth*corr;
273     Double_t crossingangle = track->GetCrossingAngle(padrow,iSector);
274     Double_t s = sin(crossingangle);
275     Double_t t = track->GetTgl();
276     sampleDEdx[iHit] *= sqrt((1-s*s)/(1+t*t));
277   }
278
279   /* Cook dEdx */
280   Int_t i;
281   Int_t swap;//stupid sorting
282   do {
283     swap=0;
284     for (i=0; i<nc-1; i++) {
285       if (sampleDEdx[i]<=sampleDEdx[i+1]) continue;
286       Float_t tmp=sampleDEdx[i];
287       sampleDEdx[i]=sampleDEdx[i+1]; sampleDEdx[i+1]=tmp;
288       swap++;
289     }
290   } while (swap);
291
292   Double_t low=0.05; Double_t up=0.7;
293   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
294   Float_t trackDEdx=0;
295   for (i=nl; i<=nu; i++) trackDEdx += sampleDEdx[i];
296   trackDEdx /= (nu-nl+1);
297
298   //  cout<<" PID: "<<nc<<" "<<nl<<" "<<nu<<" "<<trackDEdx<<" "<<track->GetPt()<<endl;
299   return trackDEdx;
300 }
301
302 struct S {Int_t lab; Int_t max;};
303 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){ 
304   //Returns the MCtrackID of the belonging clusters.
305   //If MCLabel < 0, means that track is fake.
306   //Definitions are identical to offline.
307   //Fake track means:
308   // - more than 10 percent of clusters are assigned incorrectly
309   // - more than half of the innermost 10% of clusters were assigned incorrectly.
310   
311   
312 #ifdef do_mc
313   Int_t num_of_clusters = track->GetNumberOfPoints();
314   S *s=new S[num_of_clusters];
315   Int_t i;
316   for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
317   UInt_t *hitnum = track->GetHitNumbers();  
318   UInt_t id;
319     
320   Int_t lab=123456789;
321   for (i=0; i<num_of_clusters; i++) 
322     {
323       //Tricks to get the clusters belonging to this track:
324       id = hitnum[i];
325       Int_t slice = (id>>25) & 0x7f;
326       Int_t patch = (id>>22) & 0x7;
327       UInt_t pos = id&0x3fffff;       
328       
329       AliL3SpacePointData *points = fClusters[slice][patch];
330       if(!points) continue; 
331       if(pos>=fNcl[slice][patch]) 
332         {
333           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
334             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
335           continue;
336         }
337       
338       //Get the label of the cluster:
339       lab=points[pos].fTrackID[0];
340
341       Int_t j;
342       for (j=0; j<num_of_clusters; j++)
343         if (s[j].lab==lab || s[j].max==0) break;
344       s[j].lab=lab;
345       s[j].max++;
346     }
347   
348   Int_t max=0;
349   for (i=0; i<num_of_clusters; i++) 
350     if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
351   
352   if(lab == -1)
353     return -1; //If most clusters is -1, this is a noise track.
354   if(lab < 0)
355     cerr<<"AliL3Evaluate::GetMCTrackLabel : Track label negative :"<<lab<<endl;
356   
357   delete[] s;
358   
359   for (i=0; i<num_of_clusters; i++) 
360     {
361       id = hitnum[i];
362       Int_t slice = (id>>25) & 0x7f;
363       Int_t patch = (id>>22) & 0x7;
364       UInt_t pos = id&0x3fffff;       
365       
366       AliL3SpacePointData *points = fClusters[slice][patch];
367       if(!points) continue; 
368       if(pos>=fNcl[slice][patch]) 
369         {
370           LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
371             <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
372           continue;
373         }
374       
375       if (abs(points[pos].fTrackID[1]) == lab || 
376           abs(points[pos].fTrackID[2]) == lab ) max++;
377     }
378   
379   
380   //Check if more than 10% of the clusters were assigned incorrectly:
381   if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters) 
382     {
383       return -lab;
384     }
385   else //Check if at least half of the 10% innermost clusters are assigned correctly.
386     {
387       Int_t tail=Int_t(0.10*num_of_clusters);
388       max=0;
389       for (i=1; i<=tail; i++) 
390         {
391           id = hitnum[num_of_clusters - i];
392           Int_t slice = (id>>25) & 0x7f;
393           Int_t patch = (id>>22) & 0x7;
394           UInt_t pos = id&0x3fffff;           
395           
396           AliL3SpacePointData *points = fClusters[slice][patch];
397           if(lab == abs(points[pos].fTrackID[0]) ||
398              lab == abs(points[pos].fTrackID[1]) ||
399              lab == abs(points[pos].fTrackID[2])) max++;
400         }
401       if (max < Int_t(0.5*tail)) return -lab;
402     }
403
404   return lab;
405 #else //If we are running with mc_ids or not
406   return 0;
407 #endif
408
409 }
410
411 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
412 {
413   //Get the MC id of space points in case of using the fast simulator. 
414   char fname[256];
415   sprintf(fname,"%s/point_mc.dat",path);
416   FILE *infile = fopen(fname,"r");
417   if(!infile) return;
418   Int_t hitid,hitmc,i;
419   
420   for(i=0; ; i++)
421     if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
422   rewind(infile);
423   fNFastPoints = i;
424   fMcId = new Int_t[fNFastPoints];
425   fMcIndex = new UInt_t[fNFastPoints];
426   
427   for(i=0; i<fNFastPoints; i++)
428     {
429       if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
430       fMcId[i] = hitmc;
431       fMcIndex[i] = hitid;
432     }
433   fclose(infile);
434 }
435
436 void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
437 {
438   //Create the histograms 
439   
440   LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
441     <<"Creating histograms..."<<ENDLOG;
442   
443   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
444   fNtuppel->SetDirectory(0);
445   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.); 
446   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);    
447   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
448   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
449   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
450   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
451   
452   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
453   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
454   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
455   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
456   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
457
458 }
459
460 void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange)
461 {
462   //Read the good particles from file. This file should already have been
463   //generated by macro AliTPCComparison.C.
464   
465   Char_t filename[1024];
466   if(event<0 && !padrowrange)
467     sprintf(filename,"%s/good_tracks_tpc",path);
468   else if(event>=0 && !padrowrange)
469     sprintf(filename,"%s/good_tracks_tpc_%d",path,event);
470   else
471     sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]);
472   ifstream in(filename);
473   if(!in)
474     {
475       cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
476       return;
477     }
478   Int_t MaxTracks=20000;
479   if(fGoodTracks)
480     delete [] fGoodTracks;
481   fGoodGen=0;
482   fGoodTracks = new GoodTrack[MaxTracks];
483   
484   if(fStandardComparison){
485     while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
486            fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
487            fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
488       {
489         fGoodTracks[fGoodGen].nhits=-1;
490         fGoodTracks[fGoodGen].sector=-1; 
491         fGoodGen++;
492         if (fGoodGen==MaxTracks) 
493           {
494             cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
495             break;
496           }
497       }
498   } else {
499     while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
500            fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
501            fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector) 
502       {
503         fGoodGen++;
504         if (fGoodGen==MaxTracks) 
505           {
506             cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
507             break;
508           }
509       }
510   }
511 }
512
513 //has to be modified for fakes.
514
515 void AliL3Evaluate::FillEffHistos()
516 {  
517   if(!fGoodTracks)
518     {
519       cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
520       return;
521     }
522   //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
523   for(Int_t i=0; i<fGoodGen; i++)
524     {
525       //cout<<"Checking particle "<<i<<endl;
526       if(!fStandardComparison) 
527         if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
528       Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
529       if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
530       Float_t pzg=fGoodTracks[i].pz;
531       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
532       
533       //If we are only considering tracks on one side of the TPC:
534       if(fMaxSlice <= 17)
535         if(dipangle < 0)
536           continue;
537
538       fNGoodTracksPt->Fill(ptg);
539       fNGoodTracksEta->Fill(dipangle);
540       Int_t found = 0;
541       
542       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
543         {
544           AliL3Track *track = fTracks->GetCheckedTrack(k);
545           if(!track) continue;
546           Int_t nHits = track->GetNumberOfPoints();
547           if(nHits < fMinPointsOnTrack) break;
548           Int_t tracklabel;
549           tracklabel = track->GetMCid();
550           
551           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
552           found=1;
553           Float_t pt=track->GetPt();
554           if(tracklabel == fGoodTracks[i].label) 
555             {
556               fNFoundTracksPt->Fill(ptg); 
557               fNFoundTracksEta->Fill(dipangle);
558               fNtuppel->Fill(ptg,pt,nHits);
559               fPtRes->Fill((pt-ptg)/ptg*100.);
560             }
561           else 
562             {
563               fNFakeTracksPt->Fill(ptg); 
564               fNFakeTracksEta->Fill(dipangle);
565             }
566           //fPtRes->Fill((pt-ptg)/ptg*100.);
567           //fNtuppel->Fill(ptg,pt,nHits);
568           break;
569           
570         }
571       //if(!found)
572       //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
573     }
574 }
575
576 void AliL3Evaluate::FillEffHistosNAIVE()
577 {  
578   //Fill the efficiency histograms.
579   
580   cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
581   for(Int_t i=0; i<fGoodGen; i++)
582     {
583       if(!fStandardComparison) 
584         if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
585       Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
586       if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
587       Double_t pzg=fGoodTracks[i].pz;
588       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
589       //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
590       fNGoodTracksPt->Fill(ptg);
591       fNGoodTracksEta->Fill(dipangle);
592       
593     }
594   
595   for(Int_t k=0; k<fTracks->GetNTracks(); k++)
596     {
597       AliL3Track *track = fTracks->GetCheckedTrack(k);
598       if(!track) continue;
599       Int_t nHits = track->GetNumberOfPoints();
600       if(nHits < fMinPointsOnTrack) break;
601       if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue;
602       if(fabs(track->GetPseudoRapidity())>0.9) continue;
603
604       fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
605       //Float_t pt=track->GetPt();
606       //fPtRes->Fill((pt-ptg)/ptg*100.);
607       //fNtuppel->Fill(ptg,pt,nHits);
608             
609     }
610 }
611
612 void AliL3Evaluate::CalcEffHistos()
613 {  
614
615   Stat_t ngood=fNGoodTracksPt->GetEntries();
616   Stat_t nfound=fNFoundTracksPt->GetEntries();
617   Stat_t nfake=fNFakeTracksPt->GetEntries();
618
619   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
620     <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
621   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
622     <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
623   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
624     <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
625   LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
626     <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
627   //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
628   //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
629
630   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
631   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
632   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
633   fTrackEffPt->SetMaximum(1.4);
634   fTrackEffPt->SetXTitle("P_{T} [GeV]");
635   fTrackEffPt->SetLineWidth(2);
636   fFakeTrackEffPt->SetFillStyle(3013);
637   fTrackEffPt->SetLineColor(4);
638   fFakeTrackEffPt->SetFillColor(2);
639
640   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
641   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
642   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
643   fTrackEffEta->SetMaximum(1.4);
644   fTrackEffEta->SetXTitle("#lambda [degrees]");
645   fTrackEffEta->SetLineWidth(2);
646   fFakeTrackEffEta->SetFillStyle(3013);
647   fTrackEffEta->SetLineColor(4);
648   fFakeTrackEffEta->SetFillColor(2);
649        
650 }
651
652 void AliL3Evaluate::Write2File(Char_t *outputfile)
653 {
654   //Write histograms to file:
655   
656   TFile *of = TFile::Open(outputfile,"RECREATE");
657   if(!of->IsOpen())
658     {
659       LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
660         <<"Problems opening rootfile"<<ENDLOG;
661       return;
662     }
663   
664   of->cd();
665   fNtuppel->Write();
666   fPtRes->Write();
667   fNGoodTracksPt->Write();
668   fNFoundTracksPt->Write();
669   fNFakeTracksPt->Write();
670   fTrackEffPt->Write();
671   fFakeTrackEffPt->Write();
672   fNGoodTracksEta->Write();
673   fNFoundTracksEta->Write();
674   fNFakeTracksEta->Write();
675   fTrackEffEta->Write();
676   fFakeTrackEffEta->Write();
677   
678   of->Close();
679
680 }
681
682 TNtuple *AliL3Evaluate::GetNtuple()
683 {
684   if(!fNtupleRes)
685     {
686       fNtupleRes = new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
687       fNtupleRes->SetDirectory(0);//Bug in older version of root.
688     }
689   return fNtupleRes;
690 }
691
692 void AliL3Evaluate::CalculateResiduals()
693 {
694
695   TNtuple *ntuppel = GetNtuple();
696   
697   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
698     {
699       
700       AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
701       if(!track) continue;
702       if(track->GetNHits() < fMinPointsOnTrack) break;
703       
704       track->CalculateHelix();
705       UInt_t *hitnum = track->GetHitNumbers();
706       UInt_t id;
707       
708       Float_t xyz[3];
709       Int_t padrow;
710       for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
711         {
712           id = hitnum[j];
713           Int_t slice = (id>>25) & 0x7f;
714           Int_t patch = (id>>22) & 0x7;
715           UInt_t pos = id&0x3fffff;           
716
717           //if(slice<18) continue;
718           
719           AliL3SpacePointData *points = fClusters[slice][patch];
720           
721           if(!points) 
722             {
723               LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
724                 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
725               continue;
726             }
727           if(pos>=fNcl[slice][patch]) 
728             {
729               LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
730                 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
731               continue;
732             }
733           
734           xyz[0] = points[pos].fX;
735           xyz[1] = points[pos].fY;
736           xyz[2] = points[pos].fZ;
737           padrow = points[pos].fPadRow;
738           //AliL3Transform::Global2Local(xyz,slice,kTRUE);
739           AliL3Transform::Global2LocHLT(xyz,slice);
740           
741           Float_t angle = 0;
742           AliL3Transform::Local2GlobalAngle(&angle,slice);
743           if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
744             {
745               LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Crossing point")
746                 <<"Track does not crossing padrow "<<padrow<<" in slice "<<slice<<ENDLOG;
747               continue;
748             }
749           
750           Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
751           //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);          
752           AliL3Transform::Global2LocHLT(xyz_cross,slice);
753           
754           Double_t beta = track->GetCrossingAngle(padrow,slice);
755           
756           Double_t yres = xyz_cross[1] - xyz[1];
757           Double_t zres = xyz_cross[2] - xyz[2];
758           Double_t dipangle = atan(track->GetTgl());
759           ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
760           
761         }
762     }
763 }
764
765 enum tagprimary {kPrimaryCharged = 0x4000};
766 void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp)
767 {
768   //Compare points to the exact crossing points of track and padrows.
769   //The input file to this function, contains the exact clusters calculated
770   //in AliTPC::Hits2ExactClusters.
771   
772 #ifndef do_mc
773   cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl;
774   return;
775 #else
776   cout<<"Evaluating points"<<endl;
777   TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties",
778                                  "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2");
779   ntuppel->SetDirectory(0);
780   
781   TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen");
782   ntuppel2->SetDirectory(0);
783
784   TFile *exfile = TFile::Open(rootfile);
785   if(!exfile)
786     {
787       cerr<<"Error opening rootfile "<<rootfile<<endl;
788       return;
789     }
790   gAlice = (AliRun*)exfile->Get("gAlice");
791   if (!gAlice) 
792     {
793       LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
794         <<"AliRun object non existing on file"<<ENDLOG;
795       return;
796     }
797
798   AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
799   
800   TFile *exact = TFile::Open(exactfile);
801   if(!exact)
802     {
803       cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl;
804       return;
805     }
806   
807   AliStack *astack=gAlice->Stack();
808
809   AliTPCClustersArray *arr=0;
810   for(Int_t event=0; event<nevent; event++)
811     {
812       LoadData(event,sp);   
813       exfile->cd();
814       if(arr)
815         delete arr;
816       Int_t nparticles = gAlice->GetEvent(event);
817       Int_t nprimaries = 0;//FindPrimaries(nparticles);
818       cout<<"Event "<<event<<" had "<<nparticles<<" particles and "<<nprimaries<<" primaries"<<endl;
819       exact->cd();
820       
821       //Get the exact clusters from file:
822       AliTPCClustersArray *arr = new AliTPCClustersArray;
823       arr->Setup(param);
824       arr->SetClusterType("AliComplexCluster");
825       char treeName[500];
826       sprintf(treeName,"TreeCExact_%s_%d",param->GetTitle(),event);
827       Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
828       if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return;}
829       
830       //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
831       for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
832         {
833           //Get the exact clusters for this row:
834           Int_t cursec,currow;
835           AliSegmentID *s = arr->LoadEntry(i);
836           param->AdjustSectorRow(s->GetID(),cursec,currow);
837           
838           AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
839           TClonesArray *clusters = ro->GetArray();
840           int num_of_offline=clusters->GetEntriesFast();
841           
842           //Get the found clusters:
843           Int_t slice,padrow;
844           AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
845           if(slice < fMinSlice) continue;
846           if(slice > fMaxSlice) break;
847           
848           Int_t patch = AliL3Transform::GetPatch(padrow);
849           if(sp)
850             patch=0;
851           AliL3SpacePointData *points = fClusters[slice][patch];
852           if(!points)
853             continue;
854           
855           //cout<<"Slice "<<slice<<" padrow "<<padrow<<" has "<<num_of_offline<<" clusters "<<endl;
856           Int_t clustercount=0;
857           Int_t crosscount=0;
858           for(Int_t m=0; m<num_of_offline; m++)
859             {
860               AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
861 #ifdef use_newio
862               Int_t mcId = cluster->GetTrack(0);
863 #else
864               Int_t mcId = cluster->fTracks[0];
865 #endif        
866               if(mcId <0) continue;
867         
868 #ifdef use_newio      
869               if(cluster->GetY() < 1 || cluster->GetY() > AliL3Transform::GetNPads(padrow) - 2 ||
870                  cluster->GetX() < 1 || cluster->GetX() > AliL3Transform::GetNTimeBins() - 2)
871                 continue;
872 #else
873               if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
874                  cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
875                 continue;
876 #endif        
877               Float_t xyz_ex[3];
878               
879 #ifdef use_newio
880               AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX());
881 #else         
882               AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
883 #endif        
884               //In function AliTPC::Hits2ExactClusters the time offset is not included,
885               //so we have to substract it again here.
886               if(slice<18)
887                 xyz_ex[2]-=AliL3Transform::GetZOffset();
888               else
889                 xyz_ex[2]+=AliL3Transform::GetZOffset();
890               
891               //Outside our cone:
892               if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
893                 continue;
894               
895               TParticle *part = astack->Particle(mcId);
896               crosscount++;
897               
898               if(part->Pt() < fMinGoodPt) continue;
899               
900               //Dont take secondaries, because in width calculation we assume primaries:
901               //if(!(part->TestBit(kPrimaryCharged))) continue;
902               if(part->GetFirstMother()>=0) continue;
903               
904               Int_t tempcount=0;
905               for(UInt_t c=0; c<fNcl[slice][patch]; c++)
906                 {
907                   if((Int_t)points[c].fPadRow!=padrow) continue;
908                   Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
909                   
910                   if(!offline)
911                     AliL3Transform::Global2Local(xyz_cl,cursec);
912                   tempcount++;
913                   
914                   if(points[c].fTrackID[0] != mcId &&
915                      points[c].fTrackID[1] != mcId &&
916                      points[c].fTrackID[2] != mcId)
917                     continue;
918                   
919                   //Residuals:
920                   Float_t resy = xyz_cl[1] - xyz_ex[1];
921                   Float_t resz = xyz_cl[2] - xyz_ex[2];
922                   
923                   //Cluster shape
924                   Int_t charge = (Int_t)points[c].fCharge;
925                   Float_t beta = GetCrossingAngle(part,slice,padrow,xyz_ex);
926                   Double_t tanl = xyz_ex[2]/sqrt(xyz_ex[0]*xyz_ex[0]+xyz_ex[1]*xyz_ex[1]);
927                   Float_t psigmaY2 = AliL3Transform::GetParSigmaY2(padrow,xyz_ex[2],beta);
928                   Float_t psigmaZ2 = AliL3Transform::GetParSigmaZ2(padrow,xyz_ex[2],tanl);
929                   Float_t sigmaY2 = points[c].fSigmaY2;
930                   Float_t sigmaZ2 = points[c].fSigmaZ2;
931                   ntuppel->Fill(slice,padrow,charge,resy,resz,xyz_ex[2],part->Pt(),beta,sigmaY2,sigmaZ2,psigmaY2,psigmaZ2);
932                 }
933               clustercount=tempcount;
934             }
935           ntuppel2->Fill(slice,padrow,clustercount,crosscount);
936           arr->ClearRow(cursec,currow);
937         }
938     }
939   exfile->Close();
940   exact->Close();
941
942   TFile *ofile = TFile::Open(tofile,"RECREATE");
943   ntuppel->Write();
944   ntuppel2->Write();
945   ofile->Close();
946   
947 #endif
948 }
949
950 void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp)
951 {
952   //Evaluate the cluster finder efficiency.
953   
954 #ifndef do_mc
955   cerr<<"AliL3Evaluate::GetCFeff : Compile with do_mc flag"<<endl;
956   return;
957 #else
958   TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","slice:row:ncrossings:nclusters");
959   ntuppel->SetDirectory(0);
960   
961   Char_t filename[1024];
962   sprintf(filename,"%s/alirunfile.root",path);
963   TFile *rfile = TFile::Open(filename);
964   gAlice = (AliRun*)rfile->Get("gAlice");
965
966   AliStack *astack=gAlice->Stack();
967   
968   AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName());
969       
970   Int_t zero=param->GetZeroSup();
971
972   sprintf(filename,"%s/digitfile.root",path);
973   TFile *dfile = TFile::Open(filename);
974   
975   for(Int_t event=0; event<nevent; event++)
976     {
977       LoadData(event,sp);
978       rfile->cd();
979       gAlice->GetEvent(event);
980       Int_t np = astack->GetNtrack();
981       cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl;
982       dfile->cd();
983       sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event);
984       TTree *TD=(TTree*)gDirectory->Get(filename);
985       AliSimDigits da, *digits=&da;
986       TD->GetBranch("Segment")->SetAddress(&digits);
987       
988       Int_t crossed=0,recs=0;
989       Int_t *count = new Int_t[np]; //np number of particles.
990       Int_t i;
991       Float_t xyz[3];
992       for (i=0; i<np; i++) count[i]=0;
993       
994       
995       Int_t sec,row,sl,sr;
996       for(Int_t i=0; i<(Int_t)TD->GetEntries(); i++)
997         {
998           crossed=recs=0;
999           if (!TD->GetEvent(i)) continue;
1000           param->AdjustSectorRow(digits->GetID(),sec,row);
1001           AliL3Transform::Sector2Slice(sl,sr,sec,row);
1002           if(sl < fMinSlice) continue;
1003           if(sl > fMaxSlice) break;
1004           cout<<"Processing slice "<<sl<<" row "<<sr<<endl;
1005           digits->First();
1006           do {
1007             Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
1008             Short_t dig = digits->GetDigit(it,ip);
1009             
1010             if(dig<=param->GetZeroSup()) continue;
1011             AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
1012             if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
1013               continue;
1014             
1015             Int_t idx0=digits->GetTrackID(it,ip,0); 
1016             Int_t idx1=digits->GetTrackID(it,ip,1);
1017             Int_t idx2=digits->GetTrackID(it,ip,2);
1018             
1019             if (idx0>=0 && dig>=zero) count[idx0]+=1;
1020             if (idx1>=0 && dig>=zero) count[idx1]+=1;
1021             if (idx2>=0 && dig>=zero) count[idx2]+=1;
1022           } while (digits->Next());
1023           for (Int_t j=0; j<np; j++) 
1024             {
1025               TParticle *part = astack->Particle(j);
1026               if(part->Pt() < fMinGoodPt) continue;
1027               if(part->GetFirstMother() >= 0) continue;
1028               if (count[j]>1) //at least two digits at this padrow 
1029                 {
1030                   crossed++;
1031                   count[j]=0;
1032                 }
1033             }
1034
1035           Int_t patch = AliL3Transform::GetPatch(sr);
1036           if(sp==kTRUE)
1037             patch=0;
1038           AliL3SpacePointData *points = fClusters[sl][patch];
1039           if(!points)
1040             continue;
1041           for(UInt_t k=0; k<fNcl[sl][patch]; k++)
1042             {
1043               if(points[k].fPadRow!=sr) continue;
1044               recs++;
1045             }
1046           ntuppel->Fill(sl,sr,crossed,recs);
1047         }
1048       
1049       TD->Delete();
1050       delete[] count;
1051     }
1052   TFile *file = TFile::Open(outfile,"RECREATE");
1053   ntuppel->Write();
1054   file->Close();
1055   
1056   rfile->Close();
1057   dfile->Close();
1058 #endif
1059 }
1060
1061 Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz)
1062 {
1063   //Calculate the padrow crossing angle of the particle
1064   
1065   Double_t kappa = AliL3Transform::GetBField()*AliL3Transform::GetBFact()/part->Pt();
1066   
1067   Double_t radius = 1/fabs(kappa);
1068   if(part->GetPdgCode() > 0) kappa = -kappa;
1069
1070   Float_t angl[1] = {part->Phi()};
1071   
1072   AliL3Transform::Global2LocalAngle(angl,slice);
1073   
1074   Double_t charge = -1.*kappa;
1075
1076   Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge);
1077   
1078   Double_t x0=0;
1079   Double_t y0=0;
1080   Double_t xc = x0 - radius * cos(trackPhi0);
1081   Double_t yc = y0 - radius * sin(trackPhi0);
1082
1083   Double_t tangent[2];
1084   tangent[0] = -1.*(xyz[1] - yc)/radius;
1085   tangent[1] = (xyz[0] - xc)/radius;
1086   
1087   Double_t perp_padrow[2] = {1,0}; //locally in slice
1088   
1089   Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
1090   if(cos_beta > 1) cos_beta=1;
1091   return acos(cos_beta);
1092 }
1093
1094 Int_t AliL3Evaluate::FindPrimaries(Int_t nparticles)
1095 {
1096   // cuts:
1097   Double_t vertcut = 0.001;
1098   Double_t decacut = 3.;
1099   Double_t timecut = 0.;
1100   Int_t nprch1=0;
1101   AliStack *astack=gAlice->Stack();
1102   TParticle * part = astack->Particle(0);
1103   Double_t xori = part->Vx();
1104   Double_t yori = part->Vy();
1105   Double_t zori = part->Vz();
1106   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
1107     
1108     part = astack->Particle(iprim);
1109     const char * xxx=strstr(part->GetName(),"XXX");
1110     if(xxx)continue;
1111     
1112     TParticlePDG *ppdg = part->GetPDG();
1113     if(TMath::Abs(ppdg->Charge())!=3)continue;  // only charged (no quarks)
1114     
1115     Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
1116     if(dist>vertcut)continue;  // cut on the vertex
1117     
1118     if(part->T()>timecut)continue;
1119     
1120     Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
1121     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
1122     
1123     Bool_t prmch = kTRUE;   // candidate primary track
1124     Int_t fidau=part->GetFirstDaughter();  // cut on daughters
1125     Int_t lasdau=0;
1126     Int_t ndau=0;
1127     if(fidau>=0){
1128       lasdau=part->GetLastDaughter();
1129       ndau=lasdau-fidau+1;
1130     }
1131     if(ndau>0){
1132       for(Int_t j=fidau;j<=lasdau;j++){
1133         TParticle *dau=astack->Particle(j);
1134         Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
1135         if(distd<decacut)prmch=kFALSE;  // eliminate if the decay is near the vertex
1136       }
1137     }
1138     
1139     if(prmch){
1140       nprch1++;
1141       part->SetBit(kPrimaryCharged);
1142     }
1143   }
1144
1145   return nprch1;
1146 }