Add INT7 for pA.# Please enter the commit message for your changes. s starting
[u/mrichter/AliRoot.git] / PWGPP / AliTrackComparison.cxx
CommitLineData
bdf94dc3 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
41ClassImp(AliTrackComparison)
42
43//________________________________________________________________________
44AliTrackComparison::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//________________________________________________________________________
77AliTrackComparison::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//________________________________________________________________________
109AliTrackComparison::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//________________________________________________________________________
142AliTrackComparison& AliTrackComparison::operator=(const AliTrackComparison& comp)
143{
144 //
145 //
146 //
147 if(this != &comp) {
148 TNamed::operator=(comp);
149 }
150 return *this;
151}
152
153//________________________________________________________________________
154void AliTrackComparison::Init(){
155 //
156 //Initilized the output histograms
157 //
158 MakeHistos();
159
160}
161
162//________________________________________________________________________
163AliTrackComparison::~AliTrackComparison(){
164 //
165 //
166 //
167}
168
169//________________________________________________________________________
170void AliTrackComparison::Analyze() {
171 //
172 //
173 //
174}
175
176//________________________________________________________________________
177Int_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//________________________________________________________________________
200Int_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//________________________________________________________________________
212Int_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];
f6139574 226 Double_t R = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
bdf94dc3 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//________________________________________________________________________
242Bool_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//________________________________________________________________________
255AliExternalTrackParam * 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//________________________________________________________________________
269Long64_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//________________________________________________________________________
288void AliTrackComparison::Add(AliTrackComparison *const comp){
289 //
290 // Add THnSparse
291 //
bdf94dc3 292 for (Int_t i=0;i<6;i++){
f6139574 293 if (!fResolHisto[i])
294 continue;
bdf94dc3 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//________________________________________________________________________
303void AliTrackComparison::MakeHistos()
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"};
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};
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//________________________________________________________________________
355void 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//________________________________________________________________________
387void AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
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*/