Fix for memory leak in fClusters
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
CommitLineData
3e87ef69 1// @(#) $Id$
e8541476 2
b661165c 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
3e87ef69 4//*-- Copyright &copy ALICE HLT Group
108615fc 5
118c26c3 6#include "AliL3StandardIncludes.h"
108615fc 7#include <TFile.h>
8#include <TH1.h>
9#include <TParticle.h>
10#include <TTree.h>
11#include <TClonesArray.h>
12
4499ed26 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>
b2a02bce 21#include <AliStack.h>
108615fc 22
0bd0c1ef 23#if __GNUC__ == 3
118c26c3 24#include <fstream>
25#include <iosfwd>
26#else
27#include <fstream.h>
28#endif
29
30#include "AliL3Logging.h"
108615fc 31#include "AliL3Transform.h"
32#include "AliL3SpacePointData.h"
33#include "AliL3Track.h"
34#include "AliL3FileHandler.h"
35#include "AliL3TrackArray.h"
36#include "AliL3Evaluate.h"
118c26c3 37
0bd0c1ef 38#if __GNUC__ == 3
118c26c3 39using namespace std;
40#endif
108615fc 41
3e87ef69 42/** \class AliL3Evaluate
43<pre>
b661165c 44//_____________________________________________________________
45// AliL3Evaluate
46//
47// Evaluation class for tracking. Plots efficiencies etc..
84122e15 48//
3e87ef69 49</pre>
50*/
108615fc 51
52ClassImp(AliL3Evaluate)
53
54AliL3Evaluate::AliL3Evaluate()
55{
e0f4d6b2 56 fTracks = 0;
c0217812 57 fNFastPoints = 0;
58 fMcIndex = 0;
59 fMcId = 0;
1f79afc0 60 fMinSlice=0;
61 fMaxSlice=0;
e0f4d6b2 62 fGoodTracks=0;
b2a02bce 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;
108615fc 80}
81
b2a02bce 82AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice)
108615fc 83{
3e87ef69 84
e0f4d6b2 85 if(slice)
86 {
87 fMinSlice=slice[0];
88 fMaxSlice=slice[1];
89 }
90 else
91 {
92 fMinSlice=0;
93 fMaxSlice=35;
94 }
3e87ef69 95 sprintf(fPath,"%s",datapath);
e0f4d6b2 96 fMaxFalseClusters = 0.1;
23908b9b 97 fGoodFound = 0;
98 fGoodGen = 0;
e0f4d6b2 99 fMinPointsOnTrack = min_clusters;
100 fMinHitsFromParticle = minhits;
101 fMinGoodPt = minpt;
b2a02bce 102 fMaxGoodPt = maxpt;
3e87ef69 103 memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
104 fTracks=0;
105 fGoodTracks=0;
b2a02bce 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;
3e87ef69 123}
124
125void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
126{
e0f4d6b2 127 Char_t fname[1024];
40896357 128 AliL3FileHandler *clusterfile[36][6];
108615fc 129 for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
130 {
e0f4d6b2 131 for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
108615fc 132 {
3e87ef69 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];
1f79afc0 141 fClusters[s][p] = 0;
108615fc 142 clusterfile[s][p] = new AliL3FileHandler();
3e87ef69 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);
108615fc 147 if(!clusterfile[s][p]->SetBinaryInput(fname))
148 {
149 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
150 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
bc2f4f0e 151 delete clusterfile[s][p];
152 clusterfile[s][p] = 0;
153 continue;
108615fc 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();
3e87ef69 158 if(sp==kTRUE)
159 break;
108615fc 160 }
161 }
b2a02bce 162
3e87ef69 163 sprintf(fname,"%s/tracks_%d.raw",fPath,event);
108615fc 164 AliL3FileHandler *tfile = new AliL3FileHandler();
e0f4d6b2 165 if(!tfile->SetBinaryInput(fname)){
bc2f4f0e 166 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
e0f4d6b2 167 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
bc2f4f0e 168 return;
169 }
3e87ef69 170 if(fTracks)
171 delete fTracks;
108615fc 172 fTracks = new AliL3TrackArray();
173 tfile->Binary2TrackArray(fTracks);
174 tfile->CloseBinaryInput();
175 delete tfile;
b2a02bce 176 fTracks->QSort();
bc2f4f0e 177}
178
108615fc 179
e0f4d6b2 180AliL3Evaluate::~AliL3Evaluate()
108615fc 181{
e0f4d6b2 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;
b2a02bce 199 if(fNtupleRes) delete fNtupleRes;
37a5a5a9 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];
108615fc 203}
204
de3c3890 205void 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}
8af8a456 222
108615fc 223void AliL3Evaluate::AssignIDs()
224{
225 //Assign MC id to the tracks.
e0f4d6b2 226#ifndef do_mc
227 cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
228 return;
229#endif
3e87ef69 230 fGoodFound=0;
108615fc 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);
c0217812 237 if(!track) continue;
108615fc 238 if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
23908b9b 239
23908b9b 240 fGoodFound++;
c0217812 241 Int_t tID = GetMCTrackLabel(track);
108615fc 242 track->SetMCid(tID);
108615fc 243 }
3e87ef69 244 //cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
108615fc 245}
246
de3c3890 247Float_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
37a5a5a9 298 // cout<<" PID: "<<nc<<" "<<nl<<" "<<nu<<" "<<trackDEdx<<" "<<track->GetPt()<<endl;
de3c3890 299 return trackDEdx;
300}
108615fc 301
302struct S {Int_t lab; Int_t max;};
c0217812 303Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
108615fc 304 //Returns the MCtrackID of the belonging clusters.
305 //If MCLabel < 0, means that track is fake.
3e87ef69 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
494fad94 311
312#ifdef do_mc
108615fc 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;
108615fc 319
108615fc 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];
c0217812 330 if(!points) continue;
108615fc 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:
3e87ef69 339 lab=points[pos].fTrackID[0];
340
108615fc 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
3e87ef69 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
108615fc 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];
c0217812 367 if(!points) continue;
108615fc 368 if(pos>=fNcl[slice][patch])
369 {
370 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
371 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
372 continue;
373 }
374
355debe1 375 if (abs(points[pos].fTrackID[1]) == lab ||
376 abs(points[pos].fTrackID[2]) == lab ) max++;
108615fc 377 }
378
3e87ef69 379
380 //Check if more than 10% of the clusters were assigned incorrectly:
e0f4d6b2 381 if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters)
108615fc 382 {
383 return -lab;
384 }
3e87ef69 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
108615fc 404 return lab;
494fad94 405#else //If we are running with mc_ids or not
406 return 0;
407#endif
408
108615fc 409}
410
c0217812 411void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
108615fc 412{
413 //Get the MC id of space points in case of using the fast simulator.
c0217812 414 char fname[256];
415 sprintf(fname,"%s/point_mc.dat",path);
416 FILE *infile = fopen(fname,"r");
417 if(!infile) return;
108615fc 418 Int_t hitid,hitmc,i;
419
420 for(i=0; ; i++)
421 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
108615fc 422 rewind(infile);
c0217812 423 fNFastPoints = i;
424 fMcId = new Int_t[fNFastPoints];
425 fMcIndex = new UInt_t[fNFastPoints];
108615fc 426
c0217812 427 for(i=0; i<fNFastPoints; i++)
108615fc 428 {
429 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
c0217812 430 fMcId[i] = hitmc;
431 fMcIndex[i] = hitid;
108615fc 432 }
433 fclose(infile);
c0217812 434}
435
e0f4d6b2 436void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
c0217812 437{
438 //Create the histograms
439
40896357 440 LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
441 <<"Creating histograms..."<<ENDLOG;
442
1f79afc0 443 fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
a8a3f5c2 444 fNtuppel->SetDirectory(0);
c0217812 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);
e0f4d6b2 457
458}
459
3e87ef69 460void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange)
e0f4d6b2 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];
3e87ef69 466 if(event<0 && !padrowrange)
a8a3f5c2 467 sprintf(filename,"%s/good_tracks_tpc",path);
3e87ef69 468 else if(event>=0 && !padrowrange)
469 sprintf(filename,"%s/good_tracks_tpc_%d",path,event);
a8a3f5c2 470 else
3e87ef69 471 sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]);
e0f4d6b2 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;
3e87ef69 479 if(fGoodTracks)
480 delete [] fGoodTracks;
481 fGoodGen=0;
e0f4d6b2 482 fGoodTracks = new GoodTrack[MaxTracks];
483
b2a02bce 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 }
c0217812 511}
eeddc64d 512
b2a02bce 513//has to be modified for fakes.
a8a3f5c2 514
e1a32fa5 515void AliL3Evaluate::FillEffHistos()
c0217812 516{
e0f4d6b2 517 if(!fGoodTracks)
518 {
519 cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
520 return;
521 }
3e87ef69 522 //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
e1a32fa5 523 for(Int_t i=0; i<fGoodGen; i++)
c0217812 524 {
e0f4d6b2 525 //cout<<"Checking particle "<<i<<endl;
b2a02bce 526 if(!fStandardComparison)
527 if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
e0f4d6b2 528 Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
b2a02bce 529 if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
e0f4d6b2 530 Float_t pzg=fGoodTracks[i].pz;
c0217812 531 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
b2a02bce 532
533 //If we are only considering tracks on one side of the TPC:
534 if(fMaxSlice <= 17)
535 if(dipangle < 0)
536 continue;
e0f4d6b2 537
c0217812 538 fNGoodTracksPt->Fill(ptg);
539 fNGoodTracksEta->Fill(dipangle);
540 Int_t found = 0;
e0f4d6b2 541
c0217812 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;
c0217812 548 Int_t tracklabel;
549 tracklabel = track->GetMCid();
e0f4d6b2 550
e1a32fa5 551 if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
c0217812 552 found=1;
c0217812 553 Float_t pt=track->GetPt();
3e87ef69 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);
c0217812 568 break;
569
570 }
3e87ef69 571 //if(!found)
b2a02bce 572 //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
c0217812 573 }
574}
eeddc64d 575
e0f4d6b2 576void AliL3Evaluate::FillEffHistosNAIVE()
355debe1 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 {
b2a02bce 583 if(!fStandardComparison)
584 if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
e0f4d6b2 585 Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
b2a02bce 586 if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
e0f4d6b2 587 Double_t pzg=fGoodTracks[i].pz;
355debe1 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;
b2a02bce 601 if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue;
355debe1 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}
c0217812 611
3e87ef69 612void AliL3Evaluate::CalcEffHistos()
613{
614
c0217812 615 Stat_t ngood=fNGoodTracksPt->GetEntries();
616 Stat_t nfound=fNFoundTracksPt->GetEntries();
617 Stat_t nfake=fNFakeTracksPt->GetEntries();
3e87ef69 618
c0217812 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;
3e87ef69 627 //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
628 //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
355debe1 629
c0217812 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);
3e87ef69 639
c0217812 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
108615fc 650}
651
652void AliL3Evaluate::Write2File(Char_t *outputfile)
653{
654 //Write histograms to file:
655
e0f4d6b2 656 TFile *of = TFile::Open(outputfile,"RECREATE");
108615fc 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();
1f79afc0 665 fNtuppel->Write();
108615fc 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();
108615fc 679
680}
681
b2a02bce 682TNtuple *AliL3Evaluate::GetNtuple()
1f79afc0 683{
b2a02bce 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}
1f79afc0 691
b2a02bce 692void AliL3Evaluate::CalculateResiduals()
693{
694
695 TNtuple *ntuppel = GetNtuple();
3e87ef69 696
697 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
1f79afc0 698 {
3e87ef69 699
700 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
701 if(!track) continue;
b2a02bce 702 if(track->GetNHits() < fMinPointsOnTrack) break;
3e87ef69 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++)
1f79afc0 711 {
3e87ef69 712 id = hitnum[j];
713 Int_t slice = (id>>25) & 0x7f;
714 Int_t patch = (id>>22) & 0x7;
715 UInt_t pos = id&0x3fffff;
b2a02bce 716
717 //if(slice<18) continue;
1f79afc0 718
3e87ef69 719 AliL3SpacePointData *points = fClusters[slice][patch];
1f79afc0 720
3e87ef69 721 if(!points)
1f79afc0 722 {
3e87ef69 723 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
724 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
725 continue;
1f79afc0 726 }
3e87ef69 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;
5a31e9df 738 //AliL3Transform::Global2Local(xyz,slice,kTRUE);
739 AliL3Transform::Global2LocHLT(xyz,slice);
3e87ef69 740
741 Float_t angle = 0;
742 AliL3Transform::Local2GlobalAngle(&angle,slice);
b2a02bce 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
3e87ef69 750 Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
5a31e9df 751 //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);
752 AliL3Transform::Global2LocHLT(xyz_cross,slice);
3e87ef69 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];
3e87ef69 758 Double_t dipangle = atan(track->GetTgl());
759 ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
760
1f79afc0 761 }
1f79afc0 762 }
1f79afc0 763}
764
3e87ef69 765enum tagprimary {kPrimaryCharged = 0x4000};
766void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp)
78127e35 767{
e1a32fa5 768 //Compare points to the exact crossing points of track and padrows.
1b76fd18 769 //The input file to this function, contains the exact clusters calculated
770 //in AliTPC::Hits2ExactClusters.
8af8a456 771
5bfdc6ad 772#ifndef do_mc
773 cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl;
3e87ef69 774 return;
5bfdc6ad 775#else
355debe1 776 cout<<"Evaluating points"<<endl;
3e87ef69 777 TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties",
778 "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2");
8af8a456 779 ntuppel->SetDirectory(0);
78127e35 780
3e87ef69 781 TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen");
782 ntuppel2->SetDirectory(0);
783
784 TFile *exfile = TFile::Open(rootfile);
355debe1 785 if(!exfile)
786 {
3e87ef69 787 cerr<<"Error opening rootfile "<<rootfile<<endl;
788 return;
355debe1 789 }
8af8a456 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;
3e87ef69 795 return;
8af8a456 796 }
3e87ef69 797
5e0f9911 798 AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
e1a32fa5 799
3e87ef69 800 TFile *exact = TFile::Open(exactfile);
801 if(!exact)
78127e35 802 {
3e87ef69 803 cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl;
804 return;
805 }
806
b2a02bce 807 AliStack *astack=gAlice->Stack();
808
3e87ef69 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();
78127e35 820
3e87ef69 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;}
78127e35 829
3e87ef69 830 //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
831 for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
8af8a456 832 {
3e87ef69 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;
e1a32fa5 858 for(Int_t m=0; m<num_of_offline; m++)
78127e35 859 {
e1a32fa5 860 AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
b2a02bce 861#ifdef use_newio
862 Int_t mcId = cluster->GetTrack(0);
863#else
e1a32fa5 864 Int_t mcId = cluster->fTracks[0];
b2a02bce 865#endif
294ee877 866 if(mcId <0) continue;
b2a02bce 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
3e87ef69 873 if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
874 cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
e1a32fa5 875 continue;
b2a02bce 876#endif
3e87ef69 877 Float_t xyz_ex[3];
878
b2a02bce 879#ifdef use_newio
880 AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX());
881#else
494fad94 882 AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
b2a02bce 883#endif
8af8a456 884 //In function AliTPC::Hits2ExactClusters the time offset is not included,
885 //so we have to substract it again here.
3e87ef69 886 if(slice<18)
887 xyz_ex[2]-=AliL3Transform::GetZOffset();
888 else
889 xyz_ex[2]+=AliL3Transform::GetZOffset();
8af8a456 890
3e87ef69 891 //Outside our cone:
892 if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
893 continue;
894
b2a02bce 895 TParticle *part = astack->Particle(mcId);
3e87ef69 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;
78127e35 934 }
3e87ef69 935 ntuppel2->Fill(slice,padrow,clustercount,crosscount);
936 arr->ClearRow(cursec,currow);
937 }
78127e35 938 }
3e87ef69 939 exfile->Close();
940 exact->Close();
e1a32fa5 941
3e87ef69 942 TFile *ofile = TFile::Open(tofile,"RECREATE");
943 ntuppel->Write();
944 ntuppel2->Write();
945 ofile->Close();
946
5bfdc6ad 947#endif
78127e35 948}
355debe1 949
3e87ef69 950void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp)
78127e35 951{
3e87ef69 952 //Evaluate the cluster finder efficiency.
78127e35 953
3e87ef69 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);
355debe1 960
3e87ef69 961 Char_t filename[1024];
962 sprintf(filename,"%s/alirunfile.root",path);
963 TFile *rfile = TFile::Open(filename);
964 gAlice = (AliRun*)rfile->Get("gAlice");
b2a02bce 965
966 AliStack *astack=gAlice->Stack();
355debe1 967
3e87ef69 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);
355debe1 974
3e87ef69 975 for(Int_t event=0; event<nevent; event++)
78127e35 976 {
3e87ef69 977 LoadData(event,sp);
978 rfile->cd();
979 gAlice->GetEvent(event);
b2a02bce 980 Int_t np = astack->GetNtrack();
3e87ef69 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++)
78127e35 997 {
3e87ef69 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();
355debe1 1006 do {
3e87ef69 1007 Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
1008 Short_t dig = digits->GetDigit(it,ip);
355debe1 1009
3e87ef69 1010 if(dig<=param->GetZeroSup()) continue;
494fad94 1011 AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
3e87ef69 1012 if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
355debe1 1013 continue;
1014
3e87ef69 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);
355debe1 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;
3e87ef69 1022 } while (digits->Next());
355debe1 1023 for (Int_t j=0; j<np; j++)
78127e35 1024 {
b2a02bce 1025 TParticle *part = astack->Particle(j);
3e87ef69 1026 if(part->Pt() < fMinGoodPt) continue;
1027 if(part->GetFirstMother() >= 0) continue;
355debe1 1028 if (count[j]>1) //at least two digits at this padrow
494fad94 1029 {
1030 crossed++;
1031 count[j]=0;
1032 }
355debe1 1033 }
3e87ef69 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++)
355debe1 1042 {
3e87ef69 1043 if(points[k].fPadRow!=sr) continue;
355debe1 1044 recs++;
78127e35 1045 }
3e87ef69 1046 ntuppel->Fill(sl,sr,crossed,recs);
78127e35 1047 }
494fad94 1048
3e87ef69 1049 TD->Delete();
1050 delete[] count;
78127e35 1051 }
355debe1 1052 TFile *file = TFile::Open(outfile,"RECREATE");
1053 ntuppel->Write();
1054 file->Close();
1055
3e87ef69 1056 rfile->Close();
1057 dfile->Close();
1058#endif
1059}
1060
dd7d3870 1061Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz)
3e87ef69 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);
78127e35 1092}
355debe1 1093
3e87ef69 1094Int_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;
b2a02bce 1101 AliStack *astack=gAlice->Stack();
1102 TParticle * part = astack->Particle(0);
3e87ef69 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
b2a02bce 1108 part = astack->Particle(iprim);
62bb4b3d 1109 const char * xxx=strstr(part->GetName(),"XXX");
3e87ef69 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++){
b2a02bce 1133 TParticle *dau=astack->Particle(j);
3e87ef69 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}