Updated function EvaluatePoints.
[u/mrichter/AliRoot.git] / HLT / src / AliL3Evaluate.cxx
CommitLineData
e8541476 1//$Id$
2
b661165c 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
108615fc 5
118c26c3 6#include "AliL3StandardIncludes.h"
108615fc 7#include <TFile.h>
8#include <TH1.h>
9#include <TParticle.h>
10#include <TTree.h>
11#include <TClonesArray.h>
a8a3f5c2 12#include <TNtupleD.h>
108615fc 13
4499ed26 14#include <AliRun.h>
15#include <AliSimDigits.h>
16#include <AliTPC.h>
17#include <AliTPCcluster.h>
18#include <AliTPCClustersArray.h>
19#include <AliTPCClustersRow.h>
20#include <AliTPCParam.h>
21#include <AliComplexCluster.h>
108615fc 22
118c26c3 23#if GCCVERSION == 3
24#include <fstream>
25#include <iosfwd>
26#else
27#include <fstream.h>
28#endif
29
30#include "AliL3Logging.h"
108615fc 31#include "AliL3Transform.h"
32#include "AliL3SpacePointData.h"
33#include "AliL3Track.h"
34#include "AliL3FileHandler.h"
35#include "AliL3TrackArray.h"
36#include "AliL3Evaluate.h"
118c26c3 37
38#if GCCVERSION == 3
39using namespace std;
40#endif
108615fc 41
b661165c 42//_____________________________________________________________
43// AliL3Evaluate
44//
45// Evaluation class for tracking. Plots efficiencies etc..
84122e15 46//
47
108615fc 48
49ClassImp(AliL3Evaluate)
50
51AliL3Evaluate::AliL3Evaluate()
52{
e0f4d6b2 53 fTracks = 0;
c0217812 54 fNFastPoints = 0;
55 fMcIndex = 0;
56 fMcId = 0;
1f79afc0 57 fMinSlice=0;
58 fMaxSlice=0;
e0f4d6b2 59 fGoodTracks=0;
108615fc 60}
61
e0f4d6b2 62AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Int_t *slice)
108615fc 63{
108615fc 64
e0f4d6b2 65 if(slice)
66 {
67 fMinSlice=slice[0];
68 fMaxSlice=slice[1];
69 }
70 else
71 {
72 fMinSlice=0;
73 fMaxSlice=35;
74 }
75 fMaxFalseClusters = 0.1;
23908b9b 76 fGoodFound = 0;
77 fGoodGen = 0;
e0f4d6b2 78 fMinPointsOnTrack = min_clusters;
79 fMinHitsFromParticle = minhits;
80 fMinGoodPt = minpt;
81 Char_t fname[1024];
40896357 82 AliL3FileHandler *clusterfile[36][6];
108615fc 83 for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
84 {
e0f4d6b2 85 for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
108615fc 86 {
1f79afc0 87 fClusters[s][p] = 0;
108615fc 88 clusterfile[s][p] = new AliL3FileHandler();
e0f4d6b2 89 sprintf(fname,"%s/points_%d_%d.raw",datapath,s,p);
108615fc 90 if(!clusterfile[s][p]->SetBinaryInput(fname))
91 {
92 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
93 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
bc2f4f0e 94 delete clusterfile[s][p];
95 clusterfile[s][p] = 0;
96 continue;
108615fc 97 }
98 fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
99 clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
100 clusterfile[s][p]->CloseBinaryInput();
101 }
102 }
78127e35 103
e0f4d6b2 104 sprintf(fname,"%s/tracks.raw",datapath);
108615fc 105 AliL3FileHandler *tfile = new AliL3FileHandler();
e0f4d6b2 106 if(!tfile->SetBinaryInput(fname)){
bc2f4f0e 107 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
e0f4d6b2 108 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
bc2f4f0e 109 return;
110 }
108615fc 111 fTracks = new AliL3TrackArray();
112 tfile->Binary2TrackArray(fTracks);
113 tfile->CloseBinaryInput();
114 delete tfile;
bc2f4f0e 115}
116
108615fc 117
e0f4d6b2 118AliL3Evaluate::~AliL3Evaluate()
108615fc 119{
e0f4d6b2 120 if(fGoodTracks) delete fGoodTracks;
121 if(fDigitsTree) fDigitsTree->Delete();
122 if(fTracks) delete fTracks;
123 if(fPtRes) delete fPtRes;
124 if(fNGoodTracksPt) delete fNGoodTracksPt;
125 if(fNFoundTracksPt) delete fNFoundTracksPt;
126 if(fNFakeTracksPt) delete fNFakeTracksPt;
127 if(fTrackEffPt) delete fTrackEffPt;
128 if(fFakeTrackEffPt) delete fFakeTrackEffPt;
129 if(fNGoodTracksEta) delete fNGoodTracksEta;
130 if(fNFoundTracksEta) delete fNFoundTracksEta;
131 if(fNFakeTracksEta) delete fNFakeTracksEta;
132 if(fTrackEffEta) delete fTrackEffEta;
133 if(fFakeTrackEffEta) delete fFakeTrackEffEta;
134 if(fMcIndex) delete [] fMcIndex;
135 if(fMcId) delete [] fMcId;
136 if(fNtuppel) delete fNtuppel;
108615fc 137}
138
8af8a456 139
108615fc 140void AliL3Evaluate::AssignIDs()
141{
142 //Assign MC id to the tracks.
e0f4d6b2 143#ifndef do_mc
144 cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
145 return;
146#endif
108615fc 147
108615fc 148 fTracks->QSort();
149 LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
150 <<"Assigning MC id to the found tracks...."<<ENDLOG;
151 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
152 {
153 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
c0217812 154 if(!track) continue;
108615fc 155 if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
23908b9b 156
23908b9b 157 fGoodFound++;
c0217812 158 Int_t tID = GetMCTrackLabel(track);
108615fc 159 track->SetMCid(tID);
108615fc 160 }
e0f4d6b2 161 cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
108615fc 162}
163
164
165struct S {Int_t lab; Int_t max;};
c0217812 166Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
108615fc 167 //Returns the MCtrackID of the belonging clusters.
168 //If MCLabel < 0, means that track is fake.
169 //Fake track means that more than 10 percent of clusters are assigned incorrectly.
494fad94 170
171#ifdef do_mc
108615fc 172 Int_t num_of_clusters = track->GetNumberOfPoints();
173 S *s=new S[num_of_clusters];
174 Int_t i;
175 for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
176 UInt_t *hitnum = track->GetHitNumbers();
177 UInt_t id;
108615fc 178
108615fc 179 Int_t lab=123456789;
180 for (i=0; i<num_of_clusters; i++)
181 {
182 //Tricks to get the clusters belonging to this track:
183 id = hitnum[i];
184 Int_t slice = (id>>25) & 0x7f;
185 Int_t patch = (id>>22) & 0x7;
186 UInt_t pos = id&0x3fffff;
187
188 AliL3SpacePointData *points = fClusters[slice][patch];
c0217812 189 if(!points) continue;
108615fc 190 if(pos>=fNcl[slice][patch])
191 {
192 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
193 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
194 continue;
195 }
196
197 //Get the label of the cluster:
355debe1 198
199 lab=abs(points[pos].fTrackID[0]);
e0f4d6b2 200 if(lab < 0)
201 cout<<"Track had negative id : "<<lab<<" padrow "<<(Int_t)points[pos].fPadRow<<" nhits "<<num_of_clusters<<" pt "<<track->GetPt()<<endl;
202
108615fc 203 Int_t j;
204 for (j=0; j<num_of_clusters; j++)
205 if (s[j].lab==lab || s[j].max==0) break;
206 s[j].lab=lab;
207 s[j].max++;
208 }
209
210 Int_t max=0;
211 for (i=0; i<num_of_clusters; i++)
212 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
213
214 delete[] s;
215
216 for (i=0; i<num_of_clusters; i++)
217 {
218 id = hitnum[i];
219 Int_t slice = (id>>25) & 0x7f;
220 Int_t patch = (id>>22) & 0x7;
221 UInt_t pos = id&0x3fffff;
222
223 AliL3SpacePointData *points = fClusters[slice][patch];
c0217812 224 if(!points) continue;
108615fc 225 if(pos>=fNcl[slice][patch])
226 {
227 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
228 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
229 continue;
230 }
231
355debe1 232 if (abs(points[pos].fTrackID[1]) == lab ||
233 abs(points[pos].fTrackID[2]) == lab ) max++;
108615fc 234 }
235
e0f4d6b2 236 if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters)
108615fc 237 {
238 return -lab;
239 }
240
108615fc 241 return lab;
494fad94 242#else //If we are running with mc_ids or not
243 return 0;
244#endif
245
108615fc 246}
247
c0217812 248void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
108615fc 249{
250 //Get the MC id of space points in case of using the fast simulator.
c0217812 251 char fname[256];
252 sprintf(fname,"%s/point_mc.dat",path);
253 FILE *infile = fopen(fname,"r");
254 if(!infile) return;
108615fc 255 Int_t hitid,hitmc,i;
256
257 for(i=0; ; i++)
258 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
108615fc 259 rewind(infile);
c0217812 260 fNFastPoints = i;
261 fMcId = new Int_t[fNFastPoints];
262 fMcIndex = new UInt_t[fNFastPoints];
108615fc 263
c0217812 264 for(i=0; i<fNFastPoints; i++)
108615fc 265 {
266 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
c0217812 267 fMcId[i] = hitmc;
268 fMcIndex[i] = hitid;
108615fc 269 }
270 fclose(infile);
c0217812 271}
272
e0f4d6b2 273void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
c0217812 274{
275 //Create the histograms
276
40896357 277 LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
278 <<"Creating histograms..."<<ENDLOG;
279
1f79afc0 280 fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
a8a3f5c2 281 fNtuppel->SetDirectory(0);
c0217812 282 fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
283 fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
284 fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
285 fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
286 fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
287 fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
288
289 fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
290 fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
291 fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
292 fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
293 fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
e0f4d6b2 294
295}
296
a8a3f5c2 297void AliL3Evaluate::GetGoodParticles(Char_t *path,Bool_t sector)
e0f4d6b2 298{
299 //Read the good particles from file. This file should already have been
300 //generated by macro AliTPCComparison.C.
301
302 Char_t filename[1024];
a8a3f5c2 303 if(!sector)
304 sprintf(filename,"%s/good_tracks_tpc",path);
305 else
306 sprintf(filename,"%s/good_tracks_tpc_sector",path);//Sectorwise comparison.
e0f4d6b2 307 ifstream in(filename);
308 if(!in)
309 {
310 cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
311 return;
312 }
313 Int_t MaxTracks=20000;
314 fGoodTracks = new GoodTrack[MaxTracks];
315
316 while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
317 fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
a8a3f5c2 318 fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector)
e0f4d6b2 319
320 {
321 fGoodGen++;
322 if (fGoodGen==MaxTracks)
323 {
324 cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
325 break;
326 }
327 }
c0217812 328}
eeddc64d 329
a8a3f5c2 330
e1a32fa5 331void AliL3Evaluate::FillEffHistos()
c0217812 332{
e0f4d6b2 333 if(!fGoodTracks)
334 {
335 cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
336 return;
337 }
a8a3f5c2 338 cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
e1a32fa5 339 for(Int_t i=0; i<fGoodGen; i++)
c0217812 340 {
e0f4d6b2 341 //cout<<"Checking particle "<<i<<endl;
342 if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
343 Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
a8a3f5c2 344 if(ptg < fMinGoodPt) continue;
e0f4d6b2 345 Float_t pzg=fGoodTracks[i].pz;
c0217812 346 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
e0f4d6b2 347
c0217812 348 fNGoodTracksPt->Fill(ptg);
349 fNGoodTracksEta->Fill(dipangle);
350 Int_t found = 0;
e0f4d6b2 351
c0217812 352 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
353 {
354 AliL3Track *track = fTracks->GetCheckedTrack(k);
355 if(!track) continue;
356 Int_t nHits = track->GetNumberOfPoints();
357 if(nHits < fMinPointsOnTrack) break;
c0217812 358 Int_t tracklabel;
359 tracklabel = track->GetMCid();
e0f4d6b2 360
e1a32fa5 361 if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
c0217812 362 found=1;
e1a32fa5 363 if(tracklabel == fGoodTracks[i].label) {fNFoundTracksPt->Fill(ptg); fNFoundTracksEta->Fill(dipangle);}
c0217812 364 else {fNFakeTracksPt->Fill(ptg); fNFakeTracksEta->Fill(dipangle);}
365 Float_t pt=track->GetPt();
366 fPtRes->Fill((pt-ptg)/ptg*100.);
a8a3f5c2 367 fNtuppel->Fill(ptg,pt,nHits);
c0217812 368 break;
369
370 }
a8a3f5c2 371 if(!found)
372 cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
c0217812 373 }
374}
eeddc64d 375
e0f4d6b2 376void AliL3Evaluate::FillEffHistosNAIVE()
355debe1 377{
378 //Fill the efficiency histograms.
379
380 cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
381 for(Int_t i=0; i<fGoodGen; i++)
382 {
e0f4d6b2 383 Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
384 Double_t pzg=fGoodTracks[i].pz;
355debe1 385 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
386 //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
387 fNGoodTracksPt->Fill(ptg);
388 fNGoodTracksEta->Fill(dipangle);
389
390 }
391
392 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
393 {
394 AliL3Track *track = fTracks->GetCheckedTrack(k);
395 if(!track) continue;
396 Int_t nHits = track->GetNumberOfPoints();
397 if(nHits < fMinPointsOnTrack) break;
398 if(track->GetPt()<fMinGoodPt) continue;
399 if(fabs(track->GetPseudoRapidity())>0.9) continue;
400
401 fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
402 //Float_t pt=track->GetPt();
403 //fPtRes->Fill((pt-ptg)/ptg*100.);
404 //fNtuppel->Fill(ptg,pt,nHits);
405
406 }
407}
c0217812 408
409void AliL3Evaluate::CalcEffHistos(){
410
411 Stat_t ngood=fNGoodTracksPt->GetEntries();
412 Stat_t nfound=fNFoundTracksPt->GetEntries();
413 Stat_t nfake=fNFakeTracksPt->GetEntries();
414
415 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
416 <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
417 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
418 <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
419 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
420 <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
421 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
422 <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
423
355debe1 424 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
425 <<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
426
427
c0217812 428 fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
429 fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
430 fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
431 fTrackEffPt->SetMaximum(1.4);
432 fTrackEffPt->SetXTitle("P_{T} [GeV]");
433 fTrackEffPt->SetLineWidth(2);
434 fFakeTrackEffPt->SetFillStyle(3013);
435 fTrackEffPt->SetLineColor(4);
436 fFakeTrackEffPt->SetFillColor(2);
437
438 fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
439 fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
440 fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
441 fTrackEffEta->SetMaximum(1.4);
442 fTrackEffEta->SetXTitle("#lambda [degrees]");
443 fTrackEffEta->SetLineWidth(2);
444 fFakeTrackEffEta->SetFillStyle(3013);
445 fTrackEffEta->SetLineColor(4);
446 fFakeTrackEffEta->SetFillColor(2);
447
108615fc 448}
449
450void AliL3Evaluate::Write2File(Char_t *outputfile)
451{
452 //Write histograms to file:
453
e0f4d6b2 454 TFile *of = TFile::Open(outputfile,"RECREATE");
108615fc 455 if(!of->IsOpen())
456 {
457 LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
458 <<"Problems opening rootfile"<<ENDLOG;
459 return;
460 }
461
462 of->cd();
1f79afc0 463 fNtuppel->Write();
108615fc 464 fPtRes->Write();
465 fNGoodTracksPt->Write();
466 fNFoundTracksPt->Write();
467 fNFakeTracksPt->Write();
468 fTrackEffPt->Write();
469 fFakeTrackEffPt->Write();
470 fNGoodTracksEta->Write();
471 fNFoundTracksEta->Write();
472 fNFakeTracksEta->Write();
473 fTrackEffEta->Write();
474 fFakeTrackEffEta->Write();
475
476 of->Close();
108615fc 477
478}
479
a8a3f5c2 480TNtupleD *AliL3Evaluate::CalculateResiduals(Char_t *datapath)
1f79afc0 481{
482
a8a3f5c2 483 TNtupleD *ntuppel=new TNtupleD("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
484 ntuppel->SetDirectory(0);
485
1f79afc0 486 for(int f=fMinSlice; f<=fMaxSlice; f++)
487 {
488 AliL3FileHandler *tfile = new AliL3FileHandler();
489 char fname[256];
a8a3f5c2 490 sprintf(fname,"%s/tracks_tr_%d_0.raw",datapath,f);
1f79afc0 491 if(!tfile->SetBinaryInput(fname)){
492 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
493 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
494 return 0;
495 }
496 fTracks = new AliL3TrackArray();
497 tfile->Binary2TrackArray(fTracks);
498 tfile->CloseBinaryInput();
499 delete tfile;
78127e35 500 printf("Looking in slice %d\n",f);
1f79afc0 501 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
502 {
503
504 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
505 if(!track) continue;
a8a3f5c2 506 if(track->GetNHits() < fMinPointsOnTrack) continue;
1f79afc0 507
508 track->CalculateHelix();
509 UInt_t *hitnum = track->GetHitNumbers();
510 UInt_t id;
511
512 Float_t xyz[3];
513 Int_t padrow;
78127e35 514 for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
1f79afc0 515 {
516 id = hitnum[j];
517 Int_t slice = (id>>25) & 0x7f;
518 Int_t patch = (id>>22) & 0x7;
519 UInt_t pos = id&0x3fffff;
520
78127e35 521 //if(slice!=1) continue;
522
1f79afc0 523 AliL3SpacePointData *points = fClusters[slice][patch];
524
525 if(!points)
526 {
527 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
528 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
529 continue;
530 }
531 if(pos>=fNcl[slice][patch])
532 {
533 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
534 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
535 continue;
536 }
537
538 xyz[0] = points[pos].fX;
539 xyz[1] = points[pos].fY;
540 xyz[2] = points[pos].fZ;
541 padrow = points[pos].fPadRow;
a8a3f5c2 542 AliL3Transform::Global2Local(xyz,slice);
1f79afc0 543
544 Float_t xyz_cross[3];
545 track->GetCrossingPoint(padrow,xyz_cross);
e1a32fa5 546 Double_t beta = track->GetCrossingAngle(padrow);
1f79afc0 547
548 Double_t yres = xyz_cross[1] - xyz[1];
549 Double_t zres = xyz_cross[2] - xyz[2];
550
551 Double_t dipangle = atan(track->GetTgl());
552 ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
553
554 }
555 }
556 if(fTracks)
557 delete fTracks;
558 }
559 return ntuppel;
560}
561
8af8a456 562TNtuple *AliL3Evaluate::EvaluatePoints(Char_t *path)
78127e35 563{
e1a32fa5 564 //Compare points to the exact crossing points of track and padrows.
1b76fd18 565 //The input file to this function, contains the exact clusters calculated
566 //in AliTPC::Hits2ExactClusters.
8af8a456 567
355debe1 568 cout<<"Evaluating points"<<endl;
e1a32fa5 569 TNtuple *ntuppel = new TNtuple("ntuppel","residuals","slice:padrow:resy:resz:zHit:pt");
8af8a456 570 ntuppel->SetDirectory(0);
78127e35 571
8af8a456 572 Char_t filename[1024];
573 sprintf(filename,"%s/alirunfile.root",path);
574 TFile *exfile = TFile::Open(filename);
355debe1 575 if(!exfile)
576 {
8af8a456 577 cerr<<"Error opening rootfile "<<filename<<endl;
355debe1 578 return 0;
579 }
8af8a456 580 gAlice = (AliRun*)exfile->Get("gAlice");
581 if (!gAlice)
582 {
583 LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
584 <<"AliRun object non existing on file"<<ENDLOG;
585 return false;
586 }
587
588 gAlice->GetEvent(0);
5e0f9911 589 AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
e1a32fa5 590
591 //Get the exact clusters from file:
592 AliTPCClustersArray *arr = new AliTPCClustersArray;
593 arr->Setup(param);
594 arr->SetClusterType("AliComplexCluster");
595 char treeName[500];
596 sprintf(treeName,"TreeCExact_%s",param->GetTitle());
294ee877 597 Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
e1a32fa5 598 if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return 0;}
599
355debe1 600 cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
e1a32fa5 601 for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
78127e35 602 {
e1a32fa5 603 //Get the exact clusters for this row:
604 Int_t cursec,currow;
605 AliSegmentID *s = arr->LoadEntry(i);
606 param->AdjustSectorRow(s->GetID(),cursec,currow);
78127e35 607
e1a32fa5 608 AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
609 TClonesArray *clusters = ro->GetArray();
610 int num_of_offline=clusters->GetEntriesFast();
78127e35 611
e1a32fa5 612 //Get the found clusters:
613 Int_t slice,padrow;
494fad94 614 AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
e1a32fa5 615 if(slice<fMinSlice || slice>fMaxSlice) continue;
616 AliL3SpacePointData *points = fClusters[slice][0];
8af8a456 617 if(!points)
618 {
619 cerr<<"AliL3Evaluate::EvalutePoints : Error getting clusters "<<endl;
620 return 0;
621 }
294ee877 622 printf("Checking slice %d padrow %d with %d clusters\n",slice,padrow,num_of_offline);
8af8a456 623 cout<<"There are "<<fNcl[slice][0]<<" clusters here"<<endl;
e1a32fa5 624 for(UInt_t c=0; c<fNcl[slice][0]; c++)
78127e35 625 {
8af8a456 626 if((Int_t)points[c].fPadRow!=padrow) continue;
627 Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
628 Float_t xyz_ex[3];
629 AliL3Transform::Global2Local(xyz_cl,cursec);
630
e1a32fa5 631 for(Int_t m=0; m<num_of_offline; m++)
78127e35 632 {
e1a32fa5 633 AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
634 Int_t mcId = cluster->fTracks[0];
8af8a456 635
294ee877 636 if(mcId <0) continue;
e1a32fa5 637 TParticle *part = gAlice->Particle(mcId);
8af8a456 638 if(points[c].fTrackID[0]!=mcId &&
639 points[c].fTrackID[1]!=mcId &&
640 points[c].fTrackID[2]!=mcId)
e1a32fa5 641 continue;
8af8a456 642
494fad94 643 AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
8af8a456 644
645 //In function AliTPC::Hits2ExactClusters the time offset is not included,
646 //so we have to substract it again here.
647 xyz_ex[2]-=AliL3Transform::GetZOffset();
648
294ee877 649 Float_t resy = xyz_cl[1] - xyz_ex[1];//cluster->GetY()
650 Float_t resz = xyz_cl[2] - xyz_ex[2];//cluster->GetZ()
8af8a456 651
294ee877 652 ntuppel->Fill(slice,padrow,resy,resz,xyz_ex[2],part->Pt());
78127e35 653 }
e1a32fa5 654 }
294ee877 655 arr->ClearRow(cursec,currow);
78127e35 656 }
e1a32fa5 657
78127e35 658 return ntuppel;
659}
355debe1 660
661void AliL3Evaluate::GetCFeff(Char_t *outfile)
78127e35 662{
78127e35 663
355debe1 664 TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","row:ncrossings:nclusters");
665
666 AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
667
668 TPC->SetParam(fParam);
78127e35 669
355debe1 670 Int_t ver = TPC->IsVersion();
494fad94 671 LOG(AliL3Log::kInformational,"AliL3Evaluate::GetCFeff","TPC version")
355debe1 672 <<"TPC version "<<ver<<" found on file"<<ENDLOG;
673
674 Int_t zero=TPC->GetParam()->GetZeroSup();
675
676 Int_t np = gAlice->GetNtrack();
677
678 Int_t crossed,recs;
679 Int_t *count = new Int_t[np]; //np number of particles.
680 Int_t i;
681 Float_t xyz[3];
682 for (i=0; i<np; i++) count[i]=0;
683 for(Int_t sl=fMinSlice; sl<=fMaxSlice; sl++)
78127e35 684 {
355debe1 685 for (i=0; i<=175; i++)
78127e35 686 {
355debe1 687 crossed=0;
688 recs=0;
689 Int_t index = fRowid[sl][i];
690 if (!fDigitsTree->GetEvent(index)) continue;
691 Int_t sec,row;
692 fParam->AdjustSectorRow(fDigits->GetID(),sec,row);
693 fDigits->First();
694 do {
695 Int_t it=fDigits->CurrentRow(), ip=fDigits->CurrentColumn();
696 Short_t dig = fDigits->GetDigit(it,ip);
697
698 if(dig<=fParam->GetZeroSup()) continue;
699 if(it < fParam->GetMaxTBin()-1 && it > 0)
700 if(fDigits->GetDigit(it+1,ip) <= fParam->GetZeroSup()
701 && fDigits->GetDigit(it-1,ip) <= fParam->GetZeroSup())
702 continue;
703
494fad94 704 AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
355debe1 705 if(fParam->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
706 continue;
707
708
709 Int_t idx0=fDigits->GetTrackID(it,ip,0);
710 Int_t idx1=fDigits->GetTrackID(it,ip,1);
711 Int_t idx2=fDigits->GetTrackID(it,ip,2);
712
713 if (idx0>=0 && dig>=zero) count[idx0]+=1;
714 if (idx1>=0 && dig>=zero) count[idx1]+=1;
715 if (idx2>=0 && dig>=zero) count[idx2]+=1;
716 } while (fDigits->Next());
717
718 for (Int_t j=0; j<np; j++)
78127e35 719 {
355debe1 720 if (count[j]>1) //at least two digits at this padrow
494fad94 721 {
722 crossed++;
723 count[j]=0;
724 }
355debe1 725 }
726 AliL3SpacePointData *points = fClusters[sl][0];
494fad94 727 for(UInt_t k=0; k<fNcl[sl][0]; k++)
355debe1 728 {
494fad94 729 if(points[k].fPadRow!=i) continue;
355debe1 730 recs++;
78127e35 731 }
355debe1 732 ntuppel->Fill(i,crossed,recs);
78127e35 733 }
494fad94 734
78127e35 735 }
355debe1 736 delete[] count;
737
738 TFile *file = TFile::Open(outfile,"RECREATE");
739 ntuppel->Write();
740 file->Close();
741
78127e35 742}
355debe1 743