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
28 3. Space point distortion not directly observable. Instead a derived variables
29 like DCA at vertex, local y distortion in the TPC-*TOF,TRD,ITS) matching
30 in all 5 tracking parameters are obsereved.
31 In the AliTPCcorrection fir code we calculate the derivative of given variables
35 d0 = sum{ki*dO_{i}/dp_{i}} - linear fitting of the amplitudes ki
37 Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
38 extracted from the old macros
43 #include "Riostream.h"
46 #include "TGraphErrors.h"
47 #include "AliExternalTrackParam.h"
50 #include "TMultiGraph.h"
52 #include "THnSparse.h"
60 #include "AliTPCROC.h"
61 #include "AliTPCCalROC.h"
62 #include "AliESDfriend.h"
63 #include "AliTPCcalibTime.h"
64 #include "AliSplineFit.h"
65 #include "AliCDBMetaData.h"
67 #include "AliCDBManager.h"
68 #include "AliCDBStorage.h"
69 #include "AliTPCcalibBase.h"
70 #include "AliTPCcalibDB.h"
71 #include "AliTPCcalibDButil.h"
72 #include "AliRelAlignerKalman.h"
73 #include "AliTPCParamSR.h"
74 #include "AliTPCcalibTimeGain.h"
75 #include "AliTPCcalibGainMult.h"
76 #include "AliSplineFit.h"
77 #include "AliTPCComposedCorrection.h"
78 #include "AliTPCExBTwist.h"
79 #include "AliTPCCalibGlobalMisalignment.h"
80 #include "TStatToolkit.h"
83 #include "AliTrackerBase.h"
84 #include "AliTPCCorrectionFit.h"
85 #include "AliTPCLaserTrack.h"
86 #include "TDatabasePDG.h"
87 #include "AliTPCcalibAlign.h"
89 #include "AliRieman.h"
91 ClassImp(AliTPCCorrectionFit)
93 Double_t AliTPCCorrectionFit::fgMaxChi2HelixAt=1;
95 AliTPCCorrectionFit::AliTPCCorrectionFit():
96 TNamed("TPCCorrectionFit","TPCCorrectionFit")
99 // default constructor
103 AliTPCCorrectionFit::~AliTPCCorrectionFit() {
110 Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Float_t wt, Float_t t1, Float_t t2){
112 // Evaluation distortion at point using the linear approximation
114 // refX - reference X where phi and snp is calculated
115 // evalX - ref X where the distortion is evaluated
118 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
119 pcorr->SetOmegaTauT1T2(wt,t1,t2);
121 Double_t sector = 9*phi/TMath::Pi();
122 if (sector<0) sector+=18;
123 Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr);
124 Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr);
125 if (ptype==0) return y85+(y245-y85)*(evalX-85.)/(245.-85.);
126 if (ptype==2) return (y245-y85)/(245.-85.);
133 Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps,Float_t wt, Float_t t1, Float_t t2, Bool_t *pstatus){
135 // Fit the distortion along the line with the parabolic model
136 // We assume that the track are primaries - where the vertex is at (0,0,0)
139 // phi0 - phi at the entrance of the TPC
140 // snp - local inclination angle at the entrance of the TPC
141 // refX - reference X where phi and snp is calculated
142 // evalX - ref X where the distortion is evaluated
143 // theta - inclination angle
144 // corr - internal number of the correction and distortion
145 // ptype - 0 - local y distortion
146 // //1 - local z distortion
147 // 2 - local derivative distortion
149 // 4 - distortion in the curvature
150 // nsteps - number of fit points
152 // return value - distortion at point evalX with type ptype
154 static TLinearFitter fitter(3,"pol2");
155 fitter.ClearPoints();
156 if (nsteps<3) nsteps=3;
157 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
158 if (pcorr==0) return 0;
160 if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
162 Double_t deltaX=(245-85)/(nsteps);
163 Double_t curv=2.*snp/refX;
164 Double_t dPhi0=TMath::ASin(snp);
166 for (Int_t istep=0; istep<(nsteps+1); istep++){
168 Double_t localX =85.+deltaX*istep;
169 Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
170 Double_t localPhi=phi0+dPhi;
171 Double_t sector = 9*localPhi/TMath::Pi();
172 if (sector<0) sector+=18;
173 if (sector>18) sector-=18;
174 Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
175 Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
176 Double_t x[1]={localX-dlocalX};
177 fitter.AddPoint(x,dy);
181 par[0]=fitter.GetParameter(0);
182 par[1]=fitter.GetParameter(1);
183 par[2]=fitter.GetParameter(2);
185 Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
186 if (schi2> fgMaxChi2HelixAt){
187 ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",pcorr->GetName(),schi2).Data());
188 if (pstatus) (*pstatus)=kFALSE;
191 if (pstatus) (*pstatus)=kTRUE;
192 if (ptype==0) return par[0]+par[1]*evalX+par[2]*evalX*evalX;
193 if (ptype==2) return par[1]+2*par[2]*evalX;
194 if (ptype==4) return par[2];
195 if (ptype==5) return schi2;
200 Double_t AliTPCCorrectionFit::EvalAtHelix(Double_t phi0, Double_t snp, Double_t refX, Double_t evalX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps, Float_t wt, Float_t t1, Float_t t2, Bool_t *pstatus){
202 // Fit the distortion along the line with the helix model
203 // FIXME - original trajectory to be changed - AliHelix to be used
204 // We assume that the track are primaries - where the vertex is at (0,0,0)
206 // This procedure should emulate distortions as calibrated in the AliTPCcalibTime class
207 // a.) AliTPCcalibTime::FillResHistoTPCTOF - residuals, phi0 and snp evaluated at reference refX TOF clusters (refX=evalX)
208 // b.) AliTPCcalibTime::FillResHistoTPCTRD - residuals, phi0 and snp evaluated at reference refX TPCout radius (refX=evalX)
209 // c.) AliTPCcalibTime::FillResHistoTPCITS - residuals, phi0 and snp evaluated at reference refX TPCin radius (refX=evalX)
210 // d.) AliTPCcalibTime::FillResHistoTPC (vertex) - phi0 and snp evaluated at reference refX TPCin radius, residuals at vertex (refX!=evalX)
213 // phi0 - phi at the entrance of the TPC
214 // snp - local inclination angle at the entrance of the TPC
215 // refX - reference X where phi and snp is calculated
216 // evalX - ref X where the distortion is evaluated
217 // theta - inclination angle
218 // corr - internal number of the correction and distortion
219 // ptype - 0 - local y distortion
220 // //1 - local z distortion
221 // 2 - local derivative distortion
223 // 4 - distortion in the curvature
224 // nsteps - number of fit points
226 // return value - distortion at point refX with type ptype
228 static TLinearFitter fitter(3,"pol2");
229 fitter.ClearPoints();
231 if (nsteps<4) nsteps=4;
232 Double_t deltaX=(245-85)/(nsteps);
233 AliRieman rieman(nsteps);
235 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
236 if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
238 Double_t curv=2.*snp/refX;
239 Double_t dPhi0=TMath::ASin(snp);
240 for (Int_t istep=0; istep<(nsteps+1); istep++){
242 Double_t localX =85.+deltaX*istep;
243 Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
244 Double_t localPhi=phi0+dPhi;
245 Double_t sector = 9*localPhi/TMath::Pi();
246 if (sector<0) sector+=18;
247 if (sector>18) sector-=18;
248 Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
249 Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
250 Double_t x[1]={localX-dlocalX};
251 Double_t z=theta*x[0];
252 rieman.AddPoint(x[0],dy,z,0.1,0.1);
253 fitter.AddPoint(x,dy);
258 Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
259 if (schi2>fgMaxChi2HelixAt){
260 ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("Bad chi2\t%2.2f",schi2).Data());
261 if (pstatus) (*pstatus)=kFALSE;
264 if (pstatus) (*pstatus)=kTRUE;
265 if (ptype==0) return rieman.GetYat(evalX);
266 if (ptype==2) return rieman.GetDYat(evalX);
267 if (ptype==4) return rieman.GetC();
274 void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
276 // Create cluster distortion map
278 TFile *falign = TFile::Open("CalibObjects.root");
279 TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
281 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
284 AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
286 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
289 TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
291 THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
292 //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
293 AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
294 // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
296 const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
297 for (Int_t ihis=0; ihis<4; ihis++){
298 THnSparse * hisAlign =align->GetTrackletDelta(ihis);
299 AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
309 void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
311 // Make cluster residual map from the n-dimensional histogram
312 // hisInput supposed to have given format:
318 // Vertex position assumed to be at (0,0,0)
320 //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
322 Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
323 Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
324 Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
327 new TH3F("his3D","his3D",
328 nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
329 nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
330 nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
331 hisResMap3D->GetXaxis()->SetTitle("sector");
332 hisResMap3D->GetYaxis()->SetTitle("localX");
333 hisResMap3D->GetZaxis()->SetTitle("kZ");
335 TH2F * hisResMap2D[4] ={0,0,0,0};
336 for (Int_t i=0; i<4; i++){
338 new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
339 nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
340 nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
341 hisResMap2D[i]->GetXaxis()->SetTitle("sector");
342 hisResMap2D[i]->GetYaxis()->SetTitle("localX");
348 Int_t axis0[4]={0,1,2,3};
349 Int_t axis1[4]={0,1,2,3};
351 for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
353 hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
354 THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
355 Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
357 for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
360 his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
361 THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
362 Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
365 his2->GetAxis(3)->SetRangeUser(0.01,0.3);
366 TH1 * hisA = his2->Projection(0);
367 Double_t meanA= hisA->GetMean();
368 Double_t rmsA= hisA->GetRMS();
369 Double_t entriesA= hisA->GetEntries();
372 his2->GetAxis(3)->SetRangeUser(0.01,0.3);
373 TH1 * hisC = his2->Projection(0);
374 Double_t meanC= hisC->GetMean();
375 Double_t rmsC= hisC->GetRMS();
376 Double_t entriesC= hisC->GetEntries();
378 his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
379 TH2 * hisAC = his2->Projection(0,3);
380 TProfile *profAC = hisAC->ProfileX();
382 profAC->Fit("pol1","QNR","QNR",0.05,1);
383 if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
384 Double_t offsetA=f1->GetParameter(0);
385 Double_t slopeA=f1->GetParameter(1);
386 Double_t offsetAE=f1->GetParError(0);
387 Double_t slopeAE=f1->GetParError(1);
388 Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
389 profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
390 f1=(TF1*)gROOT->FindObject("pol1");
391 Double_t offsetC=f1->GetParameter(0);
392 Double_t slopeC=f1->GetParameter(1);
393 Double_t offsetCE=f1->GetParError(0);
394 Double_t slopeCE=f1->GetParError(1);
395 Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
396 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);
398 (*pcstream)<<"deltaFit"<<
403 "entriesA="<<entriesA<<
406 "entriesC="<<entriesC<<
407 "offsetA="<<offsetA<<
409 "offsetAE="<<offsetAE<<
410 "slopeAE="<<slopeAE<<
412 "offsetC="<<offsetC<<
414 "offsetCE="<<offsetCE<<
415 "slopeCE="<<slopeCE<<
419 hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
420 hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
421 hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
422 hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
424 for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
425 Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
426 if (TMath::Abs(kZ)<0.05) continue; // crossing
427 his2->GetAxis(3)->SetRange(ibin3,ibin3);
428 if (TMath::Abs(kZ)>0.15){
429 his2->GetAxis(3)->SetRange(ibin3,ibin3);
431 TH1 * his = his2->Projection(0);
432 Double_t mean= his->GetMean();
433 Double_t rms= his->GetRMS();
434 Double_t entries= his->GetEntries();
435 //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms);
436 hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
437 Double_t phi=TMath::Pi()*sector/9;
438 if (phi>TMath::Pi()) phi+=TMath::Pi();
443 his->Fit("gaus","Q","goff");
444 fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
447 his->Fit(fgaus,"Q","goff");
448 meanG=fgaus->GetParameter(1);
449 rmsG=fgaus->GetParameter(2);
452 Double_t dsec=sector-Int_t(sector)-0.5;
453 Double_t snp=dsec*TMath::Pi()/9.;
454 (*pcstream)<<"delta"<<
468 "entries="<<entries<<
471 "entriesA="<<entriesA<<
474 "entriesC="<<entriesC<<
475 "offsetA="<<offsetA<<
478 "offsetC="<<offsetC<<
488 hisResMap3D->Write();
489 hisResMap2D[0]->Write();
490 hisResMap2D[1]->Write();
491 hisResMap2D[2]->Write();
492 hisResMap2D[3]->Write();
498 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){
500 // make a distortion map out ou fthe residual histogram
501 // Results are written to the debug streamer - pcstream
503 // his0 - input (4D) residual histogram
504 // pcstream - file to write the tree
506 // refX - track matching reference X
507 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
509 // OBJ: TAxis #Delta #Delta
510 // OBJ: TAxis tanTheta tan(#Theta)
511 // OBJ: TAxis phi #phi
512 // OBJ: TAxis snp snp
514 // marian.ivanov@cern.ch
515 const Int_t kMinEntries=10;
516 Int_t idim[4]={0,1,2,3};
520 Int_t nbins3=his0->GetAxis(3)->GetNbins();
521 Int_t first3=his0->GetAxis(3)->GetFirst();
522 Int_t last3 =his0->GetAxis(3)->GetLast();
524 for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
525 his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
526 Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
527 THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
529 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
530 Int_t first2 = his3->GetAxis(2)->GetFirst();
531 Int_t last2 = his3->GetAxis(2)->GetLast();
533 for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
534 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
535 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
536 THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
537 Int_t nbins1 = his2->GetAxis(1)->GetNbins();
538 Int_t first1 = his2->GetAxis(1)->GetFirst();
539 Int_t last1 = his2->GetAxis(1)->GetLast();
540 for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
542 Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
543 his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
544 if (TMath::Abs(x1)<0.1){
545 if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
546 if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
548 if (TMath::Abs(x1)<0.06){
549 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
551 TH1 * hisDelta = his2->Projection(0);
553 Double_t entries = hisDelta->GetEntries();
554 Double_t mean=0, rms=0;
555 if (entries>kMinEntries){
556 mean = hisDelta->GetMean();
557 rms = hisDelta->GetRMS();
559 Double_t sector = 9.*x2/TMath::Pi();
560 if (sector<0) sector+=18;
561 Double_t dsec = sector-Int_t(sector)-0.5;
570 "entries="<<entries<<
573 "refX="<<refX<< // track matching refernce plane
579 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
589 void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
591 // make a distortion map out ou fthe residual histogram
592 // Results are written to the debug streamer - pcstream
594 // his0 - input (4D) residual histogram
595 // pcstream - file to write the tree
597 // refX - track matching reference X
598 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
599 // marian.ivanov@cern.ch
602 // Collection name='TObjArray', class='TObjArray', size=16
603 // 0. OBJ: TAxis #Delta #Delta
604 // 1. OBJ: TAxis N_{cl} N_{cl}
605 // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm)
606 // 3. OBJ: TAxis z (cm) z (cm)
607 // 4. OBJ: TAxis sin(#phi) sin(#phi)
608 // 5. OBJ: TAxis tan(#theta) tan(#theta)
609 // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV)
610 // 7. OBJ: TAxis pt (GeV) pt (GeV)
611 // 8. OBJ: TAxis alpha alpha
612 const Int_t kMinEntries=10;
614 // 1. make default selections
617 Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
618 hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
619 hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
620 hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
621 hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
622 hisDelta= hisInput->Projection(0);
623 hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
625 THnSparse *his0= hisInput->Projection(4,idim0);
627 // 2. Get mean in diferent bins
629 Int_t nbins1=his0->GetAxis(1)->GetNbins();
630 Int_t first1=his0->GetAxis(1)->GetFirst();
631 Int_t last1 =his0->GetAxis(1)->GetLast();
633 Double_t bz=AliTrackerBase::GetBz();
634 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
636 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
638 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
639 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
641 THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
642 Int_t nbins3 = his1->GetAxis(3)->GetNbins();
643 Int_t first3 = his1->GetAxis(3)->GetFirst();
644 Int_t last3 = his1->GetAxis(3)->GetLast();
646 for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
647 his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
648 Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
650 his1->GetAxis(3)->SetRangeUser(-1,1);
653 THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
654 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
655 Int_t first2 = his3->GetAxis(2)->GetFirst();
656 Int_t last2 = his3->GetAxis(2)->GetLast();
658 for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
659 his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
660 Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
661 hisDelta = his3->Projection(0);
663 Double_t entries = hisDelta->GetEntries();
664 Double_t mean=0, rms=0;
665 if (entries>kMinEntries){
666 mean = hisDelta->GetMean();
667 rms = hisDelta->GetRMS();
669 Double_t sector = 9.*x2/TMath::Pi();
670 if (sector<0) sector+=18;
671 Double_t dsec = sector-Int_t(sector)-0.5;
672 Double_t snp=0; // dummy snp - equal 0
675 "bz="<<bz<< // magnetic field
676 "theta="<<x1<< // theta
677 "phi="<<x2<< // phi (alpha)
678 "z="<<x3<< // z at "vertex"
679 "snp="<<snp<< // dummy snp
680 "entries="<<entries<< // entries in bin
681 "mean="<<mean<< // mean
683 "refX="<<refX<< // track matching refernce plane
684 "type="<<type<< // parameter type
685 "sector="<<sector<< // sector
686 "dsec="<<dsec<< // dummy delta sector
689 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
700 void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
702 // make a distortion map out of the residual histogram
703 // Results are written to the debug streamer - pcstream
705 // his0 - input (4D) residual histogram
706 // pcstream - file to write the tree
708 // type - 0- y 1-z,2 -snp, 3-theta
709 // bz - magnetic field
710 // marian.ivanov@cern.ch
712 //Collection name='TObjArray', class='TObjArray', size=16
713 //0 OBJ: TAxis delta delta
714 //1 OBJ: TAxis phi phi
715 //2 OBJ: TAxis localX localX
718 //5 OBJ: TAxis is1 is1
719 //6 OBJ: TAxis is0 is0
721 //8. OBJ: TAxis IsPrimary IsPrimary
723 const Int_t kMinEntries=10;
724 THnSparse * hisSector0=0;
725 TH1 * htemp=0; // histogram to calculate mean value of parameter
726 // Double_t bz=AliTrackerBase::GetBz();
729 // Loop over pair of sector:
740 for (Int_t isec0=0; isec0<72; isec0++){
741 Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces
743 //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ?
744 hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
745 hisSector0=hisInput->Projection(7,index0);
748 for (Int_t isec1=isec0+1; isec1<72; isec1++){
749 //if (isec1!=isec0+36) continue;
750 if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
751 printf("Sectors %d\t%d\n",isec1,isec0);
752 hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
753 TH1 * hisX=hisSector0->Projection(5);
754 Double_t refX= hisX->GetMean();
756 TH1 *hisDelta=hisSector0->Projection(0);
757 Double_t dmean = hisDelta->GetMean();
758 Double_t drms = hisDelta->GetRMS();
759 hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
762 // 1. make default selections
764 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
765 THnBase *hisSector1= hisSector0->ProjectionND(5,idim0);
767 // 2. Get mean in diferent bins
769 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
771 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
772 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
773 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
775 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi
779 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
780 Double_t psec = (9*xPhi/TMath::Pi());
781 if (psec<0) psec+=18;
784 if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
785 if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
786 if (!isOK0) continue;
787 if (!isOK1) continue;
789 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
790 if (isec1!=isec0+36) {
791 hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
794 htemp = hisSector1->Projection(4);
795 xPhi=htemp->GetMean();
797 THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
798 //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
799 Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
800 Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
802 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z
806 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
807 if (isec1!=isec0+36) {
808 hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
810 htemp = hisPhi->Projection(3);
811 Double_t xZ= htemp->GetMean();
813 THnBase * hisZ= hisPhi->ProjectionND(3,idim);
814 //projected histogram according selection 3 -z
817 //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
818 Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
819 Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
820 for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
824 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
825 if (isec1!=isec0+36) {
826 hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
828 htemp = hisZ->Projection(2);
829 Double_t xSnp= htemp->GetMean();
831 THnBase * hisSnp= hisZ->ProjectionND(2,idim);
832 //projected histogram according selection 2 - snp
834 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
835 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
836 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
838 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
841 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
842 if (isec1!=isec0+36) {
843 hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
845 htemp = hisSnp->Projection(1);
846 Double_t xTheta=htemp->GetMean();
848 hisDelta = hisSnp->Projection(0);
850 Double_t entries = hisDelta->GetEntries();
851 Double_t mean=0, rms=0;
852 if (entries>kMinEntries){
853 mean = hisDelta->GetMean();
854 rms = hisDelta->GetRMS();
856 Double_t sector = 9.*xPhi/TMath::Pi();
857 if (sector<0) sector+=18;
858 Double_t dsec = sector-Int_t(sector)-0.5;
859 Int_t dtype=1; // TPC alignment type
862 "bz="<<bz<< // magnetic field
863 "ptype="<<type<< // parameter type
864 "dtype="<<dtype<< // parameter type
865 "isec0="<<isec0<< // sector 0
866 "isec1="<<isec1<< // sector 1
867 "sector="<<sector<< // sector as float
868 "dsec="<<dsec<< // delta sector
870 "theta="<<xTheta<< // theta
871 "phi="<<xPhi<< // phi (alpha)
873 "snp="<<xSnp<< // snp
876 "entries="<<entries<< // entries in bin
877 "mean="<<mean<< // mean
879 "refX="<<refX<< // track matching reference plane
882 //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);
901 // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
903 // // Make a laser fit tree for global minimization
905 // AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
906 // AliTPCCorrection * correction = calib->GetTPCComposedCorrection();
907 // if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz());
908 // correction->AddVisualCorrection(correction,0); //register correction
910 // // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
911 // //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
913 // const Double_t cutErrY=0.05;
914 // const Double_t kSigmaCut=4;
915 // // const Double_t cutErrZ=0.03;
916 // const Double_t kEpsilon=0.00000001;
917 // // const Double_t kMaxDist=1.; // max distance - space correction
918 // TVectorD *vecdY=0;
919 // TVectorD *vecdZ=0;
920 // TVectorD *veceY=0;
921 // TVectorD *veceZ=0;
922 // AliTPCLaserTrack *ltr=0;
923 // AliTPCLaserTrack::LoadTracks();
924 // tree->SetBranchAddress("dY.",&vecdY);
925 // tree->SetBranchAddress("dZ.",&vecdZ);
926 // tree->SetBranchAddress("eY.",&veceY);
927 // tree->SetBranchAddress("eZ.",&veceZ);
928 // tree->SetBranchAddress("LTr.",<r);
929 // Int_t entries= tree->GetEntries();
930 // TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
931 // Double_t bz=AliTrackerBase::GetBz();
933 // // Double_t globalXYZ[3];
934 // //Double_t globalXYZCorr[3];
935 // for (Int_t ientry=0; ientry<entries; ientry++){
936 // tree->GetEntry(ientry);
937 // if (!ltr->GetVecGX()){
938 // ltr->UpdatePoints();
941 // TVectorD fit10(5);
943 // printf("Entry\t%d\n",ientry);
944 // for (Int_t irow0=0; irow0<158; irow0+=1){
946 // TLinearFitter fitter10(4,"hyp3");
947 // TLinearFitter fitter5(2,"hyp1");
948 // Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0];
949 // if (sector<0) continue;
950 // //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue;
952 // Double_t refX= (*ltr->GetVecLX())[irow0];
953 // Int_t firstRow1 = TMath::Max(irow0-10,0);
954 // Int_t lastRow1 = TMath::Min(irow0+10,158);
955 // Double_t padWidth=(irow0<64)?0.4:0.6;
956 // // make long range fit
957 // for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){
958 // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
959 // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
960 // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
961 // Double_t idealX= (*ltr->GetVecLX())[irow1];
962 // Double_t idealY= (*ltr->GetVecLY())[irow1];
963 // // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
964 // Double_t gx= (*ltr->GetVecGX())[irow1];
965 // Double_t gy= (*ltr->GetVecGY())[irow1];
966 // Double_t gz= (*ltr->GetVecGZ())[irow1];
967 // Double_t measY=(*vecdY)[irow1]+idealY;
968 // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
969 // // deltaR = R distorted -R ideal
970 // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
971 // fitter10.AddPoint(xxx,measY,1);
973 // Bool_t isOK=kTRUE;
974 // Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
975 // Double_t mean10 =0;// fitter10.GetParameter(0);
976 // Double_t slope10 =0;// fitter10.GetParameter(0);
977 // Double_t cosPart10 = 0;// fitter10.GetParameter(2);
978 // Double_t sinPart10 =0;// fitter10.GetParameter(3);
980 // if (fitter10.GetNpoints()>10){
982 // rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4));
983 // mean10 = fitter10.GetParameter(0);
984 // slope10 = fitter10.GetParameter(1);
985 // cosPart10 = fitter10.GetParameter(2);
986 // sinPart10 = fitter10.GetParameter(3);
988 // // make short range fit
990 // for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){
991 // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue;
992 // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue;
993 // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue;
994 // Double_t idealX= (*ltr->GetVecLX())[irow1];
995 // Double_t idealY= (*ltr->GetVecLY())[irow1];
996 // // Double_t idealZ= (*ltr->GetVecLZ())[irow1];
997 // Double_t gx= (*ltr->GetVecGX())[irow1];
998 // Double_t gy= (*ltr->GetVecGY())[irow1];
999 // Double_t gz= (*ltr->GetVecGZ())[irow1];
1000 // Double_t measY=(*vecdY)[irow1]+idealY;
1001 // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0);
1002 // // deltaR = R distorted -R ideal
1003 // Double_t expY= mean10+slope10*(idealX+deltaR-refX);
1004 // if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue;
1006 // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
1007 // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)};
1008 // fitter5.AddPoint(xxx,measY-corr,1);
1013 // if (fitter5.GetNpoints()<8) isOK=kFALSE;
1015 // Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
1016 // Double_t offset5 =0;// fitter5.GetParameter(0);
1017 // Double_t slope5 =0;// fitter5.GetParameter(0);
1020 // rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
1021 // offset5 = fitter5.GetParameter(0);
1022 // slope5 = fitter5.GetParameter(0);
1025 // Double_t dtype=5;
1026 // Double_t ptype=0;
1027 // Double_t phi =(*ltr->GetVecPhi())[irow0];
1028 // Double_t theta =ltr->GetTgl();
1029 // Double_t mean=(vecdY)->GetMatrixArray()[irow0];
1030 // Double_t gx=0,gy=0,gz=0;
1031 // Double_t snp = (*ltr->GetVecP2())[irow0];
1032 // Int_t bundle= ltr->GetBundle();
1033 // Int_t id= ltr->GetId();
1034 // // Double_t rms = err->GetMatrixArray()[irow];
1036 // gx = (*ltr->GetVecGX())[irow0];
1037 // gy = (*ltr->GetVecGY())[irow0];
1038 // gz = (*ltr->GetVecGZ())[irow0];
1039 // Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0);
1040 // fitter10.GetParameters(fit10);
1041 // fitter5.GetParameters(fit5);
1042 // Double_t idealY= (*ltr->GetVecLY())[irow0];
1043 // Double_t measY=(*vecdY)[irow0]+idealY;
1044 // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth);
1045 // if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE;
1047 // (*pcstream)<<"fitFull"<< // dumpe also intermediate results
1048 // "bz="<<bz<< // magnetic filed used
1049 // "dtype="<<dtype<< // detector match type
1050 // "ptype="<<ptype<< // parameter type
1051 // "theta="<<theta<< // theta
1052 // "phi="<<phi<< // phi
1053 // "snp="<<snp<< // snp
1054 // "sector="<<sector<<
1055 // "bundle="<<bundle<<
1056 // // // "dsec="<<dsec<<
1057 // "refX="<<refX<< // reference radius
1058 // "gx="<<gx<< // global position
1059 // "gy="<<gy<< // global position
1060 // "gz="<<gz<< // global position
1061 // "dRrec="<<dRrec<< // delta Radius in reconstruction
1062 // "id="<<id<< //bundle
1063 // "rms10="<<rms10<<
1065 // "fit10.="<<&fit10<<
1066 // "fit5.="<<&fit5<<
1067 // "measY="<<measY<<
1069 // "idealY="<<idealY<<
1079 void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1082 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1083 // calculates partial distortions
1084 // Partial distortion is stored in the resulting tree
1085 // Output is storred in the file distortion_<dettype>_<partype>.root
1086 // Partial distortion is stored with the name given by correction name
1089 // Parameters of function:
1090 // input - input tree
1091 // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing
1092 // ppype - parameter type
1093 // corrArray - array with partial corrections
1094 // step - skipe entries - if 1 all entries processed - it is slow
1095 // debug 0 if debug on also space points dumped - it is slow
1097 const Double_t kMaxSnp = 0.85;
1098 const Double_t kcutSnp=0.25;
1099 const Double_t kcutTheta=1.;
1100 const Double_t kRadiusTPC=85;
1101 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1103 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1104 // const Double_t kB2C=-0.299792458e-3;
1105 const Int_t kMinEntries=20;
1106 Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
1109 tinput->SetBranchAddress("run",&run);
1110 tinput->SetBranchAddress("theta",&theta);
1111 tinput->SetBranchAddress("phi", &phi);
1112 tinput->SetBranchAddress("snp",&snp);
1113 tinput->SetBranchAddress("mean",&mean);
1114 tinput->SetBranchAddress("rms",&rms);
1115 tinput->SetBranchAddress("entries",&entries);
1116 tinput->SetBranchAddress("sector",§or);
1117 tinput->SetBranchAddress("dsec",&dsec);
1118 tinput->SetBranchAddress("refX",&refX);
1119 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1121 Int_t nentries=tinput->GetEntries();
1122 Int_t ncorr=corrArray->GetEntries();
1123 Double_t corrections[100]={0}; //
1125 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1127 if (dtype==5 || dtype==6) dtype=4;
1128 if (dtype==0) { dir=-1;}
1129 if (dtype==1) { dir=1;}
1130 if (dtype==2) { dir=-1;}
1131 if (dtype==3) { dir=1;}
1132 if (dtype==4) { dir=-1;}
1134 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1135 tinput->GetEntry(ientry);
1136 if (TMath::Abs(snp)>kMaxSnp) continue;
1139 if (dtype==2) tPar[1]=theta*kRadiusTPC;
1142 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1144 // tracks crossing CE
1145 tPar[1]=0; // track at the CE
1146 //if (TMath::Abs(theta) <0.05) continue; // deep cross
1149 if (TMath::Abs(snp) >kcutSnp) continue;
1150 if (TMath::Abs(theta) >kcutTheta) continue;
1151 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1152 Double_t bz=AliTrackerBase::GetBz();
1153 if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
1154 if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1156 if (dtype==2 && TMath::Abs(bz)>0.1 ) {
1157 tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1158 // snp at the TPC inner radius in case the vertex match used
1162 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1163 AliExternalTrackParam track(refX,phi,tPar,cov);
1167 Double_t pt=1./tPar[4];
1168 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1169 //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used
1170 Double_t refXD=refX;
1171 (*pcstream)<<"fit"<<
1172 "run="<<run<< // run number
1173 "bz="<<bz<< // magnetic filed used
1174 "dtype="<<dtype<< // detector match type
1175 "ptype="<<ptype<< // parameter type
1176 "theta="<<theta<< // theta
1177 "phi="<<phi<< // phi
1178 "snp="<<snp<< // snp
1179 "mean="<<mean<< // mean dist value
1180 "rms="<<rms<< // rms
1183 "refX="<<refXD<< // referece X as double
1184 "gx="<<xyz[0]<< // global position at reference
1185 "gy="<<xyz[1]<< // global position at reference
1186 "gz="<<xyz[2]<< // global position at reference
1187 "dRrec="<<dRrec<< // delta Radius in reconstruction
1189 "id="<<id<< // track id
1190 "entries="<<entries;// number of entries in bin
1193 if (entries<kMinEntries) isOK=kFALSE;
1195 if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1196 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1197 corrections[icorr]=0;
1198 if (entries>kMinEntries){
1199 AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1200 AliExternalTrackParam *trackOut = 0;
1201 if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1202 if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1203 if (dtype==0) {dir= -1;}
1204 if (dtype==1) {dir= 1;}
1205 if (dtype==2) {dir= -1;}
1206 if (dtype==3) {dir= 1;}
1209 if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1210 if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1211 if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1212 // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1214 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1217 corrections[icorr]=0;
1220 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
1222 (*pcstream)<<"fit"<<
1223 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1226 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1228 // special case of the TPC tracks crossing the CE
1230 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1231 corrections[icorr]=0;
1232 if (entries>kMinEntries){
1233 AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1234 AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
1235 AliExternalTrackParam *trackOut0 = 0;
1236 AliExternalTrackParam *trackOut1 = 0;
1238 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1239 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1240 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1241 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1243 if (trackOut0 && trackOut1){
1244 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1245 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1246 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1247 if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1249 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1250 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1251 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1252 if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1253 if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1255 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1256 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1258 if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1259 (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1260 (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1261 (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1262 (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1263 (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1270 corrections[icorr]=0;
1274 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1276 (*pcstream)<<"fit"<<
1277 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1280 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1289 void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1292 // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX)
1293 // calculates partial distortions
1294 // Partial distortion is stored in the resulting tree
1295 // Output is storred in the file distortion_<dettype>_<partype>.root
1296 // Partial distortion is stored with the name given by correction name
1299 // Parameters of function:
1300 // input - input tree
1301 // dtype - distortion type 10 - IROC-OROC
1302 // ppype - parameter type
1303 // corrArray - array with partial corrections
1304 // step - skipe entries - if 1 all entries processed - it is slow
1305 // debug 0 if debug on also space points dumped - it is slow
1307 const Double_t kMaxSnp = 0.8;
1308 const Int_t kMinEntries=200;
1309 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1311 const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1312 // const Double_t kB2C=-0.299792458e-3;
1313 Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1318 tinput->SetBranchAddress("run",&run);
1319 tinput->SetBranchAddress("theta",&theta);
1320 tinput->SetBranchAddress("phi", &phi);
1321 tinput->SetBranchAddress("snp",&snp);
1322 tinput->SetBranchAddress("mean",&mean);
1323 tinput->SetBranchAddress("rms",&rms);
1324 tinput->SetBranchAddress("entries",&entries);
1325 tinput->SetBranchAddress("sector",§or);
1326 tinput->SetBranchAddress("dsec",&dsec);
1327 tinput->SetBranchAddress("refX",&refXD);
1328 tinput->SetBranchAddress("z",&globalZ);
1329 tinput->SetBranchAddress("isec0",&isec0);
1330 tinput->SetBranchAddress("isec1",&isec1);
1331 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1333 Int_t nentries=tinput->GetEntries();
1334 Int_t ncorr=corrArray->GetEntries();
1335 Double_t corrections[100]={0}; //
1337 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1340 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1341 tinput->GetEntry(ientry);
1344 if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
1345 if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
1346 if (dtype==10 && id==-1) continue;
1353 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
1354 Double_t pt=1./tPar[4];
1356 printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1357 Double_t bz=AliTrackerBase::GetBz();
1358 AliExternalTrackParam track(refX,phi,tPar,cov);
1359 Double_t xyz[3],xyzIn[3],xyzOut[3];
1361 track.GetXYZAt(85,bz,xyzIn);
1362 track.GetXYZAt(245,bz,xyzOut);
1363 Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
1364 Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1365 Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1366 Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1367 Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1368 Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1371 if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
1372 if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
1373 if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
1375 if (TMath::Abs(theta)>1) isOK=kFALSE;
1377 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1379 (*pcstream)<<"fit"<<
1381 "bz="<<bz<< // magnetic filed used
1382 "dtype="<<dtype<< // detector match type
1383 "ptype="<<ptype<< // parameter type
1384 "theta="<<theta<< // theta
1385 "phi="<<phi<< // phi
1386 "snp="<<snp<< // snp
1387 "mean="<<mean<< // mean dist value
1388 "rms="<<rms<< // rms
1391 "refX="<<refXD<< // referece X
1392 "gx="<<xyz[0]<< // global position at reference
1393 "gy="<<xyz[1]<< // global position at reference
1394 "gz="<<xyz[2]<< // global position at reference
1395 "dRrec="<<dRrec<< // delta Radius in reconstruction
1397 "id="<<id<< // track id
1398 "entries="<<entries;// number of entries in bin
1400 AliExternalTrackParam *trackOut0 = 0;
1401 AliExternalTrackParam *trackOut1 = 0;
1402 AliExternalTrackParam *ptrackIn0 = 0;
1403 AliExternalTrackParam *ptrackIn1 = 0;
1405 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1407 // special case of the TPC tracks crossing the CE
1409 AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1410 corrections[icorr]=0;
1411 if (entries>kMinEntries &&isOK){
1412 AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
1413 AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
1414 ptrackIn1=&trackIn0;
1415 ptrackIn0=&trackIn1;
1417 if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1418 if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1419 if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1420 if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1422 if (trackOut0 && trackOut1){
1424 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
1425 if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1426 // rotate all tracks to the same frame
1427 if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1428 if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1429 if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1431 if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1432 if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1433 if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1435 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1436 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1437 (*pcstream)<<"fitDebug"<< // just to debug the correction
1439 "pIn0.="<<ptrackIn0<<
1440 "pIn1.="<<ptrackIn1<<
1441 "pOut0.="<<trackOut0<<
1442 "pOut1.="<<trackOut1<<
1448 corrections[icorr]=0;
1452 (*pcstream)<<"fit"<<
1453 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1456 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";