]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliTrackComparison.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
first implementation of the global alignment QA
[u/mrichter/AliRoot.git] / PWG1 / 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];
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//________________________________________________________________________
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 //
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//________________________________________________________________________
302void 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//________________________________________________________________________
354void 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//________________________________________________________________________
386void 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*/