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[3]);
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 if (!fResolHisto) return;
293 for (Int_t i=0;i<6;i++){
294 THnSparse * h0 = (THnSparse*)fResolHisto[i];
295 THnSparse * h1 = (THnSparse*)comp->GetHnSparse(i);
296 if (h0&&h1) h0->Add(h1);
301 //________________________________________________________________________
302 void AliTrackComparison::MakeHistos()
305 //Called in Init() to initialize histograms
307 Double_t xminTrack[5], xmaxTrack[5];
310 TString axisTitle[5];
311 TString hisNames[6]={"DeltaY","DeltaZ","DeltaSnp","DeltaTheta","Delta1Pt","1PtLoss"};
312 Double_t lowBins[6]={fLowBinDY, fLowBinDZ, fLowBinDSnp, fLowBinDTheta, fLowBin1Pt, fLowBin1PtLoss};
313 Double_t upBins[6]={fUpBinDY, fUpBinDZ, fUpBinDSnp, fUpBinDTheta, fUpBin1Pt, fUpBin1PtLoss};
314 Int_t nBins[6] = {fNBinsDY, fNBinsDZ, fNBinsDSnp, fNBinsDTheta, fNBins1Pt, fNBins1PtLoss};
316 axisName[0] ="#Delta";
317 axisTitle[0] ="#Delta";
320 xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
321 axisName[1] ="tanTheta";
322 axisTitle[1] ="tan(#Theta)";
325 xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
327 axisTitle[2] ="#phi";
330 xminTrack[3] =-3.; xmaxTrack[3]=3.; // 0.33 GeV cut
332 axisTitle[3] ="1/pt";
335 xminTrack[4] =0; xmaxTrack[4]=22;
336 axisName[4] ="LayerID";
337 axisTitle[4] ="LayerID";
339 for (Int_t ihis=0; ihis<6; ihis++){
340 // modify ranges for bin 0
342 binsTrack[0]=nBins[ihis];
343 xminTrack[0]=lowBins[ihis];
344 xmaxTrack[0]=upBins[ihis];
345 fResolHisto[ihis] = new THnSparseF(hisNames[ihis],hisNames[ihis], 5, binsTrack,xminTrack, xmaxTrack);
346 for (Int_t ivar=0;ivar<5;ivar++){
347 fResolHisto[ihis]->GetAxis(ivar)->SetName(axisName[ivar].Data());
348 fResolHisto[ihis]->GetAxis(ivar)->SetTitle(axisTitle[ivar].Data());
353 //________________________________________________________________________
354 void AliTrackComparison::FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
356 //Fill the THnSparse.
357 //In case of not filling all, only the 4 out of 6 THnSparse are filed
359 Double_t dY = param1->GetY()-param0->GetY(); //Residual in Y
360 Double_t dZ = param1->GetZ()-param0->GetZ(); //Residual in Z
361 Double_t dSnp = param1->GetSnp()-param0->GetSnp(); //Residual in sin(phi)
362 Double_t dTheta = param1->Theta()-param0->Theta(); //Residual in Theta
363 Double_t d1Pt = param1->GetSigned1Pt()-param0->GetSigned1Pt(); //Residual in 1/pT
364 Double_t d1PtLoss = param0->GetSigned1Pt()-tr1Pt; //Corrected energy loss
366 Double_t globalPos[3];
367 param1->GetXYZ(globalPos);
368 //printf("param1: Atan(y,x)=%5.3f | phi=%5.3f\n",TMath::ATan2(globalPos[1],globalPos[0]),param1->Phi());
369 //printf("trkP=%5.3f | param0=%5.3f | param1=%5.3f\n",trP,param0->GetP(),param1->GetP());
371 Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
373 for(Int_t j=0; j<6; j++)
375 if(!fFillAll && j<4 && j>1) continue;
376 argu[j][0] = residual[j];
377 argu[j][1] = param1->GetTgl();
378 argu[j][2] = TMath::ATan2(globalPos[1],globalPos[0]);
379 argu[j][3] = param1->GetSigned1Pt();
380 argu[j][4] = fLayerID;
381 fResolHisto[j]->Fill(argu[j]);
385 //________________________________________________________________________
386 void AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
388 // make a distortion map out of the residual histogram
389 // Results are written to the debug streamer - pcstream
391 // pcstream - file to write the tree
393 // refX - track matching reference X
394 // type - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
396 // OBJ: TAxis #Delta #Delta
397 // OBJ: TAxis tanTheta tan(#Theta)
398 // OBJ: TAxis phi #phi
399 // OBJ: TAxis 1/pt 1/pt
400 // OBJ: TAxis LayerID LayerID
401 // marian.ivanov@cern.ch
403 //Double_t refX=365; //temporary
404 Int_t run=0; //temporary
405 TString hname = Form("%s_%d",GetName(),type);
406 THnSparse * his0= GetHnSparse(type);
410 dsName+=Form("Debug_%d.root",type);
411 printf(" Create debug streamer \n");
412 dsName.ReplaceAll(" ","");
413 TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());
416 const Int_t kMinEntries=10;
417 Double_t bz=AliTrackerBase::GetBz();
418 Int_t idim[4]={0,1,2,3};
422 Int_t nbins3=his0->GetAxis(3)->GetNbins();
423 Int_t first3=his0->GetAxis(3)->GetFirst();
424 Int_t last3 =his0->GetAxis(3)->GetLast();
426 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - 1/pt
427 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-fNCombineBin,1),TMath::Min(ibin3+fNCombineBin,nbins3));
428 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
429 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
431 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
432 Int_t first2 = his3->GetAxis(2)->GetFirst();
433 Int_t last2 = his3->GetAxis(2)->GetLast();
435 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
436 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-fNCombineBin,1),TMath::Min(ibin2+fNCombineBin,nbins2));
437 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
438 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
439 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
440 Int_t first1 = his2->GetAxis(1)->GetFirst();
441 Int_t last1 = his2->GetAxis(1)->GetLast();
442 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - tan(theta)
444 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
445 Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);
447 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1+fNCombineBin,nbins1));
448 if (TMath::Abs(x1)<(fNCombineBin+0.5)*binWidth){
449 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1,nbins1));
450 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+fNCombineBin,nbins1));
452 TH1 * hisDelta = his2->Projection(0);
454 Double_t entries = hisDelta->GetEntries();
455 Double_t mean=0, rms=0;
456 if (entries>kMinEntries){
457 mean = hisDelta->GetMean();
458 rms = hisDelta->GetRMS();
460 Double_t sector = 9.*x2/TMath::Pi();
461 if (sector<0) sector+=18;
462 Double_t dsec = sector-Int_t(sector)-0.5;
463 Double_t quantiles[5]={0};
464 Double_t mean75=0, rms75=0;
465 Double_t prob[5]={0, 0.25,0.5,0.75,1};
466 if (entries>kMinEntries){
467 hisDelta->GetQuantiles(5,quantiles,prob);
468 if(x3>0)hisDelta->GetXaxis()->SetRangeUser(quantiles[0], quantiles[3]);
469 else hisDelta->GetXaxis()->SetRangeUser(quantiles[1], quantiles[4]);
470 rms75=hisDelta->GetRMS();
471 mean75=hisDelta->GetMean();
474 Double_t pt = TMath::Abs(x3);
476 (*pcstream)<<"distortion"<<
479 "theta="<<x1<< // tan(theta)
482 "mpt="<<x3<< // signed 1/pt
484 "entries="<<entries<< // entries
485 "mean="<<mean<< // normal mean
487 "q25="<<quantiles[1]<< // quantiles of distribution - importnat for asymetric distibutions
488 "q50="<<quantiles[2]<< // quantiles of distribution - importnat for asymetric distibutions
489 "q75="<<quantiles[3]<< // quantiles of distribution - importnat for asymetric distibutions
490 "mean75="<<mean75<< // mean of the truncated distribution 0-75%
491 "rms75="<<rms75<< // rms of the truncated distibution up to 0-75%
493 "rms="<<rms<< // normal rms
494 "refX="<<refX<< // track matching refernce plane
495 "type="<<type<< // tpye of residuals
496 "sector="<<sector<< // sector number according to TPC standard
500 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
507 printf(" Finished! Close debug streamer \n");
513 //TFile f("outFile.root");
514 //TList * list = f.Get("jthaeder_OutterTracking");
515 TFile f("/u/alice-st/marr/EnergyCorrection/TrainAnalysis/trunk/marr_TrackingTest.root");
516 TList * list = f.Get("marr_TrackingTest");
518 AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
519 AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
520 AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
521 AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
524 compEMCAL-> MakeDistortionMap(refX,0);
525 compEMCAL-> MakeDistortionMap(refX,1);
526 compEMCAL-> MakeDistortionMap(refX,4);
527 compEMCALEl-> MakeDistortionMap(refX,0);
528 compEMCALEl-> MakeDistortionMap(refX,1);
529 compEMCALEl-> MakeDistortionMap(refX,4);
531 compTOF->SetNCombineBin(3)
532 compTOF->MakeDistortionMap(365,4);
533 compTOF->MakeDistortionMap(365,0);
535 TFile f("emcalDistortion.root");
538 // ideal case - no mislaignment dy only due energy loss correction
539 TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
541 TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
543 TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
545 TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10"); // pions - bias +-1 %
547 TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10"); // electrons - bias +-20 %