remove obsolate classes
[u/mrichter/AliRoot.git] / PWG1 / AliTrackComparison.cxx
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 //________________________________________________________________________
44 AliTrackComparison::AliTrackComparison() 
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 //________________________________________________________________________
109 AliTrackComparison::AliTrackComparison(const AliTrackComparison& comp)
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)
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)
133 {
134   //
135   // copy constructor
136   //
137
138   for (Int_t i=0;i<6; i++) fResolHisto[i]=comp.fResolHisto[i];
139 }
140
141 //________________________________________________________________________
142 AliTrackComparison& AliTrackComparison::operator=(const AliTrackComparison& comp)
143 {
144   //
145   //
146   //
147   if(this != &comp) {
148     TNamed::operator=(comp);
149   }
150   return *this;
151 }
152
153 //________________________________________________________________________
154 void AliTrackComparison::Init(){
155   //
156   //Initilized the output histograms
157   //
158   MakeHistos();
159
160 }
161
162 //________________________________________________________________________
163 AliTrackComparison::~AliTrackComparison(){
164   //
165   //
166   //
167 }
168  
169 //________________________________________________________________________
170 void AliTrackComparison::Analyze() {
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[3]);
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 //________________________________________________________________________
269 Long64_t AliTrackComparison::Merge(TCollection *const li) {
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 //________________________________________________________________________
288 void  AliTrackComparison::Add(AliTrackComparison *const comp){
289   //
290   // Add THnSparse
291   //
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);
297   }
298 }
299
300
301 //________________________________________________________________________
302 void AliTrackComparison::MakeHistos()
303 {
304   //
305   //Called in Init() to initialize histograms
306   // 
307   Double_t xminTrack[5], xmaxTrack[5];
308   Int_t   binsTrack[5];
309   TString axisName[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};
315   //
316   axisName[0]   ="#Delta";
317   axisTitle[0]  ="#Delta";
318   //
319   binsTrack[1] =40;
320   xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
321   axisName[1]  ="tanTheta";
322   axisTitle[1]  ="tan(#Theta)";
323   //
324   binsTrack[2] =180;
325   xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi(); 
326   axisName[2]  ="phi";
327   axisTitle[2]  ="#phi";
328   //
329   binsTrack[3] =120;
330   xminTrack[3] =-3.; xmaxTrack[3]=3.;   // 0.33 GeV cut 
331   axisName[3]  ="1pt";
332   axisTitle[3]  ="1/pt";
333   //
334   binsTrack[4] =22;
335   xminTrack[4] =0; xmaxTrack[4]=22;  
336   axisName[4]  ="LayerID";
337   axisTitle[4] ="LayerID";
338
339   for (Int_t ihis=0; ihis<6; ihis++){
340     // modify ranges for bin 0
341     //
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());
349     }
350   }
351 }
352
353 //________________________________________________________________________
354 void AliTrackComparison::FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
355 {
356   //Fill the THnSparse. 
357   //In case of not filling all, only the 4 out of 6 THnSparse are filed
358   //
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
365
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());
370
371   Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
372   Double_t argu[6][5];
373   for(Int_t j=0; j<6; j++)
374     {
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]);
382     } 
383 }
384
385 //________________________________________________________________________
386 void   AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
387   //
388   // make a distortion map out of the residual histogram
389   // Results are written to the debug streamer - pcstream
390   // Parameters:
391   //   pcstream   - file to write the tree
392   //   run        - run number
393   //   refX       - track matching reference X
394   //   type       - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
395   // THnSparse axes:
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
402   
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);
407   
408   TString dsName;
409   dsName=GetName();
410   dsName+=Form("Debug_%d.root",type);
411   printf(" Create debug streamer \n");
412   dsName.ReplaceAll(" ","");
413   TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());
414
415
416   const Int_t kMinEntries=10;
417   Double_t bz=AliTrackerBase::GetBz();
418   Int_t idim[4]={0,1,2,3};
419   //
420   //
421   //
422   Int_t nbins3=his0->GetAxis(3)->GetNbins();
423   Int_t first3=his0->GetAxis(3)->GetFirst();
424   Int_t last3 =his0->GetAxis(3)->GetLast();
425   //
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
430     //
431     Int_t nbins2    = his3->GetAxis(2)->GetNbins();
432     Int_t first2    = his3->GetAxis(2)->GetFirst();
433     Int_t last2     = his3->GetAxis(2)->GetLast();
434     //
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)
443         //
444         Double_t       x1= his2->GetAxis(1)->GetBinCenter(ibin1);
445         Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);
446
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));
451         }
452         TH1 * hisDelta = his2->Projection(0);
453         //
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(); 
459         }
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();
472         }
473         
474         Double_t pt = TMath::Abs(x3);
475         Double_t z=refX*x1;
476         (*pcstream)<<"distortion"<<
477           "run="<<run<<
478           "bz="<<bz<<
479           "theta="<<x1<<          // tan(theta)
480           "phi="<<x2<<            // phi
481           "z="<<z<<               // dummy z
482           "mpt="<<x3<<            // signed 1/pt
483           "1Pt="<<pt<<
484           "entries="<<entries<<   // entries
485           "mean="<<mean<<         // normal mean
486           //
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%
492           //
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
497           "dsec="<<dsec<<
498           "\n";
499         delete hisDelta;
500         //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
501       }
502       delete his2;
503     }
504     delete his3;
505   }
506
507   printf(" Finished! Close debug streamer \n");
508   delete pcstream;
509 }
510
511 /*
512   .x ~/rootlogon.C
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");
517
518   AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
519   AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
520   AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
521   AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
522   Double_t refX=438;
523
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);
530   
531   compTOF->SetNCombineBin(3)
532   compTOF->MakeDistortionMap(365,4);
533   compTOF->MakeDistortionMap(365,0);
534   
535   TFile f("emcalDistortion.root");
536   
537
538   // ideal case - no mislaignment dy only due energy loss correction
539   TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
540   //
541   TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
542   // delta 1/pt
543   TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
544   // (delta pt)/pt
545   TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10");  // pions  - bias +-1 %
546   //
547   TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");  // electrons  - bias +-20 %
548   //
549
550
551
552
553 */