]> 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
Optimization
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
CommitLineData
108615fc 1//Author: Anders Strand Vestbo
2//Last Modified: 5.01.2001
3
4#include <stdio.h>
e170e0f9 5#include <math.h>
108615fc 6#include <TFile.h>
7#include <TH1.h>
8#include <TParticle.h>
9#include <TTree.h>
10#include <TClonesArray.h>
11
12#include "AliRun.h"
13#include "AliSimDigits.h"
14#include "AliTPC.h"
15#include "AliTPCClustersArray.h"
16#include "AliTPCClustersRow.h"
17#include "AliTPCcluster.h"
18#include "AliTPCParam.h"
e1a32fa5 19#include "AliComplexCluster.h"
108615fc 20
1f79afc0 21#include "AliL3Defs.h"
108615fc 22#include "AliL3Transform.h"
23#include "AliL3SpacePointData.h"
24#include "AliL3Track.h"
25#include "AliL3FileHandler.h"
26#include "AliL3TrackArray.h"
27#include "AliL3Evaluate.h"
28#include "AliL3Logging.h"
29
30//AliL3Evaluate
31//Class for tracking evaluation.
32
33ClassImp(AliL3Evaluate)
34
35AliL3Evaluate::AliL3Evaluate()
36{
e1a32fa5 37 fDigitsFile = NULL;
108615fc 38 fTracks = NULL;
39 fMCclusterfile = NULL;
c0217812 40 fNFastPoints = 0;
41 fMcIndex = 0;
42 fMcId = 0;
1f79afc0 43 fMinSlice=0;
44 fMaxSlice=0;
23908b9b 45 fTransform = new AliL3Transform();
108615fc 46}
47
e1a32fa5 48AliL3Evaluate::AliL3Evaluate(Char_t *mcfile,Char_t *digitsfile,Int_t *slice)
108615fc 49{
50 //Normal constructor. Input are the rootfile containing the
51 //original MC information of the simulated event.
52
e1a32fa5 53 fDigitsFile = new TFile(digitsfile,"READ");
54 fEventFile = new TFile(mcfile,"READ");
55 if(!fDigitsFile->IsOpen())
108615fc 56 {
57 LOG(AliL3Log::kError,"AliL3Evaluation::AliL3Evaluation","File Open")
58 <<"Inputfile "<<mcfile<<" does not exist"<<ENDLOG;
59 return;
60 }
61
e1a32fa5 62 fParam = (AliTPCParam*)fEventFile->Get("75x40_100x60");
108615fc 63 fTransform = new AliL3Transform();
64
65 fMinSlice = slice[0];
66 fMaxSlice = slice[1];
e170e0f9 67 fMinGoodPt = 0.1;
23908b9b 68 fNoOverlap = kFALSE;
108615fc 69}
70
1f79afc0 71AliL3Evaluate::AliL3Evaluate(Int_t *slice)
72{
78127e35 73 //ctor to use if you do not need any rootfile.
74
1f79afc0 75
76 fMinSlice = slice[0];
77 fMaxSlice = slice[1];
78 fTransform = new AliL3Transform();
79
80}
81
108615fc 82AliL3Evaluate::~AliL3Evaluate()
83{
c0217812 84 if(fDigitsTree) fDigitsTree->Delete();
e1a32fa5 85 if(fDigitsFile) {
86 fDigitsFile->Close();
87 delete fDigitsFile;
c0217812 88 }
89 if(fTransform) delete fTransform;
90 if(fTracks) delete fTracks;
91 if(fPtRes) delete fPtRes;
92 if(fNGoodTracksPt) delete fNGoodTracksPt;
93 if(fNFoundTracksPt) delete fNFoundTracksPt;
94 if(fNFakeTracksPt) delete fNFakeTracksPt;
95 if(fTrackEffPt) delete fTrackEffPt;
96 if(fFakeTrackEffPt) delete fFakeTrackEffPt;
97 if(fNGoodTracksEta) delete fNGoodTracksEta;
98 if(fNFoundTracksEta) delete fNFoundTracksEta;
99 if(fNFakeTracksEta) delete fNFakeTracksEta;
100 if(fTrackEffEta) delete fTrackEffEta;
101 if(fFakeTrackEffEta) delete fFakeTrackEffEta;
102 if(fMcIndex) delete [] fMcIndex;
103 if(fMcId) delete [] fMcId;
1f79afc0 104 if(fNtuppel) delete fNtuppel;
108615fc 105}
106
bc2f4f0e 107void AliL3Evaluate::Setup(Char_t *trackfile,Char_t *path)
108615fc 108{
108615fc 109 //Read in the hit and track information from produced files.
110
23908b9b 111 fGoodFound = 0;
112 fGoodGen = 0;
108615fc 113 Char_t fname[256];
114 AliL3FileHandler *clusterfile[36][5];
115 for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
116 {
117 for(Int_t p=0; p<5; p++)
118 {
1f79afc0 119 fClusters[s][p] = 0;
108615fc 120 clusterfile[s][p] = new AliL3FileHandler();
bc2f4f0e 121 sprintf(fname,"%s/points_%d_%d.raw",path,s,p);
108615fc 122 if(!clusterfile[s][p]->SetBinaryInput(fname))
123 {
124 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
125 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
bc2f4f0e 126 delete clusterfile[s][p];
127 clusterfile[s][p] = 0;
128 continue;
108615fc 129 }
130 fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
131 clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
132 clusterfile[s][p]->CloseBinaryInput();
133 }
134 }
78127e35 135
23908b9b 136
108615fc 137 AliL3FileHandler *tfile = new AliL3FileHandler();
bc2f4f0e 138 if(!tfile->SetBinaryInput(trackfile)){
139 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
140 <<"Inputfile "<<trackfile<<" does not exist"<<ENDLOG;
141 return;
142 }
108615fc 143 fTracks = new AliL3TrackArray();
144 tfile->Binary2TrackArray(fTracks);
145 tfile->CloseBinaryInput();
146 delete tfile;
bc2f4f0e 147}
148
149void AliL3Evaluate::SetupSlow(Char_t *trackfile,Char_t *path)
150{
151 //Setup for using the slow simulator.
152
153 fIsSlow = true;
154 Setup(trackfile,path);
108615fc 155
e1a32fa5 156 if(!InitMC())
bc2f4f0e 157 LOG(AliL3Log::kError,"AliL3Evaluation::SetupSlow","Digits Tree")
158 <<"Error setting up digits tree"<<ENDLOG;
e1a32fa5 159
108615fc 160}
161
bc2f4f0e 162void AliL3Evaluate::SetupFast(Char_t *trackfile,Char_t *mcClusterfile,Char_t *path)
108615fc 163{
164 //Setup for using the fast simulator.
165
166 fIsSlow = false;
c0217812 167 GetFastClusterIDs(path);
108615fc 168
169 fMCclusterfile = new TFile(mcClusterfile);
170 if(!fMCclusterfile->IsOpen())
171 LOG(AliL3Log::kError,"AliL3Evaluation::SetupFast","File Open")
172 <<"Inputfile "<<mcClusterfile<<" does not exist"<<ENDLOG;
173
bc2f4f0e 174 Setup(trackfile,path);
e1a32fa5 175 InitMC();
108615fc 176}
177
e1a32fa5 178Bool_t AliL3Evaluate::InitMC()
108615fc 179{
e1a32fa5 180 if(fIsSlow)
108615fc 181 {
e1a32fa5 182 fDigitsFile->cd();
183 fDigitsTree = (TTree*)fDigitsFile->Get("TreeD_75x40_100x60");
184 if(!fDigitsTree) return false;
185 fDigitsTree->GetBranch("Segment")->SetAddress(&fDigits);
186 for(Int_t i=0; i<fDigitsTree->GetEntries(); i++)
187 {
188 if(!fDigitsTree->GetEvent(i)) continue;
189 Int_t se,ro,slice,slicerow;
190 fParam->AdjustSectorRow(fDigits->GetID(),se,ro);
191 fTransform->Sector2Slice(slice,slicerow,se,ro);
192 fRowid[slice][slicerow] = i;
193 }
108615fc 194 }
195
e1a32fa5 196 fEventFile->cd();
197 AliRun *gAlice = (AliRun*)fEventFile->Get("gAlice");
108615fc 198 if (!gAlice)
199 {
200 LOG(AliL3Log::kError,"AliL3Evaluate::SetParticleArray","gAlice")
201 <<"AliRun object non existing on file"<<ENDLOG;
202 return false;
203 }
204
205 gAlice->GetEvent(0);
e1a32fa5 206
108615fc 207 return true;
e1a32fa5 208
108615fc 209}
210
211
e1a32fa5 212
213void AliL3Evaluate::DefineGoodTracks(Int_t *slice,Int_t *padrow,Int_t good_number)
108615fc 214{
215 //Loop over MC particles, and mark the good ones
216 //(which the tracker should find...)
e1a32fa5 217
108615fc 218 AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
e1a32fa5 219
108615fc 220 TPC->SetParam(fParam);
e1a32fa5 221
108615fc 222 Int_t ver = TPC->IsVersion();
223 LOG(AliL3Log::kInformational,"AliL3Evaluate::DefineGoodTracks","TPC version")
224 <<"TPC version "<<ver<<" found on file"<<ENDLOG;
225
108615fc 226 Int_t zero=TPC->GetParam()->GetZeroSup();
227
e1a32fa5 228 Int_t np = gAlice->GetNtrack();
108615fc 229 Int_t *good = new Int_t[np];
230 for(Int_t ii=0; ii<np; ii++)
231 good[ii] = 0;
232
233 if(ver==1)
234 {
235 if(fIsSlow)
236 LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
237 <<"TPC version "<<ver<<" does not match."<<ENDLOG;
238 fMCclusterfile->cd();
239 AliTPCClustersArray carray;
240 carray.Setup(fParam);
241 carray.SetClusterType("AliTPCcluster");
242 Bool_t clusterok = carray.ConnectTree("Segment Tree");
243 if(!clusterok)
244 LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
245 <<"Error loading clusters from rootfile"<<ENDLOG;
78127e35 246
bc2f4f0e 247 for(Int_t i=0; i<carray.GetTree()->GetEntries(); i++)
108615fc 248 {
249 Int_t sec,row,sl,lr;
250 AliSegmentID *s = carray.LoadEntry(i);
251 fParam->AdjustSectorRow(s->GetID(),sec,row);
252 fTransform->Sector2Slice(sl,lr,sec,row);
253
e1a32fa5 254 if(sl != slice[0]) {carray.ClearRow(sec,row); continue;}
bc2f4f0e 255 if(lr < padrow[0]) {carray.ClearRow(sec,row); continue;}
256 if(lr > padrow[1]) {carray.ClearRow(sec,row); continue;}
108615fc 257 AliTPCClustersRow *cRow = carray.GetRow(sec,row);
258 for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
259 {
260 AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
261 Int_t lab=cluster->GetLabel(0);
262 if(lab<0) continue;
263 lab=TMath::Abs(lab);
264 good[lab]++;
265 }
266 if(carray.GetRow(sec,row))
267 carray.ClearRow(sec,row);
268 }
269 }
270 else if(ver==2)
271 {
272 if(!fIsSlow)
273 LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","TPC version")
274 <<"TPC version "<<ver<<" does not match."<<ENDLOG;
275 Int_t *count = new Int_t[np]; //np number of particles.
276 Int_t i;
277 for (i=0; i<np; i++) count[i]=0;
e1a32fa5 278 for(Int_t sl=slice[0]; sl<=slice[1]; sl++)
279 {
280 for (i=padrow[0]; i<=padrow[1]; i++) {
281 Int_t index = fRowid[sl][i];
282 if (!fDigitsTree->GetEvent(index)) continue;
283 Int_t sec,row;
284 fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
285 fDigits->First();
286 do {
287 Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
288 Short_t dig = fDigits->GetDigit(it,ip);
289 Int_t idx0=fDigits->GetTrackID(it,ip,0);
290 Int_t idx1=fDigits->GetTrackID(it,ip,1);
291 Int_t idx2=fDigits->GetTrackID(it,ip,2);
e170e0f9 292
e1a32fa5 293 if (idx0>=0 && dig>=zero) count[idx0]+=1;
294 if (idx1>=0 && dig>=zero) count[idx1]+=1;
295 if (idx2>=0 && dig>=zero) count[idx2]+=1;
296 } while (fDigits->Next());
297
298 for (Int_t j=0; j<np; j++)
299 {
300 if (count[j]>1) //at least two digits at this padrow
301 good[j]++;
302
303 count[j]=0;
304 }
108615fc 305 }
108615fc 306 }
108615fc 307 delete[] count;
308 }
e1a32fa5 309
108615fc 310 else
311 {
312 LOG(AliL3Log::kError,"AliL3Evaluation::FillEffHistos","TPC version")
313 <<"No valid TPC version found"<<ENDLOG;
e1a32fa5 314 return;
108615fc 315 }
316
e1a32fa5 317
318 //The following code has been taken from offlinemacro->AliTPCComparison.C
319
320 TTree *TH=gAlice->TreeH();
321 Int_t npart=(Int_t)TH->GetEntries();
322 Int_t nt=0;
323 Int_t max = 15000;
324 while (npart--) {
325 AliTPChit *hit0=0;
326
327 TPC->ResetHits();
328 TH->GetEvent(npart);
329 AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
330 while (hit){
331 if (hit->fQ==0.) break;
332 hit = (AliTPChit*) TPC->NextHit();
333 }
334 if (hit) {
335 hit0 = new AliTPChit(*hit); //Make copy of hit
336 hit = hit0;
337 }
338 else continue;
339 AliTPChit *hit1=(AliTPChit*)TPC->NextHit();
340 if (hit1==0) continue;
341 if (hit1->fQ != 0.) continue;
342 Int_t i=hit->Track();
343 TParticle *p = (TParticle*)gAlice->Particle(i);
e170e0f9 344
345 printf("Checking particle %d with code %d\n",i,p->GetPdgCode());
e1a32fa5 346 if (p->GetFirstMother()>=0) continue; //secondary particle
347 if (good[i] < good_number) continue;
e170e0f9 348 if (p->Pt()<fMinGoodPt) continue;
e1a32fa5 349 if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
e170e0f9 350 printf("Checking particle %d, nHits %d\n",i,good[i]);
e1a32fa5 351 fGoodGen++;
352 fGoodTracks[nt].label=i;
353 fGoodTracks[nt].code=p->GetPdgCode();
354 //**** px py pz - in global coordinate system, x y z - in local !
355 fGoodTracks[nt].px=hit->X(); fGoodTracks[nt].py=hit->Y(); fGoodTracks[nt].pz=hit->Z();
356
357 nt++;
358 if (hit0) delete hit0;
359 if (nt==max) {cerr<<"Too many good tracks !n"; break;}
108615fc 360 }
e170e0f9 361
108615fc 362 delete [] good;
e1a32fa5 363
108615fc 364}
365
366void AliL3Evaluate::EvaluatePatch(Int_t slice,Int_t patch,Int_t min_points,Int_t good_number)
367{
368 //Make efficiency plots for tracking on patch level (before any merging).
369
ae97a0b9 370 Int_t row[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,175}};
e1a32fa5 371 Int_t sl[2] ={slice,slice};
372 DefineGoodTracks(sl,row[patch],good_number);
108615fc 373 SetMinPoints(min_points);
374 AssignIDs();
bc2f4f0e 375 CreateHistos();
e1a32fa5 376 FillEffHistos();
bc2f4f0e 377 CalcEffHistos();
108615fc 378}
379
380void AliL3Evaluate::EvaluateSlice(Int_t slice,Int_t min_points,Int_t good_number)
381{
382 //Make efficiency plots for tracking on a slice (after merging).
383 //min_points = minimum points on track to be considered for evaluation
384 //good_number = minimum hits (padrows) produced by simulated track for consideration.
385
ae97a0b9 386 Int_t row[2] = {0,175};
e1a32fa5 387 Int_t sl[2] ={slice,slice};
388 DefineGoodTracks(sl,row,good_number);
23908b9b 389
108615fc 390 SetMinPoints(min_points);
23908b9b 391
108615fc 392 AssignIDs();
bc2f4f0e 393 CreateHistos();
e1a32fa5 394 FillEffHistos();
bc2f4f0e 395 CalcEffHistos();
108615fc 396}
397
bc2f4f0e 398void AliL3Evaluate::EvaluateGlobal(Int_t min_points,Int_t good_number)
108615fc 399{
400 //Make efficiency plots for tracking on several slices.
bc2f4f0e 401
ae97a0b9 402 Int_t row[2] = {0,175};
bc2f4f0e 403 SetMinPoints(min_points);
e1a32fa5 404 Int_t slice[2] = {fMinSlice,fMaxSlice};
405 DefineGoodTracks(slice,row,good_number);
bc2f4f0e 406 AssignIDs();
23908b9b 407 CreateHistos(20,0,5);
e1a32fa5 408 FillEffHistos();
bc2f4f0e 409 CalcEffHistos();
108615fc 410}
411
108615fc 412void AliL3Evaluate::AssignIDs()
413{
414 //Assign MC id to the tracks.
415
108615fc 416 fTracks->QSort();
417 LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
418 <<"Assigning MC id to the found tracks...."<<ENDLOG;
419 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
420 {
421 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
c0217812 422 if(!track) continue;
108615fc 423 if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
23908b9b 424
23908b9b 425 fGoodFound++;
c0217812 426 Int_t tID = GetMCTrackLabel(track);
108615fc 427 track->SetMCid(tID);
e170e0f9 428 printf("track %i id %d nHits %d\n",i,tID,track->GetNumberOfPoints());
108615fc 429 }
108615fc 430}
431
432
433struct S {Int_t lab; Int_t max;};
c0217812 434Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
108615fc 435 //Returns the MCtrackID of the belonging clusters.
436 //If MCLabel < 0, means that track is fake.
437 //Fake track means that more than 10 percent of clusters are assigned incorrectly.
c0217812 438
108615fc 439 Int_t num_of_clusters = track->GetNumberOfPoints();
440 S *s=new S[num_of_clusters];
441 Int_t i;
442 for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
443 UInt_t *hitnum = track->GetHitNumbers();
444 UInt_t id;
108615fc 445
23908b9b 446 Int_t **trackID = GetClusterIDs(track);
108615fc 447
448 Int_t lab=123456789;
449 for (i=0; i<num_of_clusters; i++)
450 {
451 //Tricks to get the clusters belonging to this track:
452 id = hitnum[i];
453 Int_t slice = (id>>25) & 0x7f;
454 Int_t patch = (id>>22) & 0x7;
455 UInt_t pos = id&0x3fffff;
456
457 AliL3SpacePointData *points = fClusters[slice][patch];
c0217812 458 if(!points) continue;
108615fc 459 if(pos>=fNcl[slice][patch])
460 {
461 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
462 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
463 continue;
464 }
465
466 //Get the label of the cluster:
467 //printf("label %d %d %d\n",trackID[i][0],trackID[i][1],trackID[i][2]);
468 lab=abs(trackID[i][0]);
469 Int_t j;
470 for (j=0; j<num_of_clusters; j++)
471 if (s[j].lab==lab || s[j].max==0) break;
472 s[j].lab=lab;
473 s[j].max++;
474 }
475
476 Int_t max=0;
477 for (i=0; i<num_of_clusters; i++)
478 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
479
480 delete[] s;
481
482 for (i=0; i<num_of_clusters; i++)
483 {
484 id = hitnum[i];
485 Int_t slice = (id>>25) & 0x7f;
486 Int_t patch = (id>>22) & 0x7;
487 UInt_t pos = id&0x3fffff;
488
489 AliL3SpacePointData *points = fClusters[slice][patch];
c0217812 490 if(!points) continue;
108615fc 491 if(pos>=fNcl[slice][patch])
492 {
493 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
494 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
495 continue;
496 }
497
498 if (abs(trackID[i][1]) == lab ||
499 abs(trackID[i][2]) == lab ) max++;
500 }
501
502 //check if more than 10% of the clusters are incorrectly assigned (fake track):
503 if (1.-Float_t(max)/num_of_clusters > 0.10)
504 {
505 return -lab;
506 }
507
e170e0f9 508 for(Int_t j=0; j<track->GetNumberOfPoints(); j++)
509 delete [] trackID[j];
108615fc 510 delete [] trackID;
511 return lab;
512}
513
514
c0217812 515Int_t **AliL3Evaluate::GetClusterIDs(AliL3Track *track)
108615fc 516{
517 //Return the MC information of all clusters belonging to track.
c0217812 518
108615fc 519 Int_t num_of_clusters = track->GetNumberOfPoints();
520 Int_t **trackID = new Int_t*[num_of_clusters];
521
522 UInt_t *hitnum = track->GetHitNumbers();
523 UInt_t id;
108615fc 524
525 Float_t xyz[3];
c0217812 526 Int_t padrow;
108615fc 527 for(Int_t i=0; i<num_of_clusters; i++)
528 {
529 id = hitnum[i];
530 Int_t slice = (id>>25) & 0x7f;
531 Int_t patch = (id>>22) & 0x7;
532 UInt_t pos = id&0x3fffff;
533
534 AliL3SpacePointData *points = fClusters[slice][patch];
535
536 if(!points)
537 {
538 LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
539 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
540 continue;
541 }
542 if(pos>=fNcl[slice][patch])
543 {
544 LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Clusterarray")
545 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
546 continue;
547 }
548
549 xyz[0] = points[pos].fX;
550 xyz[1] = points[pos].fY;
551 xyz[2] = points[pos].fZ;
552 //sector = points[pos].fSector;
553 padrow = points[pos].fPadRow;
554 Int_t se,ro;
555 fTransform->Slice2Sector(slice,padrow,se,ro);
556 fTransform->Global2Raw(xyz,se,ro);
557
558 if(fIsSlow)
559 {
560 Int_t p = fRowid[slice][padrow];
108615fc 561 if(!fDigitsTree->GetEvent(p))
562 LOG(AliL3Log::kError,"AliL3Evaluate::GetClusterIDs","Digits Tree")
563 <<"Error reading digits tree"<<ENDLOG;
564
565 trackID[i] = new Int_t[3];
e170e0f9 566 trackID[i][0] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0);
567 trackID[i][1] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),1);
568 trackID[i][2] = fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),2);
569 //if(trackID[i][0]==6 || trackID[i][0]==32)
570 //printf("trackID %d, padrow %d pad %d time %d\n",fDigits->GetTrackID((Int_t)rint(xyz[2]),(Int_t)rint(xyz[1]),0),padrow,(int)rint(xyz[1]),(int)rint(xyz[2]));
571 /*if(trackID[i][0]<0)
572 {
573
574 printf("trackID %d, padrow %d pad %d time %d\n",trackID[i][0],padrow,(int)rint(xyz[1]),(int)rint(xyz[2]));
575 printf("on the side %d %d %d %d\n",fDigits->GetTrackID(((int)rint(xyz[2])-1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])+1),((int)rint(xyz[1])),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])-1),0),fDigits->GetTrackID(((int)rint(xyz[2])),((int)rint(xyz[1])+1),0));
576 }*/
108615fc 577 }
578 else
579 {
580 Int_t tmp_pid=0;
c0217812 581 for(Int_t ii=0; ii<fNFastPoints; ii++)
108615fc 582 {
c0217812 583 tmp_pid = fMcId[ii];
584 if(fMcIndex[ii] == id) break;
108615fc 585 }
586 trackID[i] = new Int_t[3];
587 trackID[i][0] = tmp_pid;
588 trackID[i][1] = -1;
589 trackID[i][2] = -1;
590 }
591 }
e170e0f9 592
108615fc 593 return trackID;
594}
595
c0217812 596void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
108615fc 597{
598 //Get the MC id of space points in case of using the fast simulator.
c0217812 599 char fname[256];
600 sprintf(fname,"%s/point_mc.dat",path);
601 FILE *infile = fopen(fname,"r");
602 if(!infile) return;
108615fc 603 Int_t hitid,hitmc,i;
604
605 for(i=0; ; i++)
606 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
108615fc 607 rewind(infile);
c0217812 608 fNFastPoints = i;
609 fMcId = new Int_t[fNFastPoints];
610 fMcIndex = new UInt_t[fNFastPoints];
108615fc 611
c0217812 612 for(i=0; i<fNFastPoints; i++)
108615fc 613 {
614 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
c0217812 615 fMcId[i] = hitmc;
616 fMcIndex[i] = hitid;
108615fc 617 }
618 fclose(infile);
c0217812 619}
620
621void AliL3Evaluate::CreateHistos(Int_t nbin,Int_t xlow,Int_t xup)
622{
623 //Create the histograms
624
1f79afc0 625 fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
626
c0217812 627 fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
628 fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
629 fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
630 fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
631 fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
632 fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
633
634 fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
635 fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
636 fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
637 fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
638 fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
639}
108615fc 640
e1a32fa5 641void AliL3Evaluate::FillEffHistos()
c0217812 642{
643 //Fill the efficiency histograms.
108615fc 644
e1a32fa5 645 for(Int_t i=0; i<fGoodGen; i++)
c0217812 646 {
e1a32fa5 647 Double_t ptg=fGoodTracks[i].pt,pzg=fGoodTracks[i].pz;
c0217812 648 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
649 fNGoodTracksPt->Fill(ptg);
650 fNGoodTracksEta->Fill(dipangle);
651 Int_t found = 0;
652 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
653 {
654 AliL3Track *track = fTracks->GetCheckedTrack(k);
655 if(!track) continue;
656 Int_t nHits = track->GetNumberOfPoints();
657 if(nHits < fMinPointsOnTrack) break;
658
659 Int_t tracklabel;
660 tracklabel = track->GetMCid();
661
e1a32fa5 662 if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
c0217812 663 found=1;
e1a32fa5 664 if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
c0217812 665 else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
666 Float_t pt=track->GetPt();
667 fPtRes->Fill((pt-ptg)/ptg*100.);
1f79afc0 668 fNtuppel->Fill(ptg,pt,nHits);
c0217812 669 break;
670
671 }
672 }
673}
674
675void AliL3Evaluate::CalcEffHistos(){
676
677 Stat_t ngood=fNGoodTracksPt->GetEntries();
678 Stat_t nfound=fNFoundTracksPt->GetEntries();
679 Stat_t nfake=fNFakeTracksPt->GetEntries();
680
681 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
682 <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
683 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
684 <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
685 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
686 <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
687 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
688 <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
689
690 fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
691 fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
692 fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
693 fTrackEffPt->SetMaximum(1.4);
694 fTrackEffPt->SetXTitle("P_{T} [GeV]");
695 fTrackEffPt->SetLineWidth(2);
696 fFakeTrackEffPt->SetFillStyle(3013);
697 fTrackEffPt->SetLineColor(4);
698 fFakeTrackEffPt->SetFillColor(2);
699
700 fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
701 fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
702 fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
703 fTrackEffEta->SetMaximum(1.4);
704 fTrackEffEta->SetXTitle("#lambda [degrees]");
705 fTrackEffEta->SetLineWidth(2);
706 fFakeTrackEffEta->SetFillStyle(3013);
707 fTrackEffEta->SetLineColor(4);
708 fFakeTrackEffEta->SetFillColor(2);
709
108615fc 710}
711
712void AliL3Evaluate::Write2File(Char_t *outputfile)
713{
714 //Write histograms to file:
715
716 TFile *of = new TFile(outputfile,"RECREATE");
717 if(!of->IsOpen())
718 {
719 LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
720 <<"Problems opening rootfile"<<ENDLOG;
721 return;
722 }
723
724 of->cd();
1f79afc0 725 fNtuppel->Write();
108615fc 726 fPtRes->Write();
727 fNGoodTracksPt->Write();
728 fNFoundTracksPt->Write();
729 fNFakeTracksPt->Write();
730 fTrackEffPt->Write();
731 fFakeTrackEffPt->Write();
732 fNGoodTracksEta->Write();
733 fNFoundTracksEta->Write();
734 fNFakeTracksEta->Write();
735 fTrackEffEta->Write();
736 fFakeTrackEffEta->Write();
737
738 of->Close();
739 delete of;
740
741}
742
1f79afc0 743TNtuple *AliL3Evaluate::CalculateResiduals()
744{
745
746 TNtuple *ntuppel=new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
78127e35 747
1f79afc0 748 for(int f=fMinSlice; f<=fMaxSlice; f++)
749 {
750 AliL3FileHandler *tfile = new AliL3FileHandler();
751 char fname[256];
752 sprintf(fname,"tracks_tr_%d_0.raw",f);
753 if(!tfile->SetBinaryInput(fname)){
754 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
755 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
756 return 0;
757 }
758 fTracks = new AliL3TrackArray();
759 tfile->Binary2TrackArray(fTracks);
760 tfile->CloseBinaryInput();
761 delete tfile;
78127e35 762 printf("Looking in slice %d\n",f);
1f79afc0 763 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
764 {
765
766 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
767 if(!track) continue;
768
769 track->CalculateHelix();
770 UInt_t *hitnum = track->GetHitNumbers();
771 UInt_t id;
772
773 Float_t xyz[3];
774 Int_t padrow;
78127e35 775 for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
1f79afc0 776 {
777 id = hitnum[j];
778 Int_t slice = (id>>25) & 0x7f;
779 Int_t patch = (id>>22) & 0x7;
780 UInt_t pos = id&0x3fffff;
781
78127e35 782 //if(slice!=1) continue;
783
1f79afc0 784 AliL3SpacePointData *points = fClusters[slice][patch];
785
786 if(!points)
787 {
788 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
789 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
790 continue;
791 }
792 if(pos>=fNcl[slice][patch])
793 {
794 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
795 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
796 continue;
797 }
798
799 xyz[0] = points[pos].fX;
800 xyz[1] = points[pos].fY;
801 xyz[2] = points[pos].fZ;
802 padrow = points[pos].fPadRow;
78127e35 803 //fTransform->Global2Local(xyz,slice);
1f79afc0 804
805 Float_t xyz_cross[3];
806 track->GetCrossingPoint(padrow,xyz_cross);
e1a32fa5 807 Double_t beta = track->GetCrossingAngle(padrow);
1f79afc0 808
809 Double_t yres = xyz_cross[1] - xyz[1];
810 Double_t zres = xyz_cross[2] - xyz[2];
811
812 Double_t dipangle = atan(track->GetTgl());
813 ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
814
815 }
816 }
817 if(fTracks)
818 delete fTracks;
819 }
820 return ntuppel;
821}
822
e1a32fa5 823TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *rootfile)
78127e35 824{
e1a32fa5 825 //Compare points to the exact crossing points of track and padrows.
1b76fd18 826 //The input file to this function, contains the exact clusters calculated
827 //in AliTPC::Hits2ExactClusters.
e1a32fa5 828
78127e35 829
e1a32fa5 830 TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt");
78127e35 831
e1a32fa5 832 TFile *exfile = new TFile(rootfile);
833
834 AliTPCParam *param = (AliTPCParam*)exfile->Get("75x40_100x60");
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",param->GetTitle());
842 Bool_t clusterok = arr->ConnectTree(treeName);
843 if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;}
844
845 for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
78127e35 846 {
e1a32fa5 847 //Get the exact clusters for this row:
848 Int_t cursec,currow;
849 AliSegmentID *s = arr->LoadEntry(i);
850 param->AdjustSectorRow(s->GetID(),cursec,currow);
78127e35 851
e1a32fa5 852 AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
853 TClonesArray *clusters = ro->GetArray();
854 int num_of_offline=clusters->GetEntriesFast();
78127e35 855
e1a32fa5 856 //Get the found clusters:
857 Int_t slice,padrow;
858 fTransform->Sector2Slice(slice,padrow,cursec,currow);
859 if(slice<fMinSlice || slice>fMaxSlice) continue;
860 AliL3SpacePointData *points = fClusters[slice][0];
1b76fd18 861
862 Int_t index = fRowid[slice][padrow];
863 if(!fDigitsTree->GetEvent(index))
864 printf("AliL3Evaluate::EvaluatePoints : ERROR IN DIGITSTREE\n");
865 printf("Checking padrow %d\n",padrow);
866
e1a32fa5 867 for(UInt_t c=0; c<fNcl[slice][0]; c++)
78127e35 868 {
e1a32fa5 869 if(points[c].fPadRow!=padrow) continue;
e1a32fa5 870 for(Int_t m=0; m<num_of_offline; m++)
78127e35 871 {
e1a32fa5 872 AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
873 Int_t mcId = cluster->fTracks[0];
874 TParticle *part = gAlice->Particle(mcId);
875 Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
876 fTransform->Global2Raw(xyz_cl,cursec,currow);
e170e0f9 877 if(fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),0)!=mcId &&
878 fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),1)!=mcId &&
879 fDigits->GetTrackID((Int_t)rint(xyz_cl[2]),(Int_t)rint(xyz_cl[1]),2)!=mcId)
e1a32fa5 880 continue;
881
882 Float_t resy = xyz_cl[1]-cluster->fY;
883 Float_t resz = xyz_cl[2]-cluster->fX;
884 ntuppel->Fill(slice,padrow,resy,resz,cluster->fX,part->Pt());
78127e35 885 }
e1a32fa5 886 }
78127e35 887 }
e1a32fa5 888
78127e35 889 return ntuppel;
890}
e1a32fa5 891/*
78127e35 892TNtuple *AliL3Evaluate::EvaluateGEANT()
893{
894 TNtuple *ntuppel = new TNtuple("ntuppel","residuals","resy:ptgen:padrow:zHit:slice");
895
ae97a0b9 896 Int_t good_number=175;
897 Int_t row[2] = {0,175};
78127e35 898
899 Float_t xyz_cross[3];
900 Float_t xyz_cl[3];
901
902 for(Int_t slice=fMinSlice; slice<=fMaxSlice; slice++)
903 {
904 Int_t *particle_id = new Int_t[fParticles->GetEntriesFast()];
905 TObjArray *good_particles = DefineGoodTracks(slice,row,good_number,particle_id);
906
907 if(fMCclusterfile)
908 {
909 fMCclusterfile->Close();
910 delete fMCclusterfile;
911 fMCclusterfile = new TFile("/nfs/david/subatom/alice/data/V3.04/fast/clusters/hg_8k_v0_s1-3_e0_cl.root");
912 fMCclusterfile->cd();
913 }
914 AliTPCClustersArray arr;
915 arr.Setup(fParam);
916 arr.SetClusterType("AliTPCcluster");
917 Bool_t clusterok = arr.ConnectTree("Segment Tree");
918 if(!clusterok)
919 LOG(AliL3Log::kError,"AliL3Evaluate::DefineGoodTracks","Cluster Array")
920 <<"Error loading clusters from rootfile"<<ENDLOG;
921
922 for(Int_t ci=0; ci<arr.GetTree()->GetEntries(); ci++)
923 {
924 Int_t sec,row,sl,lr;
925 AliSegmentID *s = arr.LoadEntry(ci);
926 fParam->AdjustSectorRow(s->GetID(),sec,row);
927 fTransform->Sector2Slice(sl,lr,sec,row);
928 if(sl != slice) {arr.ClearRow(sec,row); continue;}
929 //if(lr!=pr) {carray.ClearRow(sec,row); continue;}
930 AliTPCClustersRow *cRow = arr.GetRow(sec,row);
931 printf("evaluating slice %d row %d\n",slice,lr);
932 for(Int_t j=0; j<cRow->GetArray()->GetEntriesFast(); j++)
933 {
934
935 AliTPCcluster *cluster=(AliTPCcluster*)(*cRow)[j];
936 xyz_cl[1] = cluster->GetY();
937 xyz_cl[2] = cluster->GetZ();
938 Int_t lab=cluster->GetLabel(0);
939 if(lab<0) continue;
940
941 for(Int_t pa=0; pa<good_particles->GetEntriesFast(); pa++)
942 {
943 if(particle_id[pa]!=lab) continue;
944 TParticle *part = (TParticle*)good_particles->UncheckedAt(pa);
945 GetParticleCrossingPoint(part,slice,lr,xyz_cross);
946 Double_t yres = xyz_cl[1] - xyz_cross[1];
947 ntuppel->Fill(yres,part->Pt(),lr,xyz_cl[2],slice);
948 }
949 }
950 if(arr.GetRow(sec,row))
951 arr.ClearRow(sec,row);
952 }
953 delete [] particle_id;
954 delete good_particles;
955 }
956 return ntuppel;
957}
e1a32fa5 958*/
78127e35 959Bool_t AliL3Evaluate::GetParticleCrossingPoint(TParticle *part,Int_t slice,Int_t padrow,Float_t *xyz)
960{
961 //Calcluate the crossing point between a generated particle and given padrow.
962
963
964 Double_t kappa = 0.2*0.0029980/part->Pt();
965
966 Double_t radius = 1/fabs(kappa);
967 if(part->GetPdgCode() > 0) kappa = -kappa;
968
23908b9b 969 Int_t charg;
970 if(kappa>0)
971 charg=-1;
972 else
973 charg=1;
974
78127e35 975 Float_t angl[1] = {part->Phi()};
976
977 fTransform->Global2LocalAngle(angl,slice);
978
979 Double_t charge = -1.*kappa;
980 Double_t trackPhi0 = angl[0] + charge*0.5*Pi/fabs(charge);
981
982 Double_t x0=0;
983 Double_t y0=0;
984 Double_t xc = x0 - radius * cos(trackPhi0);
985 Double_t yc = y0 - radius * sin(trackPhi0);
986
987 //printf("radius %f xc %f yc %f\n",radius,xc,yc);
988
989 Double_t xHit = fTransform->Row2X(padrow);
990 xyz[0] = xHit;
991 Double_t aa = (xHit - xc)*(xHit - xc);
992 Double_t r2 = radius*radius;
993 if(aa > r2)
994 return false;
995
996 Double_t aa2 = sqrt(r2 - aa);
997 Double_t y1 = yc + aa2;
998 Double_t y2 = yc - aa2;
999 xyz[1] = y1;
1000 if(fabs(y2) < fabs(y1)) xyz[1] = y2;
1001
23908b9b 1002 Double_t tgl = part->Pz()/part->Pt();
1003
1004 Double_t yHit = xyz[1];
1005 Double_t angle1 = atan2((yHit - yc),(xHit - xc));
1006 if(angle1 < 0) angle1 += 2.*Pi;
1007 Double_t angle2 = atan2((0 - yc),(0 - xc));
1008 if(angle2 < 0) angle2 += 2.*Pi;
1009 Double_t diff_angle = angle1 - angle2;
1010 diff_angle = fmod(diff_angle,2*Pi);
1011 if((charg*diff_angle) > 0) diff_angle = diff_angle - charg*2.*Pi;
1012 Double_t s_tot = fabs(diff_angle)*radius;
1013 Double_t zHit = 0 + s_tot*tgl;
1014 xyz[2] = zHit;
1015
78127e35 1016 return true;
1017}
1018