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