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 **************************************************************************/
19 Responsible: marian.ivanov@cern.ch
20 Tools for fitting of the space point distortion parameters.
24 1. Creation of the distortion maps from the residual histograms
26 3. Parameters to calculate- dO/dp
28 Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
29 extracted from the old macros
34 #include "Riostream.h"
37 #include "TGraphErrors.h"
38 #include "AliExternalTrackParam.h"
41 #include "TMultiGraph.h"
43 #include "THnSparse.h"
51 #include "AliTPCROC.h"
52 #include "AliTPCCalROC.h"
53 #include "AliESDfriend.h"
54 #include "AliTPCcalibTime.h"
55 #include "AliSplineFit.h"
56 #include "AliCDBMetaData.h"
58 #include "AliCDBManager.h"
59 #include "AliCDBStorage.h"
60 #include "AliTPCcalibBase.h"
61 #include "AliTPCcalibDB.h"
62 #include "AliTPCcalibDButil.h"
63 #include "AliRelAlignerKalman.h"
64 #include "AliTPCParamSR.h"
65 #include "AliTPCcalibTimeGain.h"
66 #include "AliTPCcalibGainMult.h"
67 #include "AliSplineFit.h"
68 #include "AliTPCComposedCorrection.h"
69 #include "AliTPCExBTwist.h"
70 #include "AliTPCCalibGlobalMisalignment.h"
71 #include "TStatToolkit.h"
74 #include "AliTrackerBase.h"
75 #include "AliTPCCorrectionFit.h"
76 #include "AliTPCLaserTrack.h"
77 #include "TDatabasePDG.h"
78 #include "AliTPCcalibAlign.h"
81 ClassImp(AliTPCCorrectionFit)
83 AliTPCCorrectionFit::AliTPCCorrectionFit():
84 TNamed("TPCCorrectionFit","TPCCorrectionFit")
87 // default constructor
91 AliTPCCorrectionFit::~AliTPCCorrectionFit() {
98 Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){
102 Double_t sector = 9*phi/TMath::Pi();
103 if (sector<0) sector+=18;
104 Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
105 Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
106 if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.);
107 if (ptype==2) return (y245-y85)/(245.-85.);
114 Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){
116 // Fit the distortion along the line with the parabolic model
118 // phi0 - phi at the entrance of the TPC
119 // snp - local inclination angle at the entrance of the TPC
120 // refX - ref X where the distortion is evanluated
123 static TLinearFitter fitter(3,"pol2");
124 fitter.ClearPoints();
125 if (nsteps<3) nsteps=3;
126 Double_t deltaX=(245-85)/(nsteps);
127 for (Int_t istep=0; istep<(nsteps+1); istep++){
129 Double_t localX =85.+deltaX*istep;
130 Double_t localPhi=phi0+deltaX*snp*istep;
131 Double_t sector = 9*localPhi/TMath::Pi();
132 if (sector<0) sector+=18;
133 Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
134 Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
135 Double_t x[1]={localX-dlocalX};
136 fitter.AddPoint(x,y);
140 par[0]=fitter.GetParameter(0);
141 par[1]=fitter.GetParameter(1);
142 par[2]=fitter.GetParameter(2);
144 if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX;
145 if (ptype==2) return par[1]+2*par[2]*refX;
146 if (ptype==4) return par[2];
153 void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
155 // Create cluster distortion map
157 TFile *falign = TFile::Open("CalibObjects.root");
158 TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
160 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
163 AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
165 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
168 TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
170 THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
171 //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
172 AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
173 // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
175 const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
176 for (Int_t ihis=0; ihis<4; ihis++){
177 THnSparse * hisAlign =align->GetTrackletDelta(ihis);
178 AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
188 void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
190 // Make cluster residual map from the n-dimensional histogram
191 // hisInput supposed to have given format:
197 // Vertex position assumed to be at (0,0,0)
199 //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
201 Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
202 Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
203 Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
206 new TH3F("his3D","his3D",
207 nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
208 nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
209 nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
210 hisResMap3D->GetXaxis()->SetTitle("sector");
211 hisResMap3D->GetYaxis()->SetTitle("localX");
212 hisResMap3D->GetZaxis()->SetTitle("kZ");
214 TH2F * hisResMap2D[4] ={0,0,0,0};
215 for (Int_t i=0; i<4; i++){
217 new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
218 nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
219 nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
220 hisResMap2D[i]->GetXaxis()->SetTitle("sector");
221 hisResMap2D[i]->GetYaxis()->SetTitle("localX");
227 Int_t axis0[4]={0,1,2,3};
228 Int_t axis1[4]={0,1,2,3};
230 for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
232 hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
233 THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
234 Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
236 for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
239 his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
240 THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
241 Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
244 his2->GetAxis(3)->SetRangeUser(0.01,0.3);
245 TH1 * hisA = his2->Projection(0);
246 Double_t meanA= hisA->GetMean();
247 Double_t rmsA= hisA->GetRMS();
248 Double_t entriesA= hisA->GetEntries();
251 his2->GetAxis(3)->SetRangeUser(0.01,0.3);
252 TH1 * hisC = his2->Projection(0);
253 Double_t meanC= hisC->GetMean();
254 Double_t rmsC= hisC->GetRMS();
255 Double_t entriesC= hisC->GetEntries();
257 his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
258 TH2 * hisAC = his2->Projection(0,3);
259 TProfile *profAC = hisAC->ProfileX();
261 profAC->Fit("pol1","QNR","QNR",0.05,1);
262 if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
263 Double_t offsetA=f1->GetParameter(0);
264 Double_t slopeA=f1->GetParameter(1);
265 Double_t offsetAE=f1->GetParError(0);
266 Double_t slopeAE=f1->GetParError(1);
267 Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
268 profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
269 f1=(TF1*)gROOT->FindObject("pol1");
270 Double_t offsetC=f1->GetParameter(0);
271 Double_t slopeC=f1->GetParameter(1);
272 Double_t offsetCE=f1->GetParError(0);
273 Double_t slopeCE=f1->GetParError(1);
274 Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
275 if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C);
277 (*pcstream)<<"deltaFit"<<
282 "entriesA="<<entriesA<<
285 "entriesC="<<entriesC<<
286 "offsetA="<<offsetA<<
288 "offsetAE="<<offsetAE<<
289 "slopeAE="<<slopeAE<<
291 "offsetC="<<offsetC<<
293 "offsetCE="<<offsetCE<<
294 "slopeCE="<<slopeCE<<
298 hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
299 hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
300 hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
301 hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
303 for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
304 Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
305 if (TMath::Abs(kZ)<0.05) continue; // crossing
306 his2->GetAxis(3)->SetRange(ibin3,ibin3);
307 if (TMath::Abs(kZ)>0.15){
308 his2->GetAxis(3)->SetRange(ibin3,ibin3);
310 TH1 * his = his2->Projection(0);
311 Double_t mean= his->GetMean();
312 Double_t rms= his->GetRMS();
313 Double_t entries= his->GetEntries();
314 //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
315 hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
316 Double_t phi=TMath::Pi()*sector/9;
317 if (phi>TMath::Pi()) phi+=TMath::Pi();
322 his->Fit("gaus","Q","goff");
323 fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
326 his->Fit(fgaus,"Q","goff");
327 meanG=fgaus->GetParameter(1);
328 rmsG=fgaus->GetParameter(2);
331 Double_t dsec=sector-Int_t(sector)-0.5;
332 Double_t snp=dsec*TMath::Pi()/9.;
333 (*pcstream)<<"delta"<<
347 "entries="<<entries<<
350 "entriesA="<<entriesA<<
353 "entriesC="<<entriesC<<
354 "offsetA="<<offsetA<<
357 "offsetC="<<offsetC<<
367 hisResMap3D->Write();
368 hisResMap2D[0]->Write();
369 hisResMap2D[1]->Write();
370 hisResMap2D[2]->Write();
371 hisResMap2D[3]->Write();
377 void AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){
379 // make a distortion map out ou fthe residual histogram
380 // Results are written to the debug streamer - pcstream
382 // his0 - input (4D) residual histogram
383 // pcstream - file to write the tree
385 // refX - track matching reference X
386 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
388 // OBJ: TAxis #Delta #Delta
389 // OBJ: TAxis tanTheta tan(#Theta)
390 // OBJ: TAxis phi #phi
391 // OBJ: TAxis snp snp
393 // marian.ivanov@cern.ch
394 const Int_t kMinEntries=10;
395 Int_t idim[4]={0,1,2,3};
399 Int_t nbins3=his0->GetAxis(3)->GetNbins();
400 Int_t first3=his0->GetAxis(3)->GetFirst();
401 Int_t last3 =his0->GetAxis(3)->GetLast();
403 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
404 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
405 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
406 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
408 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
409 Int_t first2 = his3->GetAxis(2)->GetFirst();
410 Int_t last2 = his3->GetAxis(2)->GetLast();
412 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
413 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
414 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
415 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
416 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
417 Int_t first1 = his2->GetAxis(1)->GetFirst();
418 Int_t last1 = his2->GetAxis(1)->GetLast();
419 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
421 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
422 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
423 if (TMath::Abs(x1)<0.1){
424 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
425 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
427 if (TMath::Abs(x1)<0.06){
428 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
430 TH1 * hisDelta = his2->Projection(0);
432 Double_t entries = hisDelta->GetEntries();
433 Double_t mean=0, rms=0;
434 if (entries>kMinEntries){
435 mean = hisDelta->GetMean();
436 rms = hisDelta->GetRMS();
438 Double_t sector = 9.*x2/TMath::Pi();
439 if (sector<0) sector+=18;
440 Double_t dsec = sector-Int_t(sector)-0.5;
449 "entries="<<entries<<
452 "refX="<<refX<< // track matching refernce plane
458 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
468 void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
470 // make a distortion map out ou fthe residual histogram
471 // Results are written to the debug streamer - pcstream
473 // his0 - input (4D) residual histogram
474 // pcstream - file to write the tree
476 // refX - track matching reference X
477 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
478 // marian.ivanov@cern.ch
481 // Collection name='TObjArray', class='TObjArray', size=16
482 // 0. OBJ: TAxis #Delta #Delta
483 // 1. OBJ: TAxis N_{cl} N_{cl}
484 // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
485 // 3. OBJ: TAxis z (cm) z (cm)
486 // 4. OBJ: TAxis sin(#phi) sin(#phi)
487 // 5. OBJ: TAxis tan(#theta) tan(#theta)
488 // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
489 // 7. OBJ: TAxis pt (GeV) pt (GeV)
490 // 8. OBJ: TAxis alpha alpha
491 const Int_t kMinEntries=10;
493 // 1. make default selections
496 Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
497 hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
498 hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
499 hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
500 hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
501 hisDelta= hisInput->Projection(0);
502 hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
504 THnSparse *his0= hisInput->Projection(4,idim0);
506 // 2. Get mean in diferent bins
508 Int_t nbins1=his0->GetAxis(1)->GetNbins();
509 Int_t first1=his0->GetAxis(1)->GetFirst();
510 Int_t last1 =his0->GetAxis(1)->GetLast();
512 Double_t bz=AliTrackerBase::GetBz();
513 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
515 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
517 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
518 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
520 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
521 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
522 Int_t first3 = his1->GetAxis(3)->GetFirst();
523 Int_t last3 = his1->GetAxis(3)->GetLast();
525 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
526 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
527 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
529 his1->GetAxis(3)->SetRangeUser(-1,1);
532 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
533 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
534 Int_t first2 = his3->GetAxis(2)->GetFirst();
535 Int_t last2 = his3->GetAxis(2)->GetLast();
537 for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
538 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
539 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
540 hisDelta = his3->Projection(0);
542 Double_t entries = hisDelta->GetEntries();
543 Double_t mean=0, rms=0;
544 if (entries>kMinEntries){
545 mean = hisDelta->GetMean();
546 rms = hisDelta->GetRMS();
548 Double_t sector = 9.*x2/TMath::Pi();
549 if (sector<0) sector+=18;
550 Double_t dsec = sector-Int_t(sector)-0.5;
551 Double_t snp=0; // dummy snp - equal 0
554 "bz="<<bz<< // magnetic field
555 "theta="<<x1<< // theta
556 "phi="<<x2<< // phi (alpha)
557 "z="<<x3<< // z at "vertex"
558 "snp="<<snp<< // dummy snp
559 "entries="<<entries<< // entries in bin
560 "mean="<<mean<< // mean
562 "refX="<<refX<< // track matching refernce plane
563 "type="<<type<< // parameter type
564 "sector="<<sector<< // sector
565 "dsec="<<dsec<< // dummy delta sector
568 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
579 void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
581 // make a distortion map out of the residual histogram
582 // Results are written to the debug streamer - pcstream
584 // his0 - input (4D) residual histogram
585 // pcstream - file to write the tree
587 // type - 0- y 1-z,2 -snp, 3-theta
588 // bz - magnetic field
589 // marian.ivanov@cern.ch
591 //Collection name='TObjArray', class='TObjArray', size=16
592 //0 OBJ: TAxis delta delta
593 //1 OBJ: TAxis phi phi
594 //2 OBJ: TAxis localX localX
597 //5 OBJ: TAxis is1 is1
598 //6 OBJ: TAxis is0 is0
600 //8. OBJ: TAxis IsPrimary IsPrimary
602 const Int_t kMinEntries=10;
603 THnSparse * hisSector0=0;
604 TH1 * htemp=0; // histogram to calculate mean value of parameter
605 // Double_t bz=AliTrackerBase::GetBz();
608 // Loop over pair of sector:
619 for (Int_t isec0=0; isec0<72; isec0++){
620 Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
622 //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
623 hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
624 hisSector0=hisInput->Projection(7,index0);
627 for (Int_t isec1=isec0+1; isec1<72; isec1++){
628 //if (isec1!=isec0+36) continue;
629 if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
630 printf("Sectors %d\t%d\n",isec1,isec0);
631 hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
632 TH1 * hisX=hisSector0->Projection(5);
633 Double_t refX= hisX->GetMean();
635 TH1 *hisDelta=hisSector0->Projection(0);
636 Double_t dmean = hisDelta->GetMean();
637 Double_t drms = hisDelta->GetRMS();
638 hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
641 // 1. make default selections
643 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
644 THnBase *hisSector1= hisSector0->ProjectionND(5,idim0);
646 // 2. Get mean in diferent bins
648 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
650 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
651 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
652 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
654 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi
658 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
659 Double_t psec = (9*xPhi/TMath::Pi());
660 if (psec<0) psec+=18;
663 if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
664 if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
665 if (!isOK0) continue;
666 if (!isOK1) continue;
668 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
669 if (isec1!=isec0+36) {
670 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
673 htemp = hisSector1->Projection(4);
674 xPhi=htemp->GetMean();
676 THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
677 //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
678 Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
679 Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
681 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z
685 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
686 if (isec1!=isec0+36) {
687 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
689 htemp = hisPhi->Projection(3);
690 Double_t xZ= htemp->GetMean();
692 THnBase * hisZ= hisPhi->ProjectionND(3,idim);
693 //projected histogram according selection 3 -z
696 //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
697 Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
698 Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
699 for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
703 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
704 if (isec1!=isec0+36) {
705 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
707 htemp = hisZ->Projection(2);
708 Double_t xSnp= htemp->GetMean();
710 THnBase * hisSnp= hisZ->ProjectionND(2,idim);
711 //projected histogram according selection 2 - snp
713 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
714 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
715 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
717 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
720 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
721 if (isec1!=isec0+36) {
722 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
724 htemp = hisSnp->Projection(1);
725 Double_t xTheta=htemp->GetMean();
727 hisDelta = hisSnp->Projection(0);
729 Double_t entries = hisDelta->GetEntries();
730 Double_t mean=0, rms=0;
731 if (entries>kMinEntries){
732 mean = hisDelta->GetMean();
733 rms = hisDelta->GetRMS();
735 Double_t sector = 9.*xPhi/TMath::Pi();
736 if (sector<0) sector+=18;
737 Double_t dsec = sector-Int_t(sector)-0.5;
738 Int_t dtype=1; // TPC alignment type
741 "bz="<<bz<< // magnetic field
742 "ptype="<<type<< // parameter type
743 "dtype="<<dtype<< // parameter type
744 "isec0="<<isec0<< // sector 0
745 "isec1="<<isec1<< // sector 1
746 "sector="<<sector<< // sector as float
747 "dsec="<<dsec<< // delta sector
749 "theta="<<xTheta<< // theta
750 "phi="<<xPhi<< // phi (alpha)
752 "snp="<<xSnp<< // snp
755 "entries="<<entries<< // entries in bin
756 "mean="<<mean<< // mean
758 "refX="<<refX<< // track matching reference plane
761 //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean);
780 // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
782 // // Make a laser fit tree for global minimization
784 // AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
785 // AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
786 // if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
787 // correction->AddVisualCorrection(correction,0); //register correction
789 // // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
790 // //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
792 // const Double_t cutErrY=0.05;
793 // const Double_t kSigmaCut=4;
794 // // const Double_t cutErrZ=0.03;
795 // const Double_t kEpsilon=0.00000001;
796 // // const Double_t kMaxDist=1.; // max distance - space correction
797 // TVectorD *vecdY=0;
798 // TVectorD *vecdZ=0;
799 // TVectorD *veceY=0;
800 // TVectorD *veceZ=0;
801 // AliTPCLaserTrack *ltr=0;
802 // AliTPCLaserTrack::LoadTracks();
803 // tree->SetBranchAddress("dY.",&vecdY);
804 // tree->SetBranchAddress("dZ.",&vecdZ);
805 // tree->SetBranchAddress("eY.",&veceY);
806 // tree->SetBranchAddress("eZ.",&veceZ);
807 // tree->SetBranchAddress("LTr.",<r);
808 // Int_t entries= tree->GetEntries();
809 // TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
810 // Double_t bz=AliTrackerBase::GetBz();
812 // // Double_t globalXYZ[3];
813 // //Double_t globalXYZCorr[3];
814 // for (Int_t ientry=0; ientry<entries; ientry++){
815 // tree->GetEntry(ientry);
816 // if (!ltr->GetVecGX()){
817 // ltr->UpdatePoints();
820 // TVectorD fit10(5);
822 // printf("Entry\t%d\n",ientry);
823 // for (Int_t irow0=0; irow0<158; irow0+=1){
825 // TLinearFitter fitter10(4,"hyp3");
826 // TLinearFitter fitter5(2,"hyp1");
827 // Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
828 // if (sector<0) continue;
829 // //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
831 // Double_t refX= (*ltr->GetVecLX())[irow0];
832 // Int_t firstRow1 = TMath::Max(irow0-10,0);
833 // Int_t lastRow1 = TMath::Min(irow0+10,158);
834 // Double_t padWidth=(irow0<64)?0.4:0.6;
835 // // make long range fit
836 // for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
837 // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
838 // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
839 // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
840 // Double_t idealX= (*ltr->GetVecLX())[irow1];
841 // Double_t idealY= (*ltr->GetVecLY())[irow1];
842 // // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
843 // Double_t gx= (*ltr->GetVecGX())[irow1];
844 // Double_t gy= (*ltr->GetVecGY())[irow1];
845 // Double_t gz= (*ltr->GetVecGZ())[irow1];
846 // Double_t measY=(*vecdY)[irow1]+idealY;
847 // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
848 // // deltaR = R distorted -R ideal
849 // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
850 // fitter10.AddPoint(xxx,measY,1);
852 // Bool_t isOK=kTRUE;
853 // Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
854 // Double_t mean10 =0;// fitter10.GetParameter(0);
855 // Double_t slope10 =0;// fitter10.GetParameter(0);
856 // Double_t cosPart10 = 0;// fitter10.GetParameter(2);
857 // Double_t sinPart10 =0;// fitter10.GetParameter(3);
859 // if (fitter10.GetNpoints()>10){
861 // rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
862 // mean10 = fitter10.GetParameter(0);
863 // slope10 = fitter10.GetParameter(1);
864 // cosPart10 = fitter10.GetParameter(2);
865 // sinPart10 = fitter10.GetParameter(3);
867 // // make short range fit
869 // for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
870 // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
871 // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
872 // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
873 // Double_t idealX= (*ltr->GetVecLX())[irow1];
874 // Double_t idealY= (*ltr->GetVecLY())[irow1];
875 // // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
876 // Double_t gx= (*ltr->GetVecGX())[irow1];
877 // Double_t gy= (*ltr->GetVecGY())[irow1];
878 // Double_t gz= (*ltr->GetVecGZ())[irow1];
879 // Double_t measY=(*vecdY)[irow1]+idealY;
880 // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
881 // // deltaR = R distorted -R ideal
882 // Double_t expY= mean10+slope10*(idealX+deltaR-refX);
883 // if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
885 // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
886 // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
887 // fitter5.AddPoint(xxx,measY-corr,1);
892 // if (fitter5.GetNpoints()<8) isOK=kFALSE;
894 // Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
895 // Double_t offset5 =0;// fitter5.GetParameter(0);
896 // Double_t slope5 =0;// fitter5.GetParameter(0);
899 // rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
900 // offset5 = fitter5.GetParameter(0);
901 // slope5 = fitter5.GetParameter(0);
906 // Double_t phi =(*ltr->GetVecPhi())[irow0];
907 // Double_t theta =ltr->GetTgl();
908 // Double_t mean=(vecdY)->GetMatrixArray()[irow0];
909 // Double_t gx=0,gy=0,gz=0;
910 // Double_t snp = (*ltr->GetVecP2())[irow0];
911 // Int_t bundle= ltr->GetBundle();
912 // Int_t id= ltr->GetId();
913 // // Double_t rms = err->GetMatrixArray()[irow];
915 // gx = (*ltr->GetVecGX())[irow0];
916 // gy = (*ltr->GetVecGY())[irow0];
917 // gz = (*ltr->GetVecGZ())[irow0];
918 // Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
919 // fitter10.GetParameters(fit10);
920 // fitter5.GetParameters(fit5);
921 // Double_t idealY= (*ltr->GetVecLY())[irow0];
922 // Double_t measY=(*vecdY)[irow0]+idealY;
923 // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
924 // if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
926 // (*pcstream)<<"fitFull"<< // dumpe also intermediate results
927 // "bz="<<bz<< // magnetic filed used
928 // "dtype="<<dtype<< // detector match type
929 // "ptype="<<ptype<< // parameter type
930 // "theta="<<theta<< // theta
931 // "phi="<<phi<< // phi
932 // "snp="<<snp<< // snp
933 // "sector="<<sector<<
934 // "bundle="<<bundle<<
935 // // // "dsec="<<dsec<<
936 // "refX="<<refX<< // reference radius
937 // "gx="<<gx<< // global position
938 // "gy="<<gy<< // global position
939 // "gz="<<gz<< // global position
940 // "dRrec="<<dRrec<< // delta Radius in reconstruction
941 // "id="<<id<< //bundle
944 // "fit10.="<<&fit10<<
948 // "idealY="<<idealY<<
958 void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
961 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
962 // calculates partial distortions
963 // Partial distortion is stored in the resulting tree
964 // Output is storred in the file distortion_<dettype>_<partype>.root
965 // Partial distortion is stored with the name given by correction name
968 // Parameters of function:
969 // input - input tree
970 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
971 // ppype - parameter type
972 // corrArray - array with partial corrections
973 // step - skipe entries - if 1 all entries processed - it is slow
974 // debug 0 if debug on also space points dumped - it is slow
976 const Double_t kMaxSnp = 0.85;
977 const Double_t kcutSnp=0.25;
978 const Double_t kcutTheta=1.;
979 const Double_t kRadiusTPC=85;
980 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
982 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
983 // const Double_t kB2C=-0.299792458e-3;
984 const Int_t kMinEntries=20;
985 Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
988 tinput->SetBranchAddress("run",&run);
989 tinput->SetBranchAddress("theta",&theta);
990 tinput->SetBranchAddress("phi", &phi);
991 tinput->SetBranchAddress("snp",&snp);
992 tinput->SetBranchAddress("mean",&mean);
993 tinput->SetBranchAddress("rms",&rms);
994 tinput->SetBranchAddress("entries",&entries);
995 tinput->SetBranchAddress("sector",§or);
996 tinput->SetBranchAddress("dsec",&dsec);
997 tinput->SetBranchAddress("refX",&refX);
998 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1000 Int_t nentries=tinput->GetEntries();
1001 Int_t ncorr=corrArray->GetEntries();
1002 Double_t corrections[100]={0}; //
1004 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1006 if (dtype==5 || dtype==6) dtype=4;
1007 if (dtype==0) { dir=-1;}
1008 if (dtype==1) { dir=1;}
1009 if (dtype==2) { dir=-1;}
1010 if (dtype==3) { dir=1;}
1011 if (dtype==4) { dir=-1;}
1013 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1014 tinput->GetEntry(ientry);
1015 if (TMath::Abs(snp)>kMaxSnp) continue;
1018 if (dtype==2) tPar[1]=theta*kRadiusTPC;
1021 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1023 // tracks crossing CE
1024 tPar[1]=0; // track at the CE
1025 //if (TMath::Abs(theta) <0.05) continue; // deep cross
1028 if (TMath::Abs(snp) >kcutSnp) continue;
1029 if (TMath::Abs(theta) >kcutTheta) continue;
1030 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1031 Double_t bz=AliTrackerBase::GetBz();
1032 if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
1033 if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1035 if (dtype==2 && TMath::Abs(bz)>0.1 ) {
1036 tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1037 // snp at the TPC inner radius in case the vertex match used
1041 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1042 AliExternalTrackParam track(refX,phi,tPar,cov);
1046 Double_t pt=1./tPar[4];
1047 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1048 //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
1049 Double_t refXD=refX;
1050 (*pcstream)<<"fit"<<
1051 "run="<<run<< // run number
1052 "bz="<<bz<< // magnetic filed used
1053 "dtype="<<dtype<< // detector match type
1054 "ptype="<<ptype<< // parameter type
1055 "theta="<<theta<< // theta
1056 "phi="<<phi<< // phi
1057 "snp="<<snp<< // snp
1058 "mean="<<mean<< // mean dist value
1059 "rms="<<rms<< // rms
1062 "refX="<<refXD<< // referece X as double
1063 "gx="<<xyz[0]<< // global position at reference
1064 "gy="<<xyz[1]<< // global position at reference
1065 "gz="<<xyz[2]<< // global position at reference
1066 "dRrec="<<dRrec<< // delta Radius in reconstruction
1068 "id="<<id<< // track id
1069 "entries="<<entries;// number of entries in bin
1072 if (entries<kMinEntries) isOK=kFALSE;
1074 if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1075 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1076 corrections[icorr]=0;
1077 if (entries>kMinEntries){
1078 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1079 AliExternalTrackParam *trackOut = 0;
1080 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1081 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1082 if (dtype==0) {dir= -1;}
1083 if (dtype==1) {dir= 1;}
1084 if (dtype==2) {dir= -1;}
1085 if (dtype==3) {dir= 1;}
1088 if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1089 if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1090 if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1091 // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1093 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1096 corrections[icorr]=0;
1099 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
1101 (*pcstream)<<"fit"<<
1102 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1105 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1107 // special case of the TPC tracks crossing the CE
1109 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1110 corrections[icorr]=0;
1111 if (entries>kMinEntries){
1112 AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1113 AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
1114 AliExternalTrackParam *trackOut0 = 0;
1115 AliExternalTrackParam *trackOut1 = 0;
1117 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1118 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1119 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1120 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1122 if (trackOut0 && trackOut1){
1123 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1124 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1125 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1126 if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1128 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1129 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1130 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1131 if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1132 if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1134 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1135 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1137 if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1138 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1139 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1140 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1141 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1142 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1149 corrections[icorr]=0;
1153 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1155 (*pcstream)<<"fit"<<
1156 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1159 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1168 void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1171 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1172 // calculates partial distortions
1173 // Partial distortion is stored in the resulting tree
1174 // Output is storred in the file distortion_<dettype>_<partype>.root
1175 // Partial distortion is stored with the name given by correction name
1178 // Parameters of function:
1179 // input - input tree
1180 // dtype - distortion type 10 - IROC-OROC
1181 // ppype - parameter type
1182 // corrArray - array with partial corrections
1183 // step - skipe entries - if 1 all entries processed - it is slow
1184 // debug 0 if debug on also space points dumped - it is slow
1186 const Double_t kMaxSnp = 0.8;
1187 const Int_t kMinEntries=200;
1188 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1190 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1191 // const Double_t kB2C=-0.299792458e-3;
1192 Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1197 tinput->SetBranchAddress("run",&run);
1198 tinput->SetBranchAddress("theta",&theta);
1199 tinput->SetBranchAddress("phi", &phi);
1200 tinput->SetBranchAddress("snp",&snp);
1201 tinput->SetBranchAddress("mean",&mean);
1202 tinput->SetBranchAddress("rms",&rms);
1203 tinput->SetBranchAddress("entries",&entries);
1204 tinput->SetBranchAddress("sector",§or);
1205 tinput->SetBranchAddress("dsec",&dsec);
1206 tinput->SetBranchAddress("refX",&refXD);
1207 tinput->SetBranchAddress("z",&globalZ);
1208 tinput->SetBranchAddress("isec0",&isec0);
1209 tinput->SetBranchAddress("isec1",&isec1);
1210 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1212 Int_t nentries=tinput->GetEntries();
1213 Int_t ncorr=corrArray->GetEntries();
1214 Double_t corrections[100]={0}; //
1216 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1219 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1220 tinput->GetEntry(ientry);
1223 if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
1224 if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
1225 if (dtype==10 && id==-1) continue;
1232 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
1233 Double_t pt=1./tPar[4];
1235 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1236 Double_t bz=AliTrackerBase::GetBz();
1237 AliExternalTrackParam track(refX,phi,tPar,cov);
1238 Double_t xyz[3],xyzIn[3],xyzOut[3];
1240 track.GetXYZAt(85,bz,xyzIn);
1241 track.GetXYZAt(245,bz,xyzOut);
1242 Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
1243 Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1244 Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1245 Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1246 Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1247 Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1250 if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
1251 if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
1252 if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
1254 if (TMath::Abs(theta)>1) isOK=kFALSE;
1256 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1258 (*pcstream)<<"fit"<<
1260 "bz="<<bz<< // magnetic filed used
1261 "dtype="<<dtype<< // detector match type
1262 "ptype="<<ptype<< // parameter type
1263 "theta="<<theta<< // theta
1264 "phi="<<phi<< // phi
1265 "snp="<<snp<< // snp
1266 "mean="<<mean<< // mean dist value
1267 "rms="<<rms<< // rms
1270 "refX="<<refXD<< // referece X
1271 "gx="<<xyz[0]<< // global position at reference
1272 "gy="<<xyz[1]<< // global position at reference
1273 "gz="<<xyz[2]<< // global position at reference
1274 "dRrec="<<dRrec<< // delta Radius in reconstruction
1276 "id="<<id<< // track id
1277 "entries="<<entries;// number of entries in bin
1279 AliExternalTrackParam *trackOut0 = 0;
1280 AliExternalTrackParam *trackOut1 = 0;
1281 AliExternalTrackParam *ptrackIn0 = 0;
1282 AliExternalTrackParam *ptrackIn1 = 0;
1284 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1286 // special case of the TPC tracks crossing the CE
1288 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1289 corrections[icorr]=0;
1290 if (entries>kMinEntries &&isOK){
1291 AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
1292 AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
1293 ptrackIn1=&trackIn0;
1294 ptrackIn0=&trackIn1;
1296 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1297 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1298 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1299 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1301 if (trackOut0 && trackOut1){
1303 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
1304 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1305 // rotate all tracks to the same frame
1306 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1307 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1308 if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1310 if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1311 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1312 if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1314 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1315 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1316 (*pcstream)<<"fitDebug"<< // just to debug the correction
1318 "pIn0.="<<ptrackIn0<<
1319 "pIn1.="<<ptrackIn1<<
1320 "pOut0.="<<trackOut0<<
1321 "pOut1.="<<trackOut1<<
1327 corrections[icorr]=0;
1331 (*pcstream)<<"fit"<<
1332 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1335 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";