1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //-------------------------------------------------------------------
17 // Base class to test the extrapolation performance from TPC to outer
18 // detectors. Several member functions AddTracks() with different
19 // arguments can be called to execute extrapolation and several
20 // THnSparse are filled with residuls and basic track information
22 // Anthor: R.Ma, M.Ivanov 02/04/2011
23 //--------------------------------------------------------------------
26 #include "AliTrackComparison.h"
27 #include "AliExternalTrackParam.h"
28 #include "AliTrackReference.h"
29 #include "AliTracker.h"
30 #include "AliTrackPointArray.h"
31 #include "THnSparse.h"
32 #include "TParticle.h"
34 #include "TParticlePDG.h"
35 #include "TParticle.h"
36 #include "TTreeStream.h"
39 #include "TTreeStream.h"
41 ClassImp(AliTrackComparison)
43 //________________________________________________________________________
44 AliTrackComparison::AliTrackComparison()
70 // Default constructor
72 for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
76 //________________________________________________________________________
77 AliTrackComparison::AliTrackComparison(const Text_t *name, const Text_t *title)
103 // Non default cosntructor
105 for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
108 //________________________________________________________________________
109 AliTrackComparison::AliTrackComparison(const AliTrackComparison& comp)
112 , fLowBinDY(comp.fLowBinDY)
113 , fUpBinDY(comp.fUpBinDY)
114 , fLowBinDZ(comp.fLowBinDZ)
115 , fUpBinDZ(comp.fUpBinDZ)
116 , fLowBinDSnp(comp.fLowBinDSnp)
117 , fUpBinDSnp(comp.fUpBinDSnp)
118 , fLowBinDTheta(comp.fLowBinDTheta)
119 , fUpBinDTheta(comp.fUpBinDTheta)
120 , fLowBin1Pt(comp.fLowBin1Pt)
121 , fUpBin1Pt(comp.fUpBin1Pt)
122 , fLowBin1PtLoss(comp.fLowBin1PtLoss)
123 , fUpBin1PtLoss(comp.fUpBin1PtLoss)
124 , fNBinsDY(comp.fNBinsDY)
125 , fNBinsDZ(comp.fNBinsDZ)
126 , fNBinsDSnp(comp.fNBinsDSnp)
127 , fNBinsDTheta(comp.fNBinsDTheta)
128 , fNBins1Pt(comp.fNBins1Pt)
129 , fNBins1PtLoss(comp.fNBins1PtLoss)
130 , fLayerID(comp.fLayerID)
131 , fFillAll(comp.fFillAll)
132 , fNCombineBin(comp.fNCombineBin)
138 for (Int_t i=0;i<6; i++) fResolHisto[i]=comp.fResolHisto[i];
141 //________________________________________________________________________
142 AliTrackComparison& AliTrackComparison::operator=(const AliTrackComparison& comp)
148 TNamed::operator=(comp);
153 //________________________________________________________________________
154 void AliTrackComparison::Init(){
156 //Initilized the output histograms
162 //________________________________________________________________________
163 AliTrackComparison::~AliTrackComparison(){
169 //________________________________________________________________________
170 void AliTrackComparison::Analyze() {
176 //________________________________________________________________________
177 Int_t AliTrackComparison::AddTracks(AliTrackReference *ref0, AliTrackReference *ref1, Double_t mass, Int_t charge)
179 // Make track param out of track reference
180 // Test propagation from ref0 to ref1
181 // Fill the THnSparse with ref0 and ref1
183 AliExternalTrackParam *param0 = 0;
184 AliExternalTrackParam *param1 = 0;
186 param0=MakeTrack(ref0,charge);
187 param1=MakeTrack(ref1,charge);
189 if (!param0 || !param1) return 0;
191 Double_t tr1Pt = param0->GetSigned1Pt();
193 if(!PropagateToPoint(param0,param1,mass)) return 0;
195 FillHistos(param0,param1,tr1Pt);
199 //________________________________________________________________________
200 Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
202 //Test propagation from param0 to param1
203 //Fill the THnSparse with param0 and param1
205 Double_t tr1Pt=param0->GetSigned1Pt();
206 if( !PropagateToPoint(param0,param1,mass)) return 0;
207 FillHistos(param0,param1,tr1Pt);
211 //________________________________________________________________________
212 Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0, AliTrackPoint *point1, Double_t mass, Double_t energy, Double_t *vxyz)
214 //Test propagation from param0 to point1
215 //This function is usually called in real data analysis when only the
216 //position of the track point is known. In this case, the angles
217 //in the track param are not the angles of the track momentum, but position.
218 //Only the first two and last two THnSparse are meaninglful and should be
219 //filled. Set this via SetFillAll(kFALSE)
221 Double_t gPos[3]= {point1->GetX(),point1->GetY(),point1->GetZ()};
223 Double_t pos[3], pxyz[3];
224 for(Int_t i=0; i<3; i++)
225 pos[i] = gPos[i]-vxyz[i];
226 Double_t R = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
227 for(Int_t i=0; i<3; i++) pxyz[i]= energy*pos[i]/R;
230 for (Int_t i=0; i<21;i++) cv[i]=0;
231 AliExternalTrackParam * param1 = new AliExternalTrackParam(gPos,pxyz,cv,param0->Charge());
233 if(!param1) return 0;
234 Double_t tr1Pt = param0->GetSigned1Pt();
235 if(!PropagateToPoint(param0,param1,mass)) return 0;
237 FillHistos(param0,param1,tr1Pt);
241 //________________________________________________________________________
242 Bool_t AliTrackComparison::PropagateToPoint(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
245 //Extrapolate is performed
247 Double_t radius = param1->GetX();
248 param0->Rotate(param1->GetAlpha());
249 AliTracker::PropagateTrackToBxByBz(param0, radius+fStep, mass, fStep, kFALSE,0.99,-1);
250 Bool_t isOK = param0->PropagateTo(radius,AliTracker::GetBz());
254 //________________________________________________________________________
255 AliExternalTrackParam * AliTrackComparison::MakeTrack(const AliTrackReference* ref, Int_t charge)
258 // Make track out of the track ref
259 // the covariance matrix - equal 0 - starting from ideal MC position
260 Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
261 Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
263 for (Int_t i=0; i<21;i++) cv[i]=0;
264 AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
268 //________________________________________________________________________
269 Long64_t AliTrackComparison::Merge(TCollection *const li) {
271 //Merge the comparison objects
273 TIterator* iter = li->MakeIterator();
274 AliTrackComparison* comp = 0;
275 TString strName(GetName());
276 while ((comp = (AliTrackComparison*)iter->Next())) {
277 if (!comp->InheritsFrom(AliTrackComparison::Class())) {
280 if (strName.CompareTo(comp->GetName())!=0) return -1;
281 // add histograms here...
287 //________________________________________________________________________
288 void AliTrackComparison::Add(AliTrackComparison *const comp){
292 for (Int_t i=0;i<6;i++){
295 THnSparse * h0 = (THnSparse*)fResolHisto[i];
296 THnSparse * h1 = (THnSparse*)comp->GetHnSparse(i);
297 if (h0&&h1) h0->Add(h1);
302 //________________________________________________________________________
303 void AliTrackComparison::MakeHistos()
306 //Called in Init() to initialize histograms
308 Double_t xminTrack[5], xmaxTrack[5];
311 TString axisTitle[5];
312 TString hisNames[6]={"DeltaY","DeltaZ","DeltaSnp","DeltaTheta","Delta1Pt","1PtLoss"};
313 Double_t lowBins[6]={fLowBinDY, fLowBinDZ, fLowBinDSnp, fLowBinDTheta, fLowBin1Pt, fLowBin1PtLoss};
314 Double_t upBins[6]={fUpBinDY, fUpBinDZ, fUpBinDSnp, fUpBinDTheta, fUpBin1Pt, fUpBin1PtLoss};
315 Int_t nBins[6] = {fNBinsDY, fNBinsDZ, fNBinsDSnp, fNBinsDTheta, fNBins1Pt, fNBins1PtLoss};
317 axisName[0] ="#Delta";
318 axisTitle[0] ="#Delta";
321 xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
322 axisName[1] ="tanTheta";
323 axisTitle[1] ="tan(#Theta)";
326 xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
328 axisTitle[2] ="#phi";
331 xminTrack[3] =-3.; xmaxTrack[3]=3.; // 0.33 GeV cut
333 axisTitle[3] ="1/pt";
336 xminTrack[4] =0; xmaxTrack[4]=22;
337 axisName[4] ="LayerID";
338 axisTitle[4] ="LayerID";
340 for (Int_t ihis=0; ihis<6; ihis++){
341 // modify ranges for bin 0
343 binsTrack[0]=nBins[ihis];
344 xminTrack[0]=lowBins[ihis];
345 xmaxTrack[0]=upBins[ihis];
346 fResolHisto[ihis] = new THnSparseF(hisNames[ihis],hisNames[ihis], 5, binsTrack,xminTrack, xmaxTrack);
347 for (Int_t ivar=0;ivar<5;ivar++){
348 fResolHisto[ihis]->GetAxis(ivar)->SetName(axisName[ivar].Data());
349 fResolHisto[ihis]->GetAxis(ivar)->SetTitle(axisTitle[ivar].Data());
354 //________________________________________________________________________
355 void AliTrackComparison::FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
357 //Fill the THnSparse.
358 //In case of not filling all, only the 4 out of 6 THnSparse are filed
360 Double_t dY = param1->GetY()-param0->GetY(); //Residual in Y
361 Double_t dZ = param1->GetZ()-param0->GetZ(); //Residual in Z
362 Double_t dSnp = param1->GetSnp()-param0->GetSnp(); //Residual in sin(phi)
363 Double_t dTheta = param1->Theta()-param0->Theta(); //Residual in Theta
364 Double_t d1Pt = param1->GetSigned1Pt()-param0->GetSigned1Pt(); //Residual in 1/pT
365 Double_t d1PtLoss = param0->GetSigned1Pt()-tr1Pt; //Corrected energy loss
367 Double_t globalPos[3];
368 param1->GetXYZ(globalPos);
369 //printf("param1: Atan(y,x)=%5.3f | phi=%5.3f\n",TMath::ATan2(globalPos[1],globalPos[0]),param1->Phi());
370 //printf("trkP=%5.3f | param0=%5.3f | param1=%5.3f\n",trP,param0->GetP(),param1->GetP());
372 Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
374 for(Int_t j=0; j<6; j++)
376 if(!fFillAll && j<4 && j>1) continue;
377 argu[j][0] = residual[j];
378 argu[j][1] = param1->GetTgl();
379 argu[j][2] = TMath::ATan2(globalPos[1],globalPos[0]);
380 argu[j][3] = param1->GetSigned1Pt();
381 argu[j][4] = fLayerID;
382 fResolHisto[j]->Fill(argu[j]);
386 //________________________________________________________________________
387 void AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
389 // make a distortion map out of the residual histogram
390 // Results are written to the debug streamer - pcstream
392 // pcstream - file to write the tree
394 // refX - track matching reference X
395 // type - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
397 // OBJ: TAxis #Delta #Delta
398 // OBJ: TAxis tanTheta tan(#Theta)
399 // OBJ: TAxis phi #phi
400 // OBJ: TAxis 1/pt 1/pt
401 // OBJ: TAxis LayerID LayerID
402 // marian.ivanov@cern.ch
404 //Double_t refX=365; //temporary
405 Int_t run=0; //temporary
406 TString hname = Form("%s_%d",GetName(),type);
407 THnSparse * his0= GetHnSparse(type);
411 dsName+=Form("Debug_%d.root",type);
412 printf(" Create debug streamer \n");
413 dsName.ReplaceAll(" ","");
414 TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());
417 const Int_t kMinEntries=10;
418 Double_t bz=AliTrackerBase::GetBz();
419 Int_t idim[4]={0,1,2,3};
423 Int_t nbins3=his0->GetAxis(3)->GetNbins();
424 Int_t first3=his0->GetAxis(3)->GetFirst();
425 Int_t last3 =his0->GetAxis(3)->GetLast();
427 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - 1/pt
428 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-fNCombineBin,1),TMath::Min(ibin3+fNCombineBin,nbins3));
429 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
430 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
432 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
433 Int_t first2 = his3->GetAxis(2)->GetFirst();
434 Int_t last2 = his3->GetAxis(2)->GetLast();
436 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
437 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-fNCombineBin,1),TMath::Min(ibin2+fNCombineBin,nbins2));
438 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
439 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
440 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
441 Int_t first1 = his2->GetAxis(1)->GetFirst();
442 Int_t last1 = his2->GetAxis(1)->GetLast();
443 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - tan(theta)
445 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
446 Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);
448 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1+fNCombineBin,nbins1));
449 if (TMath::Abs(x1)<(fNCombineBin+0.5)*binWidth){
450 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1,nbins1));
451 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+fNCombineBin,nbins1));
453 TH1 * hisDelta = his2->Projection(0);
455 Double_t entries = hisDelta->GetEntries();
456 Double_t mean=0, rms=0;
457 if (entries>kMinEntries){
458 mean = hisDelta->GetMean();
459 rms = hisDelta->GetRMS();
461 Double_t sector = 9.*x2/TMath::Pi();
462 if (sector<0) sector+=18;
463 Double_t dsec = sector-Int_t(sector)-0.5;
464 Double_t quantiles[5]={0};
465 Double_t mean75=0, rms75=0;
466 Double_t prob[5]={0, 0.25,0.5,0.75,1};
467 if (entries>kMinEntries){
468 hisDelta->GetQuantiles(5,quantiles,prob);
469 if(x3>0)hisDelta->GetXaxis()->SetRangeUser(quantiles[0], quantiles[3]);
470 else hisDelta->GetXaxis()->SetRangeUser(quantiles[1], quantiles[4]);
471 rms75=hisDelta->GetRMS();
472 mean75=hisDelta->GetMean();
475 Double_t pt = TMath::Abs(x3);
477 (*pcstream)<<"distortion"<<
480 "theta="<<x1<< // tan(theta)
483 "mpt="<<x3<< // signed 1/pt
485 "entries="<<entries<< // entries
486 "mean="<<mean<< // normal mean
488 "q25="<<quantiles[1]<< // quantiles of distribution - importnat for asymetric distibutions
489 "q50="<<quantiles[2]<< // quantiles of distribution - importnat for asymetric distibutions
490 "q75="<<quantiles[3]<< // quantiles of distribution - importnat for asymetric distibutions
491 "mean75="<<mean75<< // mean of the truncated distribution 0-75%
492 "rms75="<<rms75<< // rms of the truncated distibution up to 0-75%
494 "rms="<<rms<< // normal rms
495 "refX="<<refX<< // track matching refernce plane
496 "type="<<type<< // tpye of residuals
497 "sector="<<sector<< // sector number according to TPC standard
501 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
508 printf(" Finished! Close debug streamer \n");
514 //TFile f("outFile.root");
515 //TList * list = f.Get("jthaeder_OutterTracking");
516 TFile f("/u/alice-st/marr/EnergyCorrection/TrainAnalysis/trunk/marr_TrackingTest.root");
517 TList * list = f.Get("marr_TrackingTest");
519 AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
520 AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
521 AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
522 AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
525 compEMCAL-> MakeDistortionMap(refX,0);
526 compEMCAL-> MakeDistortionMap(refX,1);
527 compEMCAL-> MakeDistortionMap(refX,4);
528 compEMCALEl-> MakeDistortionMap(refX,0);
529 compEMCALEl-> MakeDistortionMap(refX,1);
530 compEMCALEl-> MakeDistortionMap(refX,4);
532 compTOF->SetNCombineBin(3)
533 compTOF->MakeDistortionMap(365,4);
534 compTOF->MakeDistortionMap(365,0);
536 TFile f("emcalDistortion.root");
539 // ideal case - no mislaignment dy only due energy loss correction
540 TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
542 TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
544 TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
546 TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10"); // pions - bias +-1 %
548 TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10"); // electrons - bias +-20 %