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