]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/TPCcalib/AliTPCCorrectionFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / TPC / TPCcalib / AliTPCCorrectionFit.cxx
CommitLineData
93096a07 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
18/*
19 Responsible: marian.ivanov@cern.ch
20 Tools for fitting of the space point distortion parameters.
21 Functionality
22
23
241. Creation of the distortion maps from the residual histograms
1cadb1b2 25
93096a07 262. Making fit trees
1cadb1b2 27
283. 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
32 dO_{i}/dp_{i}
33
344. Global fit - later
35 d0 = sum{ki*dO_{i}/dp_{i}} - linear fitting of the amplitudes ki
93096a07 36
37 Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
38 extracted from the old macros
39
40
41*/
42
43#include "Riostream.h"
44#include <fstream>
45#include "TMap.h"
46#include "TGraphErrors.h"
47#include "AliExternalTrackParam.h"
48#include "TFile.h"
49#include "TGraph.h"
50#include "TMultiGraph.h"
51#include "TCanvas.h"
52#include "THnSparse.h"
53#include "THnBase.h"
54#include "TProfile.h"
55#include "TROOT.h"
56#include "TLegend.h"
57#include "TPad.h"
58#include "TH2D.h"
59#include "TH3D.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"
66#include "AliCDBId.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"
81#include "TChain.h"
82#include "TCut.h"
83#include "AliTrackerBase.h"
84#include "AliTPCCorrectionFit.h"
85#include "AliTPCLaserTrack.h"
86#include "TDatabasePDG.h"
87#include "AliTPCcalibAlign.h"
4425a337 88#include "AliLog.h"
1cadb1b2 89#include "AliRieman.h"
93096a07 90
91ClassImp(AliTPCCorrectionFit)
92
1636786f 93Double_t AliTPCCorrectionFit::fgMaxChi2HelixAt=1;
94
93096a07 95AliTPCCorrectionFit::AliTPCCorrectionFit():
96 TNamed("TPCCorrectionFit","TPCCorrectionFit")
97{
98 //
99 // default constructor
100 //
101}
102
103AliTPCCorrectionFit::~AliTPCCorrectionFit() {
104 //
105 // Destructor
106 //
107}
108
4f1f7ba7 109/*
110Double_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){
93096a07 111 //
4f1f7ba7 112 // Evaluation distortion at point using the linear approximation
93096a07 113 //
4f1f7ba7 114 // refX - reference X where phi and snp is calculated
115 // evalX - ref X where the distortion is evaluated
116
1636786f 117 if (wt<50){
118 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
119 pcorr->SetOmegaTauT1T2(wt,t1,t2);
120 }
93096a07 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);
4f1f7ba7 125 if (ptype==0) return y85+(y245-y85)*(evalX-85.)/(245.-85.);
93096a07 126 if (ptype==2) return (y245-y85)/(245.-85.);
127 return 0;
128}
4f1f7ba7 129*/
93096a07 130
131
132
155e9dd0 133Double_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){
93096a07 134 //
135 // Fit the distortion along the line with the parabolic model
1cadb1b2 136 // We assume that the track are primaries - where the vertex is at (0,0,0)
137 //
93096a07 138 // Parameters:
1cadb1b2 139 // phi0 - phi at the entrance of the TPC
140 // snp - local inclination angle at the entrance of the TPC
4f1f7ba7 141 // refX - reference X where phi and snp is calculated
142 // evalX - ref X where the distortion is evaluated
1cadb1b2 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
148 // //3
149 // 4 - distortion in the curvature
150 // nsteps - number of fit points
151 //
4f1f7ba7 152 // return value - distortion at point evalX with type ptype
1cadb1b2 153 //
93096a07 154 static TLinearFitter fitter(3,"pol2");
155 fitter.ClearPoints();
156 if (nsteps<3) nsteps=3;
1636786f 157 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
4f1f7ba7 158 if (pcorr==0) return 0;
1636786f 159 if (wt<50){
155e9dd0 160 if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
1636786f 161 }
93096a07 162 Double_t deltaX=(245-85)/(nsteps);
4f1f7ba7 163 Double_t curv=2.*snp/refX;
164 Double_t dPhi0=TMath::ASin(snp);
165
93096a07 166 for (Int_t istep=0; istep<(nsteps+1); istep++){
167 //
168 Double_t localX =85.+deltaX*istep;
4f1f7ba7 169 Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
170 Double_t localPhi=phi0+dPhi;
93096a07 171 Double_t sector = 9*localPhi/TMath::Pi();
172 if (sector<0) sector+=18;
4f1f7ba7 173 if (sector>18) sector-=18;
1cadb1b2 174 Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
93096a07 175 Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
176 Double_t x[1]={localX-dlocalX};
1cadb1b2 177 fitter.AddPoint(x,dy);
93096a07 178 }
179 fitter.Eval();
180 Double_t par[3];
181 par[0]=fitter.GetParameter(0);
182 par[1]=fitter.GetParameter(1);
183 par[2]=fitter.GetParameter(2);
4f1f7ba7 184
1636786f 185 Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
186 if (schi2> fgMaxChi2HelixAt){
4f1f7ba7 187 ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",pcorr->GetName(),schi2).Data());
155e9dd0 188 if (pstatus) (*pstatus)=kFALSE;
1636786f 189 return 0;
190 }
155e9dd0 191 if (pstatus) (*pstatus)=kTRUE;
4f1f7ba7 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;
93096a07 194 if (ptype==4) return par[2];
1636786f 195 if (ptype==5) return schi2;
93096a07 196 return 0;
197}
198
199
155e9dd0 200Double_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){
1cadb1b2 201 //
202 // Fit the distortion along the line with the helix model
2cc73087 203 // FIXME - original trajectory to be changed - AliHelix to be used
1cadb1b2 204 // We assume that the track are primaries - where the vertex is at (0,0,0)
205 //
4f1f7ba7 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)
211 // *
1cadb1b2 212 // Parameters:
213 // phi0 - phi at the entrance of the TPC
214 // snp - local inclination angle at the entrance of the TPC
4f1f7ba7 215 // refX - reference X where phi and snp is calculated
216 // evalX - ref X where the distortion is evaluated
1cadb1b2 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
222 // //3
223 // 4 - distortion in the curvature
224 // nsteps - number of fit points
225 //
226 // return value - distortion at point refX with type ptype
227 //
1636786f 228 static TLinearFitter fitter(3,"pol2");
229 fitter.ClearPoints();
4f1f7ba7 230 //
1636786f 231 if (nsteps<4) nsteps=4;
1cadb1b2 232 Double_t deltaX=(245-85)/(nsteps);
233 AliRieman rieman(nsteps);
1636786f 234 if (wt<50){
155e9dd0 235 AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
236 if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
1636786f 237 }
4f1f7ba7 238 Double_t curv=2.*snp/refX;
239 Double_t dPhi0=TMath::ASin(snp);
1cadb1b2 240 for (Int_t istep=0; istep<(nsteps+1); istep++){
241 //
242 Double_t localX =85.+deltaX*istep;
4f1f7ba7 243 Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
244 Double_t localPhi=phi0+dPhi;
1cadb1b2 245 Double_t sector = 9*localPhi/TMath::Pi();
246 if (sector<0) sector+=18;
4f1f7ba7 247 if (sector>18) sector-=18;
1cadb1b2 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);
1636786f 253 fitter.AddPoint(x,dy);
1cadb1b2 254 }
1636786f 255 fitter.Eval();
1cadb1b2 256 rieman.Update();
257 //
1636786f 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());
155e9dd0 261 if (pstatus) (*pstatus)=kFALSE;
1636786f 262 return 0;
263 }
155e9dd0 264 if (pstatus) (*pstatus)=kTRUE;
4f1f7ba7 265 if (ptype==0) return rieman.GetYat(evalX);
266 if (ptype==2) return rieman.GetDYat(evalX);
1cadb1b2 267 if (ptype==4) return rieman.GetC();
268 return 0;
269}
270
271
93096a07 272
273
274void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
275 //
276 // Create cluster distortion map
277 //
278 TFile *falign = TFile::Open("CalibObjects.root");
279 TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
63d14e61 280 if (!arrayAlign) {
4425a337 281 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
63d14e61 282 return;
283 }
93096a07 284 AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
63d14e61 285 if (!align) {
4425a337 286 AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
63d14e61 287 return;
288 }
93096a07 289 TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
290
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);
295
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);
300 }
301 delete pcstream;
302 delete falign;
303}
304
305
306
307
308
309void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){
310 //
311 // Make cluster residual map from the n-dimensional histogram
312 // hisInput supposed to have given format:
313 // 4 Dim:
314 // delta,
315 // sector
316 // localX
317 // kZ
318 // Vertex position assumed to be at (0,0,0)
319 //
320 //TTreeSRedirector *pcstream=new TTreeSRedirector(sname);
321 //
322 Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
323 Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
324 Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
325 TF1 *fgaus=0;
326 TH3F * hisResMap3D =
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");
334
335 TH2F * hisResMap2D[4] ={0,0,0,0};
336 for (Int_t i=0; i<4; i++){
337 hisResMap2D[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");
343 }
344 //
345 //
346 //
347 TF1 * f1= 0;
348 Int_t axis0[4]={0,1,2,3};
349 Int_t axis1[4]={0,1,2,3};
350 Int_t counter=0;
351 for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
352 // phi- sector range
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);
356 //
357 for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
358 // local x range
359 // kz fits
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);
363 //
364 //A side
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();
370 delete hisA;
371 //C side
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();
377 delete hisC;
378 his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
379 TH2 * hisAC = his2->Projection(0,3);
380 TProfile *profAC = hisAC->ProfileX();
381 delete hisAC;
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);
397 counter++;
398 (*pcstream)<<"deltaFit"<<
399 "sector="<<sector<<
400 "localX="<<localX<<
401 "meanA="<<meanA<<
402 "rmsA="<<rmsA<<
403 "entriesA="<<entriesA<<
404 "meanC="<<meanC<<
405 "rmsC="<<rmsC<<
406 "entriesC="<<entriesC<<
407 "offsetA="<<offsetA<<
408 "slopeA="<<slopeA<<
409 "offsetAE="<<offsetAE<<
410 "slopeAE="<<slopeAE<<
411 "chi2A="<<chi2A<<
412 "offsetC="<<offsetC<<
413 "slopeC="<<slopeC<<
414 "offsetCE="<<offsetCE<<
415 "slopeCE="<<slopeCE<<
416 "chi2C="<<chi2C<<
417 "\n";
418 //
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);
423
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);
430 }
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();
439 Double_t meanG=0;
440 Double_t rmsG=0;
441 if (entries>50){
442 if (!fgaus) {
443 his->Fit("gaus","Q","goff");
444 fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
445 }
446 if (fgaus) {
447 his->Fit(fgaus,"Q","goff");
448 meanG=fgaus->GetParameter(1);
449 rmsG=fgaus->GetParameter(2);
450 }
451 }
452 Double_t dsec=sector-Int_t(sector)-0.5;
453 Double_t snp=dsec*TMath::Pi()/9.;
454 (*pcstream)<<"delta"<<
455 "ptype="<<ptype<<
456 "dtype="<<dtype<<
457 "sector="<<sector<<
458 "dsec="<<dsec<<
459 "snp="<<snp<<
460 "phi="<<phi<<
461 "localX="<<localX<<
462 "kZ="<<kZ<<
463 "theta="<<kZ<<
464 "mean="<<mean<<
465 "rms="<<rms<<
466 "meanG="<<meanG<<
467 "rmsG="<<rmsG<<
468 "entries="<<entries<<
469 "meanA="<<meanA<<
470 "rmsA="<<rmsA<<
471 "entriesA="<<entriesA<<
472 "meanC="<<meanC<<
473 "rmsC="<<rmsC<<
474 "entriesC="<<entriesC<<
475 "offsetA="<<offsetA<<
476 "slopeA="<<slopeA<<
477 "chi2A="<<chi2A<<
478 "offsetC="<<offsetC<<
479 "slopeC="<<slopeC<<
480 "chi2C="<<chi2C<<
481 "\n";
482 delete his;
483 }
484 delete his2;
485 }
486 delete his1;
487 }
488 hisResMap3D->Write();
489 hisResMap2D[0]->Write();
490 hisResMap2D[1]->Write();
491 hisResMap2D[2]->Write();
492 hisResMap2D[3]->Write();
493 // delete pcstream;
494}
495
496
497
498void 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){
499 //
500 // make a distortion map out ou fthe residual histogram
501 // Results are written to the debug streamer - pcstream
502 // Parameters:
503 // his0 - input (4D) residual histogram
504 // pcstream - file to write the tree
505 // run - run number
506 // refX - track matching reference X
507 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
508 // THnSparse axes:
509 // OBJ: TAxis #Delta #Delta
510 // OBJ: TAxis tanTheta tan(#Theta)
511 // OBJ: TAxis phi #phi
512 // OBJ: TAxis snp snp
513
514 // marian.ivanov@cern.ch
515 const Int_t kMinEntries=10;
516 Int_t idim[4]={0,1,2,3};
517 //
518 //
519 //
520 Int_t nbins3=his0->GetAxis(3)->GetNbins();
521 Int_t first3=his0->GetAxis(3)->GetFirst();
522 Int_t last3 =his0->GetAxis(3)->GetLast();
523 //
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
528 //
529 Int_t nbins2 = his3->GetAxis(2)->GetNbins();
530 Int_t first2 = his3->GetAxis(2)->GetFirst();
531 Int_t last2 = his3->GetAxis(2)->GetLast();
532 //
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
541 //
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));
547 }
548 if (TMath::Abs(x1)<0.06){
549 his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
550 }
551 TH1 * hisDelta = his2->Projection(0);
552 //
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();
558 }
559 Double_t sector = 9.*x2/TMath::Pi();
560 if (sector<0) sector+=18;
561 Double_t dsec = sector-Int_t(sector)-0.5;
562 Double_t z=refX*x1;
563 (*pcstream)<<hname<<
564 "run="<<run<<
565 "bz="<<bz<<
566 "theta="<<x1<<
567 "phi="<<x2<<
568 "z="<<z<< // dummy z
569 "snp="<<x3<<
570 "entries="<<entries<<
571 "mean="<<mean<<
572 "rms="<<rms<<
573 "refX="<<refX<< // track matching refernce plane
574 "type="<<type<< //
575 "sector="<<sector<<
576 "dsec="<<dsec<<
577 "\n";
578 delete hisDelta;
579 //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
580 }
581 delete his2;
582 }
583 delete his3;
584 }
585}
586
587
588
589void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){
590 //
591 // make a distortion map out ou fthe residual histogram
592 // Results are written to the debug streamer - pcstream
593 // Parameters:
594 // his0 - input (4D) residual histogram
595 // pcstream - file to write the tree
596 // run - run number
597 // refX - track matching reference X
598 // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt
599 // marian.ivanov@cern.ch
600 //
601 // Histo axeses
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;
613 //
614 // 1. make default selections
615 //
616 TH1 * hisDelta=0;
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());
624 delete hisDelta;
625 THnSparse *his0= hisInput->Projection(4,idim0);
626 //
627 // 2. Get mean in diferent bins
628 //
629 Int_t nbins1=his0->GetAxis(1)->GetNbins();
630 Int_t first1=his0->GetAxis(1)->GetFirst();
631 Int_t last1 =his0->GetAxis(1)->GetLast();
632 //
633 Double_t bz=AliTrackerBase::GetBz();
634 Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
635 //
636 for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
637 //
638 Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
639 his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
640 //
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();
645 //
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);
649 if (ibin3<first3) {
650 his1->GetAxis(3)->SetRangeUser(-1,1);
651 x3=0;
652 }
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();
657 //
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);
662 //
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();
668 }
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
673 (*pcstream)<<hname<<
674 "run="<<run<<
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
682 "rms="<<rms<<
683 "refX="<<refX<< // track matching refernce plane
684 "type="<<type<< // parameter type
685 "sector="<<sector<< // sector
686 "dsec="<<dsec<< // dummy delta sector
687 "\n";
688 delete hisDelta;
689 printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
690 }
691 delete his3;
692 }
693 delete his1;
694 }
695 delete his0;
696}
697
698
699
700void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){
701 //
702 // make a distortion map out of the residual histogram
703 // Results are written to the debug streamer - pcstream
704 // Parameters:
705 // his0 - input (4D) residual histogram
706 // pcstream - file to write the tree
707 // run - run number
708 // type - 0- y 1-z,2 -snp, 3-theta
709 // bz - magnetic field
710 // marian.ivanov@cern.ch
711
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
716 //3 OBJ: TAxis kY kY
717 //4 OBJ: TAxis kZ kZ
718 //5 OBJ: TAxis is1 is1
719 //6 OBJ: TAxis is0 is0
720 //7. OBJ: TAxis z z
721 //8. OBJ: TAxis IsPrimary IsPrimary
722
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();
727
728 //
729 // Loop over pair of sector:
730 // isPrim - 8 ==> 8
731 // isec0 - 6 ==> 7
732 // isec1 - 5 ==> 6
733 // refX - 2 ==> 5
734 //
735 // phi - 1 ==> 4
736 // z - 7 ==> 3
737 // snp - 3 ==> 2
738 // theta- 4 ==> 1
739 // 0 ==> 0;
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
742 //
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);
746 //
747 //
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();
755 delete hisX;
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);
760 delete hisDelta;
761 //
762 // 1. make default selections
763 //
764 Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
765 THnBase *hisSector1= hisSector0->ProjectionND(5,idim0);
766 //
767 // 2. Get mean in diferent bins
768 //
769 Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4}
770 //
771 // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins();
772 Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
773 Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
774 //
775 for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi
776 //
777 // Phi loop
778 //
779 Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
780 Double_t psec = (9*xPhi/TMath::Pi());
781 if (psec<0) psec+=18;
782 Bool_t isOK0=kFALSE;
783 Bool_t isOK1=kFALSE;
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;
788 //
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));
792 }
793 //
794 htemp = hisSector1->Projection(4);
795 xPhi=htemp->GetMean();
796 delete htemp;
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();
801 //
802 for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z
803 //
804 // Z loop
805 //
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));
809 }
810 htemp = hisPhi->Projection(3);
811 Double_t xZ= htemp->GetMean();
812 delete htemp;
813 THnBase * hisZ= hisPhi->ProjectionND(3,idim);
814 //projected histogram according selection 3 -z
815 //
816 //
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
821 //
822 // Snp loop
823 //
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));
827 }
828 htemp = hisZ->Projection(2);
829 Double_t xSnp= htemp->GetMean();
830 delete htemp;
831 THnBase * hisSnp= hisZ->ProjectionND(2,idim);
832 //projected histogram according selection 2 - snp
833
834 //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
835 Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
836 Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
837 //
838 for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
839
840
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));
844 }
845 htemp = hisSnp->Projection(1);
846 Double_t xTheta=htemp->GetMean();
847 delete htemp;
848 hisDelta = hisSnp->Projection(0);
849 //
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();
855 }
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
860 (*pcstream)<<hname<<
861 "run="<<run<<
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
869 //
870 "theta="<<xTheta<< // theta
871 "phi="<<xPhi<< // phi (alpha)
872 "z="<<xZ<< // z
873 "snp="<<xSnp<< // snp
874
875 //
876 "entries="<<entries<< // entries in bin
877 "mean="<<mean<< // mean
878 "rms="<<rms<< // rms
879 "refX="<<refX<< // track matching reference plane
880 "\n";
881 delete hisDelta;
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);
883 //
884 }//ibinTheta
885 delete hisSnp;
886 } //ibinSnp
887 delete hisZ;
888 }//ibinZ
889 delete hisPhi;
890 }//ibinPhi
891 delete hisSector1;
892 }//isec1
893 delete hisSector0;
894 }//isec0
895}
896
897
898
899
900
901// void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){
902// //
903// // Make a laser fit tree for global minimization
904// //
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
909
910// // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
911// //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
912// //
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.",&ltr);
929// Int_t entries= tree->GetEntries();
930// TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root");
931// Double_t bz=AliTrackerBase::GetBz();
932// //
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();
939// }
940// //
941// TVectorD fit10(5);
942// TVectorD fit5(5);
943// printf("Entry\t%d\n",ientry);
944// for (Int_t irow0=0; irow0<158; irow0+=1){
945// //
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;
951
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);
972// }
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);
979
980// if (fitter10.GetNpoints()>10){
981// fitter10.Eval();
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);
987// //
988// // make short range fit
989// //
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;
1005// //
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);
1009// }
1010// }else{
1011// isOK=kFALSE;
1012// }
1013// if (fitter5.GetNpoints()<8) isOK=kFALSE;
1014
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);
1018// if (isOK){
1019// fitter5.Eval();
1020// rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4));
1021// offset5 = fitter5.GetParameter(0);
1022// slope5 = fitter5.GetParameter(0);
1023// }
1024// //
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];
1035// //
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;
1046// //
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<<
1064// "rms5="<<rms5<<
1065// "fit10.="<<&fit10<<
1066// "fit5.="<<&fit5<<
1067// "measY="<<measY<<
1068// "mean="<<mean<<
1069// "idealY="<<idealY<<
1070// "corr="<<corr<<
1071// "isOK="<<isOK<<
1072// "\n";
1073// }
1074// }
1075// delete pcstream;
1076// }
1077
1078
1079void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1080 //
1081 // Make a fit tree:
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
1087 //
1088 //
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
1096
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();
1102 //
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;
1107 Float_t refX;
1108 Int_t run;
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",&sector);
1117 tinput->SetBranchAddress("dsec",&dsec);
1118 tinput->SetBranchAddress("refX",&refX);
1119 TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1120 //
1121 Int_t nentries=tinput->GetEntries();
1122 Int_t ncorr=corrArray->GetEntries();
1123 Double_t corrections[100]={0}; //
1124 Double_t tPar[5];
1125 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1126 Int_t dir=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;}
1133 //
1134 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1135 tinput->GetEntry(ientry);
1136 if (TMath::Abs(snp)>kMaxSnp) continue;
1137 tPar[0]=0;
1138 tPar[1]=theta*refX;
1139 if (dtype==2) tPar[1]=theta*kRadiusTPC;
1140 tPar[2]=snp;
1141 tPar[3]=theta;
1142 tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1143 if (dtype==4){
1144 // tracks crossing CE
1145 tPar[1]=0; // track at the CE
1146 //if (TMath::Abs(theta) <0.05) continue; // deep cross
1147 }
1148
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);
1155
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
1159 }
1160 }
1161 //
1162 tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1163 AliExternalTrackParam track(refX,phi,tPar,cov);
1164 Double_t xyz[3];
1165 track.GetXYZ(xyz);
1166 Int_t id=0;
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
1181 "sector="<<sector<<
1182 "dsec="<<dsec<<
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
1188 "pt="<<pt<< // pt
1189 "id="<<id<< // track id
1190 "entries="<<entries;// number of entries in bin
1191 //
1192 Bool_t isOK=kTRUE;
1193 if (entries<kMinEntries) isOK=kFALSE;
1194 //
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;}
1207 //
1208 if (trackOut){
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());
1213 //
1214 corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1215 delete trackOut;
1216 }else{
1217 corrections[icorr]=0;
1218 isOK=kFALSE;
1219 }
1220 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
1221 }
1222 (*pcstream)<<"fit"<<
1223 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1224 }
1225
1226 if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1227 //
1228 // special case of the TPC tracks crossing the CE
1229 //
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;
1237 //
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);
1242 //
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;
1248 //
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;
1254 //
1255 corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1256 corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1257 if (isOK)
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)
1264 ){
1265 isOK=kFALSE;
1266 }
1267 delete trackOut0;
1268 delete trackOut1;
1269 }else{
1270 corrections[icorr]=0;
1271 isOK=kFALSE;
1272 }
1273 //
1274 //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1275 }
1276 (*pcstream)<<"fit"<<
1277 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1278 }
1279 //
1280 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1281 }
1282
1283
1284 delete pcstream;
1285}
1286
1287
1288
1289void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){
1290 //
1291 // Make a fit tree:
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
1297 //
1298 //
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
1306
1307 const Double_t kMaxSnp = 0.8;
1308 const Int_t kMinEntries=200;
1309 // AliTPCROC *tpcRoc =AliTPCROC::Instance();
1310 //
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;
1314 Int_t isec1, isec0;
1315 Double_t refXD;
1316 Float_t refX;
1317 Int_t run;
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",&sector);
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));
1332 //
1333 Int_t nentries=tinput->GetEntries();
1334 Int_t ncorr=corrArray->GetEntries();
1335 Double_t corrections[100]={0}; //
1336 Double_t tPar[5];
1337 Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1338 Int_t dir=0;
1339 //
1340 for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1341 tinput->GetEntry(ientry);
1342 refX=refXD;
1343 Int_t id=-1;
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;
1347 //
1348 dir=-1;
1349 tPar[0]=0;
1350 tPar[1]=globalZ;
1351 tPar[2]=snp;
1352 tPar[3]=theta;
1353 tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
1354 Double_t pt=1./tPar[4];
1355 //
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];
1360 track.GetXYZ(xyz);
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);
1369 //
1370 Bool_t isOK=kTRUE;
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
1374 // Do downscale
1375 if (TMath::Abs(theta)>1) isOK=kFALSE;
1376 //
1377 Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1378 //
1379 (*pcstream)<<"fit"<<
1380 "run="<<run<< //run
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
1389 "sector="<<sector<<
1390 "dsec="<<dsec<<
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
1396 "pt="<<pt<< //pt
1397 "id="<<id<< // track id
1398 "entries="<<entries;// number of entries in bin
1399 //
1400 AliExternalTrackParam *trackOut0 = 0;
1401 AliExternalTrackParam *trackOut1 = 0;
1402 AliExternalTrackParam *ptrackIn0 = 0;
1403 AliExternalTrackParam *ptrackIn1 = 0;
1404
1405 for (Int_t icorr=0; icorr<ncorr; icorr++) {
1406 //
1407 // special case of the TPC tracks crossing the CE
1408 //
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;
1416 //
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);
1421 //
1422 if (trackOut0 && trackOut1){
1423 //
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;
1430 //
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;
1434 //
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
1438 "mean="<<mean<<
1439 "pIn0.="<<ptrackIn0<<
1440 "pIn1.="<<ptrackIn1<<
1441 "pOut0.="<<trackOut0<<
1442 "pOut1.="<<trackOut1<<
1443 "refX="<<refXD<<
1444 "\n";
1445 delete trackOut0;
1446 delete trackOut1;
1447 }else{
1448 corrections[icorr]=0;
1449 isOK=kFALSE;
1450 }
1451 }
1452 (*pcstream)<<"fit"<<
1453 Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1454 }
1455 //
1456 (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1457 }
1458 delete pcstream;
1459}
1460
1461