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