3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliL3StandardIncludes.h"
11 #include <TClonesArray.h>
14 #include <AliSimDigits.h>
16 #include <AliTPCcluster.h>
17 #include <AliTPCClustersArray.h>
18 #include <AliTPCClustersRow.h>
19 #include <AliTPCParam.h>
20 #include <AliComplexCluster.h>
30 #include "AliL3Logging.h"
31 #include "AliL3Transform.h"
32 #include "AliL3SpacePointData.h"
33 #include "AliL3Track.h"
34 #include "AliL3FileHandler.h"
35 #include "AliL3TrackArray.h"
36 #include "AliL3Evaluate.h"
42 /** \class AliL3Evaluate
44 //_____________________________________________________________
47 // Evaluation class for tracking. Plots efficiencies etc..
52 ClassImp(AliL3Evaluate)
54 AliL3Evaluate::AliL3Evaluate()
79 fStandardComparison=kTRUE;
82 AliL3Evaluate::AliL3Evaluate(Char_t *datapath,Int_t min_clusters,Int_t minhits,Double_t minpt,Double_t maxpt,Int_t *slice)
95 sprintf(fPath,"%s",datapath);
96 fMaxFalseClusters = 0.1;
99 fMinPointsOnTrack = min_clusters;
100 fMinHitsFromParticle = minhits;
103 memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
122 fStandardComparison=kTRUE;
125 void AliL3Evaluate::LoadData(Int_t event,Bool_t sp)
128 AliL3FileHandler *clusterfile[36][6];
129 for(Int_t s=fMinSlice; s<=fMaxSlice; s++)
131 for(Int_t p=0; p<AliL3Transform::GetNPatches(); p++)
140 delete fClusters[s][p];
142 clusterfile[s][p] = new AliL3FileHandler();
144 sprintf(fname,"%s/points_%d_%d.raw",fPath,s,patch);
146 sprintf(fname,"%s/points_%d_%d_%d.raw",fPath,event,s,patch);
147 if(!clusterfile[s][p]->SetBinaryInput(fname))
149 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
150 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
151 delete clusterfile[s][p];
152 clusterfile[s][p] = 0;
155 fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
156 clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
157 clusterfile[s][p]->CloseBinaryInput();
163 sprintf(fname,"%s/tracks_%d.raw",fPath,event);
164 AliL3FileHandler *tfile = new AliL3FileHandler();
165 if(!tfile->SetBinaryInput(fname)){
166 LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
167 <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
172 fTracks = new AliL3TrackArray();
173 tfile->Binary2TrackArray(fTracks);
174 tfile->CloseBinaryInput();
180 AliL3Evaluate::~AliL3Evaluate()
182 if(fGoodTracks) delete fGoodTracks;
183 if(fDigitsTree) fDigitsTree->Delete();
184 if(fTracks) delete fTracks;
185 if(fPtRes) delete fPtRes;
186 if(fNGoodTracksPt) delete fNGoodTracksPt;
187 if(fNFoundTracksPt) delete fNFoundTracksPt;
188 if(fNFakeTracksPt) delete fNFakeTracksPt;
189 if(fTrackEffPt) delete fTrackEffPt;
190 if(fFakeTrackEffPt) delete fFakeTrackEffPt;
191 if(fNGoodTracksEta) delete fNGoodTracksEta;
192 if(fNFoundTracksEta) delete fNFoundTracksEta;
193 if(fNFakeTracksEta) delete fNFakeTracksEta;
194 if(fTrackEffEta) delete fTrackEffEta;
195 if(fFakeTrackEffEta) delete fFakeTrackEffEta;
196 if(fMcIndex) delete [] fMcIndex;
197 if(fMcId) delete [] fMcId;
198 if(fNtuppel) delete fNtuppel;
199 if(fNtupleRes) delete fNtupleRes;
202 void AliL3Evaluate::AssignPIDs()
205 LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignPIDs","Track Loop")
206 <<"Assigning pid to the found tracks...."<<ENDLOG;
207 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
209 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
211 if(track->GetNumberOfPoints() < fMinPointsOnTrack)
214 Float_t pid = GetTrackPID(track);
220 void AliL3Evaluate::AssignIDs()
222 //Assign MC id to the tracks.
224 cerr<<"AliL3Evaluate::AssignIDs() : You need to compile with the do_mc flag!"<<endl;
229 LOG(AliL3Log::kDebug,"AliL3Evaluate::AssignIDs","Track Loop")
230 <<"Assigning MC id to the found tracks...."<<ENDLOG;
231 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
233 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
235 if(track->GetNumberOfPoints() < fMinPointsOnTrack) break;
238 Int_t tID = GetMCTrackLabel(track);
241 //cout<<"Found "<<fGoodFound<<" good tracks "<<endl;
244 Float_t AliL3Evaluate::GetTrackPID(AliL3Track *track)
246 track->CalculateHelix();
248 Int_t nc=track->GetNHits();
249 UInt_t *hits = track->GetHitNumbers();
250 Float_t sampleDEdx[159];
251 for (Int_t iHit = 0; iHit < nc; iHit++) {
252 UInt_t hitID = hits[iHit];
253 Int_t iSector = (hitID>>25) & 0x7f;
254 Int_t patch = (hitID>>22) & 0x7;
255 UInt_t position = hitID&0x3fffff;
256 AliL3SpacePointData *points = fClusters[iSector][patch];
257 if(!points) continue;
258 if(position>=fNcl[iSector][patch])
260 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
261 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
264 UChar_t padrow = points[position].fPadRow;
265 Float_t pWidth = AliL3Transform::GetPadPitchWidthLow();
267 pWidth = AliL3Transform::GetPadPitchWidthUp();
268 Float_t corr=1.; if (padrow>63) corr=0.67;
269 sampleDEdx[iHit] = points[position].fCharge/pWidth*corr;
270 Double_t crossingangle = track->GetCrossingAngle(padrow,iSector);
271 Double_t s = sin(crossingangle);
272 Double_t t = track->GetTgl();
273 sampleDEdx[iHit] *= sqrt((1-s*s)/(1+t*t));
278 Int_t swap;//stupid sorting
281 for (i=0; i<nc-1; i++) {
282 if (sampleDEdx[i]<=sampleDEdx[i+1]) continue;
283 Float_t tmp=sampleDEdx[i];
284 sampleDEdx[i]=sampleDEdx[i+1]; sampleDEdx[i+1]=tmp;
289 Double_t low=0.05; Double_t up=0.7;
290 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
292 for (i=nl; i<=nu; i++) trackDEdx += sampleDEdx[i];
293 trackDEdx /= (nu-nl+1);
295 cout<<" PID: "<<nc<<" "<<nl<<" "<<nu<<" "<<trackDEdx<<endl;
299 struct S {Int_t lab; Int_t max;};
300 Int_t AliL3Evaluate::GetMCTrackLabel(AliL3Track *track){
301 //Returns the MCtrackID of the belonging clusters.
302 //If MCLabel < 0, means that track is fake.
303 //Definitions are identical to offline.
305 // - more than 10 percent of clusters are assigned incorrectly
306 // - more than half of the innermost 10% of clusters were assigned incorrectly.
310 Int_t num_of_clusters = track->GetNumberOfPoints();
311 S *s=new S[num_of_clusters];
313 for (i=0; i<num_of_clusters; i++) s[i].lab=s[i].max=0;
314 UInt_t *hitnum = track->GetHitNumbers();
318 for (i=0; i<num_of_clusters; i++)
320 //Tricks to get the clusters belonging to this track:
322 Int_t slice = (id>>25) & 0x7f;
323 Int_t patch = (id>>22) & 0x7;
324 UInt_t pos = id&0x3fffff;
326 AliL3SpacePointData *points = fClusters[slice][patch];
327 if(!points) continue;
328 if(pos>=fNcl[slice][patch])
330 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
331 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
335 //Get the label of the cluster:
336 lab=points[pos].fTrackID[0];
339 for (j=0; j<num_of_clusters; j++)
340 if (s[j].lab==lab || s[j].max==0) break;
346 for (i=0; i<num_of_clusters; i++)
347 if (s[i].max>max) {max=s[i].max; lab=s[i].lab;}
350 return -1; //If most clusters is -1, this is a noise track.
352 cerr<<"AliL3Evaluate::GetMCTrackLabel : Track label negative :"<<lab<<endl;
356 for (i=0; i<num_of_clusters; i++)
359 Int_t slice = (id>>25) & 0x7f;
360 Int_t patch = (id>>22) & 0x7;
361 UInt_t pos = id&0x3fffff;
363 AliL3SpacePointData *points = fClusters[slice][patch];
364 if(!points) continue;
365 if(pos>=fNcl[slice][patch])
367 LOG(AliL3Log::kError,"AliL3Evaluate::GetMCTrackLabel","Clusterarray")
368 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
372 if (abs(points[pos].fTrackID[1]) == lab ||
373 abs(points[pos].fTrackID[2]) == lab ) max++;
377 //Check if more than 10% of the clusters were assigned incorrectly:
378 if (1.-Float_t(max)/num_of_clusters > fMaxFalseClusters)
382 else //Check if at least half of the 10% innermost clusters are assigned correctly.
384 Int_t tail=Int_t(0.10*num_of_clusters);
386 for (i=1; i<=tail; i++)
388 id = hitnum[num_of_clusters - i];
389 Int_t slice = (id>>25) & 0x7f;
390 Int_t patch = (id>>22) & 0x7;
391 UInt_t pos = id&0x3fffff;
393 AliL3SpacePointData *points = fClusters[slice][patch];
394 if(lab == abs(points[pos].fTrackID[0]) ||
395 lab == abs(points[pos].fTrackID[1]) ||
396 lab == abs(points[pos].fTrackID[2])) max++;
398 if (max < Int_t(0.5*tail)) return -lab;
402 #else //If we are running with mc_ids or not
408 void AliL3Evaluate::GetFastClusterIDs(Char_t *path)
410 //Get the MC id of space points in case of using the fast simulator.
412 sprintf(fname,"%s/point_mc.dat",path);
413 FILE *infile = fopen(fname,"r");
418 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
421 fMcId = new Int_t[fNFastPoints];
422 fMcIndex = new UInt_t[fNFastPoints];
424 for(i=0; i<fNFastPoints; i++)
426 if(fscanf(infile,"%d %d",&hitid,&hitmc)==EOF) break;
433 void AliL3Evaluate::CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
435 //Create the histograms
437 LOG(AliL3Log::kInformational,"AliL3Evaluate::CreateHistos","Allocating")
438 <<"Creating histograms..."<<ENDLOG;
440 fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
441 fNtuppel->SetDirectory(0);
442 fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
443 fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
444 fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
445 fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
446 fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
447 fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
449 fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",20,-50,50);
450 fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",20,-50,50);
451 fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",20,-50,50);
452 fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",20,-50,50);
453 fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",20,-50,50);
457 void AliL3Evaluate::GetGoodParticles(Char_t *path,Int_t event,Int_t *padrowrange)
459 //Read the good particles from file. This file should already have been
460 //generated by macro AliTPCComparison.C.
462 Char_t filename[1024];
463 if(event<0 && !padrowrange)
464 sprintf(filename,"%s/good_tracks_tpc",path);
465 else if(event>=0 && !padrowrange)
466 sprintf(filename,"%s/good_tracks_tpc_%d",path,event);
468 sprintf(filename,"%s/good_tracks_tpc_%d_%d_%d",path,event,padrowrange[0],padrowrange[1]);
469 ifstream in(filename);
472 cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
475 Int_t MaxTracks=20000;
477 delete [] fGoodTracks;
479 fGoodTracks = new GoodTrack[MaxTracks];
481 if(fStandardComparison){
482 while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
483 fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
484 fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
486 fGoodTracks[fGoodGen].nhits=-1;
487 fGoodTracks[fGoodGen].sector=-1;
489 if (fGoodGen==MaxTracks)
491 cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
496 while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
497 fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
498 fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y >>fGoodTracks[fGoodGen].z>>fGoodTracks[fGoodGen].nhits>>fGoodTracks[fGoodGen].sector)
501 if (fGoodGen==MaxTracks)
503 cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
510 //has to be modified for fakes.
512 void AliL3Evaluate::FillEffHistos()
516 cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
519 //cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
520 for(Int_t i=0; i<fGoodGen; i++)
522 //cout<<"Checking particle "<<i<<endl;
523 if(!fStandardComparison)
524 if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
525 Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
526 if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
527 Float_t pzg=fGoodTracks[i].pz;
528 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
530 //If we are only considering tracks on one side of the TPC:
535 fNGoodTracksPt->Fill(ptg);
536 fNGoodTracksEta->Fill(dipangle);
539 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
541 AliL3Track *track = fTracks->GetCheckedTrack(k);
543 Int_t nHits = track->GetNumberOfPoints();
544 if(nHits < fMinPointsOnTrack) break;
546 tracklabel = track->GetMCid();
548 if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
550 Float_t pt=track->GetPt();
551 if(tracklabel == fGoodTracks[i].label)
553 fNFoundTracksPt->Fill(ptg);
554 fNFoundTracksEta->Fill(dipangle);
555 fNtuppel->Fill(ptg,pt,nHits);
556 fPtRes->Fill((pt-ptg)/ptg*100.);
560 fNFakeTracksPt->Fill(ptg);
561 fNFakeTracksEta->Fill(dipangle);
563 //fPtRes->Fill((pt-ptg)/ptg*100.);
564 //fNtuppel->Fill(ptg,pt,nHits);
569 //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
573 void AliL3Evaluate::FillEffHistosNAIVE()
575 //Fill the efficiency histograms.
577 cout<<endl<<"Note: Doing NAIVE evaluation "<<endl;
578 for(Int_t i=0; i<fGoodGen; i++)
580 if(!fStandardComparison)
581 if(fGoodTracks[i].nhits < fMinHitsFromParticle) continue;
582 Double_t ptg=TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
583 if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
584 Double_t pzg=fGoodTracks[i].pz;
585 Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
586 //printf("filling particle with pt %f and dipangle %f\n",ptg,dipangle);
587 fNGoodTracksPt->Fill(ptg);
588 fNGoodTracksEta->Fill(dipangle);
592 for(Int_t k=0; k<fTracks->GetNTracks(); k++)
594 AliL3Track *track = fTracks->GetCheckedTrack(k);
596 Int_t nHits = track->GetNumberOfPoints();
597 if(nHits < fMinPointsOnTrack) break;
598 if(track->GetPt()<fMinGoodPt || track->GetPt() > fMaxGoodPt) continue;
599 if(fabs(track->GetPseudoRapidity())>0.9) continue;
601 fNFoundTracksPt->Fill(track->GetPt()); fNFoundTracksEta->Fill(track->GetPseudoRapidity());
602 //Float_t pt=track->GetPt();
603 //fPtRes->Fill((pt-ptg)/ptg*100.);
604 //fNtuppel->Fill(ptg,pt,nHits);
609 void AliL3Evaluate::CalcEffHistos()
612 Stat_t ngood=fNGoodTracksPt->GetEntries();
613 Stat_t nfound=fNFoundTracksPt->GetEntries();
614 Stat_t nfake=fNFakeTracksPt->GetEntries();
616 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
617 <<AliL3Log::kDec<<"There was "<<ngood<<" generated good tracks"<<ENDLOG;
618 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
619 <<AliL3Log::kDec<<"Found "<<nfound<<" tracks"<<ENDLOG;
620 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
621 <<AliL3Log::kDec<<"Integral efficiency is about "<<nfound/ngood*100<<ENDLOG;
622 LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
623 <<AliL3Log::kDec<<"Fake tracks relative is about "<<nfake/ngood*100<<ENDLOG;
624 //LOG(AliL3Log::kInformational,"AliL3Evaluate::FillEffHistos","Efficiency")
625 //<<AliL3Log::kDec<<"Naive efficiency "<<(Double_t)fGoodFound/(Double_t)fGoodGen<<ENDLOG;
627 fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
628 fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
629 fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
630 fTrackEffPt->SetMaximum(1.4);
631 fTrackEffPt->SetXTitle("P_{T} [GeV]");
632 fTrackEffPt->SetLineWidth(2);
633 fFakeTrackEffPt->SetFillStyle(3013);
634 fTrackEffPt->SetLineColor(4);
635 fFakeTrackEffPt->SetFillColor(2);
637 fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
638 fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
639 fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
640 fTrackEffEta->SetMaximum(1.4);
641 fTrackEffEta->SetXTitle("#lambda [degrees]");
642 fTrackEffEta->SetLineWidth(2);
643 fFakeTrackEffEta->SetFillStyle(3013);
644 fTrackEffEta->SetLineColor(4);
645 fFakeTrackEffEta->SetFillColor(2);
649 void AliL3Evaluate::Write2File(Char_t *outputfile)
651 //Write histograms to file:
653 TFile *of = TFile::Open(outputfile,"RECREATE");
656 LOG(AliL3Log::kError,"AliL3Evaluate::Write2File","File Open")
657 <<"Problems opening rootfile"<<ENDLOG;
664 fNGoodTracksPt->Write();
665 fNFoundTracksPt->Write();
666 fNFakeTracksPt->Write();
667 fTrackEffPt->Write();
668 fFakeTrackEffPt->Write();
669 fNGoodTracksEta->Write();
670 fNFoundTracksEta->Write();
671 fNFakeTracksEta->Write();
672 fTrackEffEta->Write();
673 fFakeTrackEffEta->Write();
679 TNtuple *AliL3Evaluate::GetNtuple()
683 fNtupleRes = new TNtuple("ntuppel","Residuals","residual_trans:residual_long:zHit:pt:dipangle:beta:padrow:nHits");
684 fNtupleRes->SetDirectory(0);//Bug in older version of root.
689 void AliL3Evaluate::CalculateResiduals()
692 TNtuple *ntuppel = GetNtuple();
694 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
697 AliL3Track *track = (AliL3Track*)fTracks->GetCheckedTrack(i);
699 if(track->GetNHits() < fMinPointsOnTrack) break;
701 track->CalculateHelix();
702 UInt_t *hitnum = track->GetHitNumbers();
707 for(Int_t j=0; j<track->GetNumberOfPoints()-1; j++)
710 Int_t slice = (id>>25) & 0x7f;
711 Int_t patch = (id>>22) & 0x7;
712 UInt_t pos = id&0x3fffff;
714 //if(slice<18) continue;
716 AliL3SpacePointData *points = fClusters[slice][patch];
720 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
721 <<"No points at slice "<<slice<<" patch "<<patch<<" pos "<<pos<<ENDLOG;
724 if(pos>=fNcl[slice][patch])
726 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Clusterarray")
727 <<AliL3Log::kDec<<"ERROR"<<ENDLOG;
731 xyz[0] = points[pos].fX;
732 xyz[1] = points[pos].fY;
733 xyz[2] = points[pos].fZ;
734 padrow = points[pos].fPadRow;
735 //AliL3Transform::Global2Local(xyz,slice,kTRUE);
736 AliL3Transform::Global2LocHLT(xyz,slice);
739 AliL3Transform::Local2GlobalAngle(&angle,slice);
740 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
742 LOG(AliL3Log::kError,"AliL3Evaluate::CalculateResiduals","Crossing point")
743 <<"Track does not crossing padrow "<<padrow<<" in slice "<<slice<<ENDLOG;
747 Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
748 //AliL3Transform::Global2Local(xyz_cross,slice,kTRUE);
749 AliL3Transform::Global2LocHLT(xyz_cross,slice);
751 Double_t beta = track->GetCrossingAngle(padrow,slice);
753 Double_t yres = xyz_cross[1] - xyz[1];
754 Double_t zres = xyz_cross[2] - xyz[2];
755 Double_t dipangle = atan(track->GetTgl());
756 ntuppel->Fill(yres,zres,xyz_cross[2],track->GetPt(),dipangle,beta,padrow,track->GetNumberOfPoints());
762 enum tagprimary {kPrimaryCharged = 0x4000};
763 void AliL3Evaluate::EvaluatePoints(Char_t *rootfile,Char_t *exactfile,Char_t *tofile,Int_t nevent,Bool_t offline,Bool_t sp)
765 //Compare points to the exact crossing points of track and padrows.
766 //The input file to this function, contains the exact clusters calculated
767 //in AliTPC::Hits2ExactClusters.
770 cerr<<"AliL3Evaluate::EvaluatePoints : Compile with do_mc flag!"<<endl;
773 cout<<"Evaluating points"<<endl;
774 TNtuple *ntuppel = new TNtuple("ntuppel_res","Cluster properties",
775 "slice:padrow:charge:resy:resz:zHit:pt:beta:sigmaY2:sigmaZ2:psigmaY2:psigmaZ2");
776 ntuppel->SetDirectory(0);
778 TNtuple *ntuppel2 = new TNtuple("ntuppel_eff","Efficiency","slice:padrow:nfound:ngen");
779 ntuppel2->SetDirectory(0);
781 TFile *exfile = TFile::Open(rootfile);
784 cerr<<"Error opening rootfile "<<rootfile<<endl;
787 gAlice = (AliRun*)exfile->Get("gAlice");
790 LOG(AliL3Log::kError,"AliL3Evaluate::InitMC","gAlice")
791 <<"AliRun object non existing on file"<<ENDLOG;
795 AliTPCParam *param = (AliTPCParam*)exfile->Get(AliL3Transform::GetParamName());
797 TFile *exact = TFile::Open(exactfile);
800 cerr<<"AliL3Evaluate::EvaluatePoints : Problems opening file :"<<exactfile<<endl;
804 AliStack *astack=gAlice->Stack();
806 AliTPCClustersArray *arr=0;
807 for(Int_t event=0; event<nevent; event++)
813 Int_t nparticles = gAlice->GetEvent(event);
814 Int_t nprimaries = 0;//FindPrimaries(nparticles);
815 cout<<"Event "<<event<<" had "<<nparticles<<" particles and "<<nprimaries<<" primaries"<<endl;
818 //Get the exact clusters from file:
819 AliTPCClustersArray *arr = new AliTPCClustersArray;
821 arr->SetClusterType("AliComplexCluster");
823 sprintf(treeName,"TreeCExact_%s_%d",param->GetTitle(),event);
824 Bool_t clusterok = arr->ConnectTree(treeName);//Segment Tree (for offline clusters)
825 if(!clusterok) {printf("AliL3Evaluate::EvaluatePoints : Error in clusterloading\n"); return;}
827 //cout<<"Entering loop with "<<(Int_t)arr->GetTree()->GetEntries()<<endl;
828 for(Int_t i=0; i<arr->GetTree()->GetEntries(); i++)
830 //Get the exact clusters for this row:
832 AliSegmentID *s = arr->LoadEntry(i);
833 param->AdjustSectorRow(s->GetID(),cursec,currow);
835 AliTPCClustersRow *ro = (AliTPCClustersRow *)arr->GetRow(cursec,currow);
836 TClonesArray *clusters = ro->GetArray();
837 int num_of_offline=clusters->GetEntriesFast();
839 //Get the found clusters:
841 AliL3Transform::Sector2Slice(slice,padrow,cursec,currow);
842 if(slice < fMinSlice) continue;
843 if(slice > fMaxSlice) break;
845 Int_t patch = AliL3Transform::GetPatch(padrow);
848 AliL3SpacePointData *points = fClusters[slice][patch];
852 //cout<<"Slice "<<slice<<" padrow "<<padrow<<" has "<<num_of_offline<<" clusters "<<endl;
853 Int_t clustercount=0;
855 for(Int_t m=0; m<num_of_offline; m++)
857 AliComplexCluster *cluster = (AliComplexCluster *)clusters->UncheckedAt(m);
859 Int_t mcId = cluster->GetTrack(0);
861 Int_t mcId = cluster->fTracks[0];
863 if(mcId <0) continue;
866 if(cluster->GetY() < 1 || cluster->GetY() > AliL3Transform::GetNPads(padrow) - 2 ||
867 cluster->GetX() < 1 || cluster->GetX() > AliL3Transform::GetNTimeBins() - 2)
870 if(cluster->fY < 1 || cluster->fY > AliL3Transform::GetNPads(padrow) - 2 ||
871 cluster->fX < 1 || cluster->fX > AliL3Transform::GetNTimeBins() - 2)
877 AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->GetY(),cluster->GetX());
879 AliL3Transform::Raw2Local(xyz_ex,cursec,currow,cluster->fY,cluster->fX);
881 //In function AliTPC::Hits2ExactClusters the time offset is not included,
882 //so we have to substract it again here.
884 xyz_ex[2]-=AliL3Transform::GetZOffset();
886 xyz_ex[2]+=AliL3Transform::GetZOffset();
889 if(param->GetPadRowRadii(cursec,currow)<230./250.*fabs(xyz_ex[2]))
892 TParticle *part = astack->Particle(mcId);
895 if(part->Pt() < fMinGoodPt) continue;
897 //Dont take secondaries, because in width calculation we assume primaries:
898 //if(!(part->TestBit(kPrimaryCharged))) continue;
899 if(part->GetFirstMother()>=0) continue;
902 for(UInt_t c=0; c<fNcl[slice][patch]; c++)
904 if((Int_t)points[c].fPadRow!=padrow) continue;
905 Float_t xyz_cl[3] = {points[c].fX,points[c].fY,points[c].fZ};
908 AliL3Transform::Global2Local(xyz_cl,cursec);
911 if(points[c].fTrackID[0] != mcId &&
912 points[c].fTrackID[1] != mcId &&
913 points[c].fTrackID[2] != mcId)
917 Float_t resy = xyz_cl[1] - xyz_ex[1];
918 Float_t resz = xyz_cl[2] - xyz_ex[2];
921 Int_t charge = (Int_t)points[c].fCharge;
922 Float_t beta = GetCrossingAngle(part,slice,padrow,xyz_ex);
923 Double_t tanl = xyz_ex[2]/sqrt(xyz_ex[0]*xyz_ex[0]+xyz_ex[1]*xyz_ex[1]);
924 Float_t psigmaY2 = AliL3Transform::GetParSigmaY2(padrow,xyz_ex[2],beta);
925 Float_t psigmaZ2 = AliL3Transform::GetParSigmaZ2(padrow,xyz_ex[2],tanl);
926 Float_t sigmaY2 = points[c].fSigmaY2;
927 Float_t sigmaZ2 = points[c].fSigmaZ2;
928 ntuppel->Fill(slice,padrow,charge,resy,resz,xyz_ex[2],part->Pt(),beta,sigmaY2,sigmaZ2,psigmaY2,psigmaZ2);
930 clustercount=tempcount;
932 ntuppel2->Fill(slice,padrow,clustercount,crosscount);
933 arr->ClearRow(cursec,currow);
939 TFile *ofile = TFile::Open(tofile,"RECREATE");
947 void AliL3Evaluate::GetCFeff(Char_t *path,Char_t *outfile,Int_t nevent,Bool_t sp)
949 //Evaluate the cluster finder efficiency.
952 cerr<<"AliL3Evaluate::GetCFeff : Compile with do_mc flag"<<endl;
955 TNtuple *ntuppel = new TNtuple("ntuppel","Cluster finder efficiency","slice:row:ncrossings:nclusters");
956 ntuppel->SetDirectory(0);
958 Char_t filename[1024];
959 sprintf(filename,"%s/alirunfile.root",path);
960 TFile *rfile = TFile::Open(filename);
961 gAlice = (AliRun*)rfile->Get("gAlice");
963 AliStack *astack=gAlice->Stack();
965 AliTPCParam *param = (AliTPCParam*)rfile->Get(AliL3Transform::GetParamName());
967 Int_t zero=param->GetZeroSup();
969 sprintf(filename,"%s/digitfile.root",path);
970 TFile *dfile = TFile::Open(filename);
972 for(Int_t event=0; event<nevent; event++)
976 gAlice->GetEvent(event);
977 Int_t np = astack->GetNtrack();
978 cout<<"Processing event "<<event<<" with "<<np<<" particles "<<endl;
980 sprintf(filename,"TreeD_75x40_100x60_150x60_%d",event);
981 TTree *TD=(TTree*)gDirectory->Get(filename);
982 AliSimDigits da, *digits=&da;
983 TD->GetBranch("Segment")->SetAddress(&digits);
985 Int_t crossed=0,recs=0;
986 Int_t *count = new Int_t[np]; //np number of particles.
989 for (i=0; i<np; i++) count[i]=0;
993 for(Int_t i=0; i<(Int_t)TD->GetEntries(); i++)
996 if (!TD->GetEvent(i)) continue;
997 param->AdjustSectorRow(digits->GetID(),sec,row);
998 AliL3Transform::Sector2Slice(sl,sr,sec,row);
999 if(sl < fMinSlice) continue;
1000 if(sl > fMaxSlice) break;
1001 cout<<"Processing slice "<<sl<<" row "<<sr<<endl;
1004 Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
1005 Short_t dig = digits->GetDigit(it,ip);
1007 if(dig<=param->GetZeroSup()) continue;
1008 AliL3Transform::Raw2Local(xyz,sec,row,ip,it);
1009 if(param->GetPadRowRadii(sec,row)<230./250.*fabs(xyz[2]))
1012 Int_t idx0=digits->GetTrackID(it,ip,0);
1013 Int_t idx1=digits->GetTrackID(it,ip,1);
1014 Int_t idx2=digits->GetTrackID(it,ip,2);
1016 if (idx0>=0 && dig>=zero) count[idx0]+=1;
1017 if (idx1>=0 && dig>=zero) count[idx1]+=1;
1018 if (idx2>=0 && dig>=zero) count[idx2]+=1;
1019 } while (digits->Next());
1020 for (Int_t j=0; j<np; j++)
1022 TParticle *part = astack->Particle(j);
1023 if(part->Pt() < fMinGoodPt) continue;
1024 if(part->GetFirstMother() >= 0) continue;
1025 if (count[j]>1) //at least two digits at this padrow
1032 Int_t patch = AliL3Transform::GetPatch(sr);
1035 AliL3SpacePointData *points = fClusters[sl][patch];
1038 for(UInt_t k=0; k<fNcl[sl][patch]; k++)
1040 if(points[k].fPadRow!=sr) continue;
1043 ntuppel->Fill(sl,sr,crossed,recs);
1049 TFile *file = TFile::Open(outfile,"RECREATE");
1058 Float_t AliL3Evaluate::GetCrossingAngle(TParticle *part,Int_t slice,Int_t /*padrow*/,Float_t *xyz)
1060 //Calculate the padrow crossing angle of the particle
1062 Double_t kappa = AliL3Transform::GetBField()*AliL3Transform::GetBFact()/part->Pt();
1064 Double_t radius = 1/fabs(kappa);
1065 if(part->GetPdgCode() > 0) kappa = -kappa;
1067 Float_t angl[1] = {part->Phi()};
1069 AliL3Transform::Global2LocalAngle(angl,slice);
1071 Double_t charge = -1.*kappa;
1073 Double_t trackPhi0 = angl[0] + charge*0.5*AliL3Transform::Pi()/fabs(charge);
1077 Double_t xc = x0 - radius * cos(trackPhi0);
1078 Double_t yc = y0 - radius * sin(trackPhi0);
1080 Double_t tangent[2];
1081 tangent[0] = -1.*(xyz[1] - yc)/radius;
1082 tangent[1] = (xyz[0] - xc)/radius;
1084 Double_t perp_padrow[2] = {1,0}; //locally in slice
1086 Double_t cos_beta = fabs(tangent[0]*perp_padrow[0] + tangent[1]*perp_padrow[1]);
1087 if(cos_beta > 1) cos_beta=1;
1088 return acos(cos_beta);
1091 Int_t AliL3Evaluate::FindPrimaries(Int_t nparticles)
1094 Double_t vertcut = 0.001;
1095 Double_t decacut = 3.;
1096 Double_t timecut = 0.;
1098 AliStack *astack=gAlice->Stack();
1099 TParticle * part = astack->Particle(0);
1100 Double_t xori = part->Vx();
1101 Double_t yori = part->Vy();
1102 Double_t zori = part->Vz();
1103 for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
1105 part = astack->Particle(iprim);
1106 const char * xxx=strstr(part->GetName(),"XXX");
1109 TParticlePDG *ppdg = part->GetPDG();
1110 if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
1112 Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
1113 if(dist>vertcut)continue; // cut on the vertex
1115 if(part->T()>timecut)continue;
1117 Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
1118 if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
1120 Bool_t prmch = kTRUE; // candidate primary track
1121 Int_t fidau=part->GetFirstDaughter(); // cut on daughters
1125 lasdau=part->GetLastDaughter();
1126 ndau=lasdau-fidau+1;
1129 for(Int_t j=fidau;j<=lasdau;j++){
1130 TParticle *dau=astack->Particle(j);
1131 Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
1132 if(distd<decacut)prmch=kFALSE; // eliminate if the decay is near the vertex
1138 part->SetBit(kPrimaryCharged);