AliPhysics  5364b50 (5364b50)
AliTrackComparison.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
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  **************************************************************************/
15 
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
21 //
22 // Anthor: R.Ma, M.Ivanov 02/04/2011
23 //--------------------------------------------------------------------
24 /* $Id:$ */
25 
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"
33 #include "TMatrix.h"
34 #include "TParticlePDG.h"
35 #include "TParticle.h"
36 #include "TTreeStream.h"
37 #include "TAxis.h"
38 #include "TH1D.h"
39 #include "TTreeStream.h"
40 
41 ClassImp(AliTrackComparison)
42 
43 //________________________________________________________________________
45  :TNamed()
46  , fStep(1)
47  , fLowBinDY(-10)
48  , fUpBinDY(10)
49  , fLowBinDZ(-10)
50  , fUpBinDZ(10)
51  , fLowBinDSnp(-0.5)
52  , fUpBinDSnp(0.5)
53  , fLowBinDTheta(-0.5)
54  , fUpBinDTheta(0.5)
55  , fLowBin1Pt(-3)
56  , fUpBin1Pt(3)
57  , fLowBin1PtLoss(-2)
58  , fUpBin1PtLoss(2)
59  , fNBinsDY(50)
60  , fNBinsDZ(50)
61  , fNBinsDSnp(50)
62  , fNBinsDTheta(50)
63  , fNBins1Pt(50)
64  , fNBins1PtLoss(100)
65  , fLayerID(-1)
66  , fFillAll(kTRUE)
67  , fNCombineBin(1)
68 {
69  //
70  // Default constructor
71  //
72  for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
73 }
74 
75 
76 //________________________________________________________________________
77 AliTrackComparison::AliTrackComparison(const Text_t *name, const Text_t *title)
78  :TNamed(name,title)
79  , fStep(1)
80  , fLowBinDY(-10)
81  , fUpBinDY(10)
82  , fLowBinDZ(-10)
83  , fUpBinDZ(10)
84  , fLowBinDSnp(-0.5)
85  , fUpBinDSnp(0.5)
86  , fLowBinDTheta(-0.5)
87  , fUpBinDTheta(0.5)
88  , fLowBin1Pt(-3)
89  , fUpBin1Pt(3)
90  , fLowBin1PtLoss(-2)
91  , fUpBin1PtLoss(2)
92  , fNBinsDY(50)
93  , fNBinsDZ(50)
94  , fNBinsDSnp(50)
95  , fNBinsDTheta(50)
96  , fNBins1Pt(50)
97  , fNBins1PtLoss(100)
98  , fLayerID(-1)
99  , fFillAll(kTRUE)
100  , fNCombineBin(1)
101 {
102  //
103  // Non default cosntructor
104  //
105  for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
106 }
107 
108 //________________________________________________________________________
110  :TNamed(comp)
111  , fStep(comp.fStep)
112  , fLowBinDY(comp.fLowBinDY)
113  , fUpBinDY(comp.fUpBinDY)
114  , fLowBinDZ(comp.fLowBinDZ)
115  , fUpBinDZ(comp.fUpBinDZ)
116  , fLowBinDSnp(comp.fLowBinDSnp)
117  , fUpBinDSnp(comp.fUpBinDSnp)
119  , fUpBinDTheta(comp.fUpBinDTheta)
120  , fLowBin1Pt(comp.fLowBin1Pt)
121  , fUpBin1Pt(comp.fUpBin1Pt)
124  , fNBinsDY(comp.fNBinsDY)
125  , fNBinsDZ(comp.fNBinsDZ)
126  , fNBinsDSnp(comp.fNBinsDSnp)
127  , fNBinsDTheta(comp.fNBinsDTheta)
128  , fNBins1Pt(comp.fNBins1Pt)
130  , fLayerID(comp.fLayerID)
131  , fFillAll(comp.fFillAll)
132  , fNCombineBin(comp.fNCombineBin)
133 {
134  //
135  // copy constructor
136  //
137 
138  for (Int_t i=0;i<6; i++) fResolHisto[i]=comp.fResolHisto[i];
139 }
140 
141 //________________________________________________________________________
143 {
144  //
145  //
146  //
147  if(this != &comp) {
148  TNamed::operator=(comp);
149  }
150  return *this;
151 }
152 
153 //________________________________________________________________________
155  //
156  //Initilized the output histograms
157  //
158  MakeHistos();
159 
160 }
161 
162 //________________________________________________________________________
164  //
165  //
166  //
167 }
168 
169 //________________________________________________________________________
171  //
172  //
173  //
174 }
175 
176 //________________________________________________________________________
177 Int_t AliTrackComparison::AddTracks(AliTrackReference *ref0, AliTrackReference *ref1, Double_t mass, Int_t charge)
178 {
179  // Make track param out of track reference
180  // Test propagation from ref0 to ref1
181  // Fill the THnSparse with ref0 and ref1
182 
183  AliExternalTrackParam *param0 = 0;
184  AliExternalTrackParam *param1 = 0;
185 
186  param0=MakeTrack(ref0,charge);
187  param1=MakeTrack(ref1,charge);
188 
189  if (!param0 || !param1) return 0;
190 
191  Double_t tr1Pt = param0->GetSigned1Pt();
192 
193  if(!PropagateToPoint(param0,param1,mass)) return 0;
194 
195  FillHistos(param0,param1,tr1Pt);
196  return 1;
197 }
198 
199 //________________________________________________________________________
200 Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
201 {
202  //Test propagation from param0 to param1
203  //Fill the THnSparse with param0 and param1
204  //
205  Double_t tr1Pt=param0->GetSigned1Pt();
206  if( !PropagateToPoint(param0,param1,mass)) return 0;
207  FillHistos(param0,param1,tr1Pt);
208  return 1;
209 }
210 
211 //________________________________________________________________________
212 Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0, AliTrackPoint *point1, Double_t mass, Double_t energy, Double_t *vxyz)
213 {
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)
220 
221  Double_t gPos[3]= {point1->GetX(),point1->GetY(),point1->GetZ()};
222 
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;
228 
229  Double_t cv[21];
230  for (Int_t i=0; i<21;i++) cv[i]=0;
231  AliExternalTrackParam * param1 = new AliExternalTrackParam(gPos,pxyz,cv,param0->Charge());
232 
233  if(!param1) return 0;
234  Double_t tr1Pt = param0->GetSigned1Pt();
235  if(!PropagateToPoint(param0,param1,mass)) return 0;
236 
237  FillHistos(param0,param1,tr1Pt);
238  return 1;
239 }
240 
241 //________________________________________________________________________
242 Bool_t AliTrackComparison::PropagateToPoint(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
243 {
244  //
245  //Extrapolate is performed
246  //
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());
251  return isOK;
252 }
253 
254 //________________________________________________________________________
255 AliExternalTrackParam * AliTrackComparison::MakeTrack(const AliTrackReference* ref, Int_t charge)
256 {
257  //
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()};
262  Double_t cv[21];
263  for (Int_t i=0; i<21;i++) cv[i]=0;
264  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
265  return param;
266 }
267 
268 //________________________________________________________________________
270  //
271  //Merge the comparison objects
272  //
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())) {
278  return -1;
279  }
280  if (strName.CompareTo(comp->GetName())!=0) return -1;
281  // add histograms here...
282  Add(comp);
283  }
284  return 0;
285 }
286 
287 //________________________________________________________________________
289  //
290  // Add THnSparse
291  //
292  for (Int_t i=0;i<6;i++){
293  if (!fResolHisto[i])
294  continue;
295  THnSparse * h0 = (THnSparse*)fResolHisto[i];
296  THnSparse * h1 = (THnSparse*)comp->GetHnSparse(i);
297  if (h0&&h1) h0->Add(h1);
298  }
299 }
300 
301 
302 //________________________________________________________________________
304 {
305  //
306  //Called in Init() to initialize histograms
307  //
308  Double_t xminTrack[5], xmaxTrack[5];
309  Int_t binsTrack[5];
310  TString axisName[5];
311  TString axisTitle[5];
312  TString hisNames[6]={"DeltaY","DeltaZ","DeltaSnp","DeltaTheta","Delta1Pt","1PtLoss"};
316  //
317  axisName[0] ="#Delta";
318  axisTitle[0] ="#Delta";
319  //
320  binsTrack[1] =40;
321  xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
322  axisName[1] ="tanTheta";
323  axisTitle[1] ="tan(#Theta)";
324  //
325  binsTrack[2] =180;
326  xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
327  axisName[2] ="phi";
328  axisTitle[2] ="#phi";
329  //
330  binsTrack[3] =120;
331  xminTrack[3] =-3.; xmaxTrack[3]=3.; // 0.33 GeV cut
332  axisName[3] ="1pt";
333  axisTitle[3] ="1/pt";
334  //
335  binsTrack[4] =22;
336  xminTrack[4] =0; xmaxTrack[4]=22;
337  axisName[4] ="LayerID";
338  axisTitle[4] ="LayerID";
339 
340  for (Int_t ihis=0; ihis<6; ihis++){
341  // modify ranges for bin 0
342  //
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());
350  }
351  }
352 }
353 
354 //________________________________________________________________________
355 void AliTrackComparison::FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
356 {
357  //Fill the THnSparse.
358  //In case of not filling all, only the 4 out of 6 THnSparse are filed
359  //
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
366 
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());
371 
372  Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
373  Double_t argu[6][5];
374  for(Int_t j=0; j<6; j++)
375  {
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]);
383  }
384 }
385 
386 //________________________________________________________________________
388  //
389  // make a distortion map out of the residual histogram
390  // Results are written to the debug streamer - pcstream
391  // Parameters:
392  // pcstream - file to write the tree
393  // run - run number
394  // refX - track matching reference X
395  // type - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
396  // THnSparse axes:
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
403 
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);
408 
409  TString dsName;
410  dsName=GetName();
411  dsName+=Form("Debug_%d.root",type);
412  printf(" Create debug streamer \n");
413  dsName.ReplaceAll(" ","");
414  TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());
415 
416 
417  const Int_t kMinEntries=10;
418  Double_t bz=AliTrackerBase::GetBz();
419  Int_t idim[4]={0,1,2,3};
420  //
421  //
422  //
423  Int_t nbins3=his0->GetAxis(3)->GetNbins();
424  Int_t first3=his0->GetAxis(3)->GetFirst();
425  Int_t last3 =his0->GetAxis(3)->GetLast();
426  //
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
431  //
432  Int_t nbins2 = his3->GetAxis(2)->GetNbins();
433  Int_t first2 = his3->GetAxis(2)->GetFirst();
434  Int_t last2 = his3->GetAxis(2)->GetLast();
435  //
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)
444  //
445  Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
446  Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);
447 
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));
452  }
453  TH1 * hisDelta = his2->Projection(0);
454  //
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();
460  }
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();
473  }
474 
475  Double_t pt = TMath::Abs(x3);
476  Double_t z=refX*x1;
477  (*pcstream)<<"distortion"<<
478  "run="<<run<<
479  "bz="<<bz<<
480  "theta="<<x1<< // tan(theta)
481  "phi="<<x2<< // phi
482  "z="<<z<< // dummy z
483  "mpt="<<x3<< // signed 1/pt
484  "1Pt="<<pt<<
485  "entries="<<entries<< // entries
486  "mean="<<mean<< // normal mean
487  //
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%
493  //
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
498  "dsec="<<dsec<<
499  "\n";
500  delete hisDelta;
501  //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
502  }
503  delete his2;
504  }
505  delete his3;
506  }
507 
508  printf(" Finished! Close debug streamer \n");
509  delete pcstream;
510 }
511 
512 /*
513  .x ~/rootlogon.C
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");
518 
519  AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
520  AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
521  AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
522  AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
523  Double_t refX=438;
524 
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);
531 
532  compTOF->SetNCombineBin(3)
533  compTOF->MakeDistortionMap(365,4);
534  compTOF->MakeDistortionMap(365,0);
535 
536  TFile f("emcalDistortion.root");
537 
538 
539  // ideal case - no mislaignment dy only due energy loss correction
540  TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
541  //
542  TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
543  // delta 1/pt
544  TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
545  // (delta pt)/pt
546  TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10"); // pions - bias +-1 %
547  //
548  TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10"); // electrons - bias +-20 %
549  //
550 
551 
552 
553 
554 */
Int_t charge
THnSparse * GetHnSparse(Int_t index)
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
virtual Long64_t Merge(TCollection *const li)
long long Long64_t
Definition: External.C:43
Double_t mass
energy
Definition: HFPtSpectrum.C:44
Double_t bz
Int_t AddTracks(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
int Int_t
Definition: External.C:63
AliTrackComparison & operator=(const AliTrackComparison &comp)
AliExternalTrackParam * MakeTrack(const AliTrackReference *ref, Int_t charge)
Bool_t PropagateToPoint(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
void FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
void MakeDistortionMap(Double_t refX, Int_t type)
virtual void Add(AliTrackComparison *const comp)
bool Bool_t
Definition: External.C:53
THnSparse * fResolHisto[6]
Definition: External.C:196