Line data Source code
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 :
24 : 1. Creation of the distortion maps from the residual histograms
25 :
26 : 2. Making fit trees
27 :
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
32 : dO_{i}/dp_{i}
33 :
34 : 4. Global fit - later
35 : d0 = sum{ki*dO_{i}/dp_{i}} - linear fitting of the amplitudes ki
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"
88 : #include "AliLog.h"
89 : #include "AliRieman.h"
90 :
91 6 : ClassImp(AliTPCCorrectionFit)
92 :
93 : Double_t AliTPCCorrectionFit::fgMaxChi2HelixAt=1;
94 :
95 : AliTPCCorrectionFit::AliTPCCorrectionFit():
96 0 : TNamed("TPCCorrectionFit","TPCCorrectionFit")
97 0 : {
98 : //
99 : // default constructor
100 : //
101 0 : }
102 :
103 0 : AliTPCCorrectionFit::~AliTPCCorrectionFit() {
104 : //
105 : // Destructor
106 : //
107 0 : }
108 :
109 : /*
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){
111 : //
112 : // Evaluation distortion at point using the linear approximation
113 : //
114 : // refX - reference X where phi and snp is calculated
115 : // evalX - ref X where the distortion is evaluated
116 :
117 : if (wt<50){
118 : AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
119 : pcorr->SetOmegaTauT1T2(wt,t1,t2);
120 : }
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.);
127 : return 0;
128 : }
129 : */
130 :
131 :
132 :
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){
134 : //
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)
137 : //
138 : // Parameters:
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
148 : // //3
149 : // 4 - distortion in the curvature
150 : // nsteps - number of fit points
151 : //
152 : // return value - distortion at point evalX with type ptype
153 : //
154 0 : static TLinearFitter fitter(3,"pol2");
155 0 : fitter.ClearPoints();
156 0 : if (nsteps<3) nsteps=3;
157 0 : AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
158 0 : if (pcorr==0) return 0;
159 0 : if (wt<50){
160 0 : if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
161 : }
162 0 : Double_t deltaX=(245-85)/(nsteps);
163 0 : Double_t curv=2.*snp/refX;
164 0 : Double_t dPhi0=TMath::ASin(snp);
165 :
166 0 : for (Int_t istep=0; istep<(nsteps+1); istep++){
167 : //
168 0 : Double_t localX =85.+deltaX*istep;
169 0 : Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
170 0 : Double_t localPhi=phi0+dPhi;
171 0 : Double_t sector = 9*localPhi/TMath::Pi();
172 0 : if (sector<0) sector+=18;
173 0 : if (sector>18) sector-=18;
174 0 : Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
175 0 : Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
176 0 : Double_t x[1]={localX-dlocalX};
177 0 : fitter.AddPoint(x,dy);
178 0 : }
179 0 : fitter.Eval();
180 : Double_t par[3];
181 0 : par[0]=fitter.GetParameter(0);
182 0 : par[1]=fitter.GetParameter(1);
183 0 : par[2]=fitter.GetParameter(2);
184 :
185 0 : Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
186 0 : if (schi2> fgMaxChi2HelixAt){
187 0 : ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("%s\tbad chi2:\t%2.2f",pcorr->GetName(),schi2).Data());
188 0 : if (pstatus) (*pstatus)=kFALSE;
189 0 : return 0;
190 : }
191 0 : if (pstatus) (*pstatus)=kTRUE;
192 0 : if (ptype==0) return par[0]+par[1]*evalX+par[2]*evalX*evalX;
193 0 : if (ptype==2) return par[1]+2*par[2]*evalX;
194 0 : if (ptype==4) return par[2];
195 0 : if (ptype==5) return schi2;
196 0 : return 0;
197 0 : }
198 :
199 :
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){
201 : //
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)
205 : //
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 : // *
212 : // Parameters:
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
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 : //
228 0 : static TLinearFitter fitter(3,"pol2");
229 0 : fitter.ClearPoints();
230 : //
231 0 : if (nsteps<4) nsteps=4;
232 0 : Double_t deltaX=(245-85)/(nsteps);
233 0 : AliRieman rieman(nsteps);
234 0 : if (wt<50){
235 0 : AliTPCCorrection* pcorr = AliTPCCorrection::GetVisualCorrection(corr);
236 0 : if (pcorr) pcorr->SetOmegaTauT1T2(wt,t1,t2);
237 0 : }
238 0 : Double_t curv=2.*snp/refX;
239 0 : Double_t dPhi0=TMath::ASin(snp);
240 0 : for (Int_t istep=0; istep<(nsteps+1); istep++){
241 : //
242 0 : Double_t localX =85.+deltaX*istep;
243 0 : Double_t dPhi =TMath::ASin(0.5*localX*curv)-dPhi0;
244 0 : Double_t localPhi=phi0+dPhi;
245 0 : Double_t sector = 9*localPhi/TMath::Pi();
246 0 : if (sector<0) sector+=18;
247 0 : if (sector>18) sector-=18;
248 0 : Double_t dy=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr);
249 0 : Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr);
250 0 : Double_t x[1]={localX-dlocalX};
251 0 : Double_t z=theta*x[0];
252 0 : rieman.AddPoint(x[0],dy,z,0.1,0.1);
253 0 : fitter.AddPoint(x,dy);
254 0 : }
255 0 : fitter.Eval();
256 0 : rieman.Update();
257 : //
258 0 : Double_t schi2= TMath::Sqrt(fitter.GetChisquare()/nsteps);
259 0 : if (schi2>fgMaxChi2HelixAt){
260 0 : ::Error("AliTPCCorrectionFit::EvalAtHelix",TString::Format("Bad chi2\t%2.2f",schi2).Data());
261 0 : if (pstatus) (*pstatus)=kFALSE;
262 0 : return 0;
263 : }
264 0 : if (pstatus) (*pstatus)=kTRUE;
265 0 : if (ptype==0) return rieman.GetYat(evalX);
266 0 : if (ptype==2) return rieman.GetDYat(evalX);
267 0 : if (ptype==4) return rieman.GetC();
268 0 : return 0;
269 0 : }
270 :
271 :
272 :
273 :
274 : void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){
275 : //
276 : // Create cluster distortion map
277 : //
278 0 : TFile *falign = TFile::Open("CalibObjects.root");
279 0 : TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign");
280 0 : if (!arrayAlign) {
281 0 : AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
282 0 : return;
283 : }
284 0 : AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC");
285 0 : if (!align) {
286 0 : AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task");
287 0 : return;
288 : }
289 0 : TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root");
290 :
291 0 : THnBase * hdY = (THnBase*)align->GetClusterDelta(0);
292 : //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1);
293 0 : AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz);
294 : // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz);
295 :
296 0 : const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"};
297 0 : for (Int_t ihis=0; ihis<4; ihis++){
298 0 : THnSparse * hisAlign =align->GetTrackletDelta(ihis);
299 0 : AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz);
300 : }
301 0 : delete pcstream;
302 0 : delete falign;
303 0 : }
304 :
305 :
306 :
307 :
308 :
309 : void 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 0 : Int_t nbins1=hisInput->GetAxis(1)->GetNbins();
323 0 : Int_t nbins2=hisInput->GetAxis(2)->GetNbins();
324 0 : Int_t nbins3=hisInput->GetAxis(3)->GetNbins();
325 : TF1 *fgaus=0;
326 : TH3F * hisResMap3D =
327 0 : new TH3F("his3D","his3D",
328 0 : nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
329 0 : nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(),
330 0 : nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax());
331 0 : hisResMap3D->GetXaxis()->SetTitle("sector");
332 0 : hisResMap3D->GetYaxis()->SetTitle("localX");
333 0 : hisResMap3D->GetZaxis()->SetTitle("kZ");
334 :
335 0 : TH2F * hisResMap2D[4] ={0,0,0,0};
336 0 : for (Int_t i=0; i<4; i++){
337 0 : hisResMap2D[i]=
338 0 : new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i),
339 0 : nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(),
340 0 : nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax());
341 0 : hisResMap2D[i]->GetXaxis()->SetTitle("sector");
342 0 : hisResMap2D[i]->GetYaxis()->SetTitle("localX");
343 : }
344 : //
345 : //
346 : //
347 : TF1 * f1= 0;
348 0 : Int_t axis0[4]={0,1,2,3};
349 0 : Int_t axis1[4]={0,1,2,3};
350 : Int_t counter=0;
351 0 : for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){
352 : // phi- sector range
353 0 : hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1);
354 0 : THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0);
355 0 : Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1);
356 : //
357 0 : for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){
358 : // local x range
359 : // kz fits
360 0 : his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1);
361 0 : THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1);
362 0 : Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2);
363 : //
364 : //A side
365 0 : his2->GetAxis(3)->SetRangeUser(0.01,0.3);
366 0 : TH1 * hisA = his2->Projection(0);
367 0 : Double_t meanA= hisA->GetMean();
368 0 : Double_t rmsA= hisA->GetRMS();
369 0 : Double_t entriesA= hisA->GetEntries();
370 0 : delete hisA;
371 : //C side
372 0 : his2->GetAxis(3)->SetRangeUser(0.01,0.3);
373 0 : TH1 * hisC = his2->Projection(0);
374 0 : Double_t meanC= hisC->GetMean();
375 0 : Double_t rmsC= hisC->GetRMS();
376 0 : Double_t entriesC= hisC->GetEntries();
377 0 : delete hisC;
378 0 : his2->GetAxis(3)->SetRangeUser(-1.2,1.2);
379 0 : TH2 * hisAC = his2->Projection(0,3);
380 0 : TProfile *profAC = hisAC->ProfileX();
381 0 : delete hisAC;
382 0 : profAC->Fit("pol1","QNR","QNR",0.05,1);
383 0 : if (!f1) f1=(TF1*)gROOT->FindObject("pol1");
384 0 : Double_t offsetA=f1->GetParameter(0);
385 0 : Double_t slopeA=f1->GetParameter(1);
386 0 : Double_t offsetAE=f1->GetParError(0);
387 0 : Double_t slopeAE=f1->GetParError(1);
388 0 : Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters();
389 0 : profAC->Fit("pol1","QNR","QNR",-1.1,-0.1);
390 0 : f1=(TF1*)gROOT->FindObject("pol1");
391 0 : Double_t offsetC=f1->GetParameter(0);
392 0 : Double_t slopeC=f1->GetParameter(1);
393 0 : Double_t offsetCE=f1->GetParError(0);
394 0 : Double_t slopeCE=f1->GetParError(1);
395 0 : Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters();
396 0 : 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 0 : counter++;
398 0 : (*pcstream)<<"deltaFit"<<
399 0 : "sector="<<sector<<
400 0 : "localX="<<localX<<
401 0 : "meanA="<<meanA<<
402 0 : "rmsA="<<rmsA<<
403 0 : "entriesA="<<entriesA<<
404 0 : "meanC="<<meanC<<
405 0 : "rmsC="<<rmsC<<
406 0 : "entriesC="<<entriesC<<
407 0 : "offsetA="<<offsetA<<
408 0 : "slopeA="<<slopeA<<
409 0 : "offsetAE="<<offsetAE<<
410 0 : "slopeAE="<<slopeAE<<
411 0 : "chi2A="<<chi2A<<
412 0 : "offsetC="<<offsetC<<
413 0 : "slopeC="<<slopeC<<
414 0 : "offsetCE="<<offsetCE<<
415 0 : "slopeCE="<<slopeCE<<
416 0 : "chi2C="<<chi2C<<
417 : "\n";
418 : //
419 0 : hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA);
420 0 : hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA);
421 0 : hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC);
422 0 : hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC);
423 :
424 0 : for (Int_t ibin3=1; ibin3<nbins3; ibin3++){
425 0 : Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3);
426 0 : if (TMath::Abs(kZ)<0.05) continue; // crossing
427 0 : his2->GetAxis(3)->SetRange(ibin3,ibin3);
428 0 : if (TMath::Abs(kZ)>0.15){
429 0 : his2->GetAxis(3)->SetRange(ibin3,ibin3);
430 0 : }
431 0 : TH1 * his = his2->Projection(0);
432 0 : Double_t mean= his->GetMean();
433 0 : Double_t rms= his->GetRMS();
434 0 : 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 0 : hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean);
437 0 : Double_t phi=TMath::Pi()*sector/9;
438 0 : if (phi>TMath::Pi()) phi+=TMath::Pi();
439 0 : Double_t meanG=0;
440 0 : Double_t rmsG=0;
441 0 : if (entries>50){
442 0 : if (!fgaus) {
443 0 : his->Fit("gaus","Q","goff");
444 0 : fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone());
445 0 : }
446 0 : if (fgaus) {
447 0 : his->Fit(fgaus,"Q","goff");
448 0 : meanG=fgaus->GetParameter(1);
449 0 : rmsG=fgaus->GetParameter(2);
450 0 : }
451 : }
452 0 : Double_t dsec=sector-Int_t(sector)-0.5;
453 0 : Double_t snp=dsec*TMath::Pi()/9.;
454 0 : (*pcstream)<<"delta"<<
455 0 : "ptype="<<ptype<<
456 0 : "dtype="<<dtype<<
457 0 : "sector="<<sector<<
458 0 : "dsec="<<dsec<<
459 0 : "snp="<<snp<<
460 0 : "phi="<<phi<<
461 0 : "localX="<<localX<<
462 0 : "kZ="<<kZ<<
463 0 : "theta="<<kZ<<
464 0 : "mean="<<mean<<
465 0 : "rms="<<rms<<
466 0 : "meanG="<<meanG<<
467 0 : "rmsG="<<rmsG<<
468 0 : "entries="<<entries<<
469 0 : "meanA="<<meanA<<
470 0 : "rmsA="<<rmsA<<
471 0 : "entriesA="<<entriesA<<
472 0 : "meanC="<<meanC<<
473 0 : "rmsC="<<rmsC<<
474 0 : "entriesC="<<entriesC<<
475 0 : "offsetA="<<offsetA<<
476 0 : "slopeA="<<slopeA<<
477 0 : "chi2A="<<chi2A<<
478 0 : "offsetC="<<offsetC<<
479 0 : "slopeC="<<slopeC<<
480 0 : "chi2C="<<chi2C<<
481 : "\n";
482 0 : delete his;
483 0 : }
484 0 : delete his2;
485 0 : }
486 0 : delete his1;
487 0 : }
488 0 : hisResMap3D->Write();
489 0 : hisResMap2D[0]->Write();
490 0 : hisResMap2D[1]->Write();
491 0 : hisResMap2D[2]->Write();
492 0 : hisResMap2D[3]->Write();
493 : // delete pcstream;
494 0 : }
495 :
496 :
497 :
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){
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 0 : Int_t idim[4]={0,1,2,3};
517 : //
518 : //
519 : //
520 0 : Int_t nbins3=his0->GetAxis(3)->GetNbins();
521 0 : Int_t first3=his0->GetAxis(3)->GetFirst();
522 0 : Int_t last3 =his0->GetAxis(3)->GetLast();
523 : //
524 0 : for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle
525 0 : his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3));
526 0 : Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
527 0 : THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
528 : //
529 0 : Int_t nbins2 = his3->GetAxis(2)->GetNbins();
530 0 : Int_t first2 = his3->GetAxis(2)->GetFirst();
531 0 : Int_t last2 = his3->GetAxis(2)->GetLast();
532 : //
533 0 : for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
534 0 : his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2));
535 0 : Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
536 0 : THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
537 0 : Int_t nbins1 = his2->GetAxis(1)->GetNbins();
538 0 : Int_t first1 = his2->GetAxis(1)->GetFirst();
539 0 : Int_t last1 = his2->GetAxis(1)->GetLast();
540 0 : for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta
541 : //
542 0 : Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
543 0 : his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
544 0 : if (TMath::Abs(x1)<0.1){
545 0 : if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1));
546 0 : if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1));
547 : }
548 0 : if (TMath::Abs(x1)<0.06){
549 0 : his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1));
550 0 : }
551 0 : TH1 * hisDelta = his2->Projection(0);
552 : //
553 0 : Double_t entries = hisDelta->GetEntries();
554 0 : Double_t mean=0, rms=0;
555 0 : if (entries>kMinEntries){
556 0 : mean = hisDelta->GetMean();
557 0 : rms = hisDelta->GetRMS();
558 0 : }
559 0 : Double_t sector = 9.*x2/TMath::Pi();
560 0 : if (sector<0) sector+=18;
561 0 : Double_t dsec = sector-Int_t(sector)-0.5;
562 0 : Double_t z=refX*x1;
563 0 : (*pcstream)<<hname<<
564 0 : "run="<<run<<
565 0 : "bz="<<bz<<
566 0 : "theta="<<x1<<
567 0 : "phi="<<x2<<
568 0 : "z="<<z<< // dummy z
569 0 : "snp="<<x3<<
570 0 : "entries="<<entries<<
571 0 : "mean="<<mean<<
572 0 : "rms="<<rms<<
573 0 : "refX="<<refX<< // track matching refernce plane
574 0 : "type="<<type<< //
575 0 : "sector="<<sector<<
576 0 : "dsec="<<dsec<<
577 : "\n";
578 0 : delete hisDelta;
579 : //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
580 0 : }
581 0 : delete his2;
582 0 : }
583 0 : delete his3;
584 0 : }
585 0 : }
586 :
587 :
588 :
589 : void 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 0 : Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z
618 0 : hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks
619 0 : hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe
620 0 : hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance
621 0 : hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks"
622 0 : hisDelta= hisInput->Projection(0);
623 0 : hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS());
624 0 : delete hisDelta;
625 0 : THnSparse *his0= hisInput->Projection(4,idim0);
626 : //
627 : // 2. Get mean in diferent bins
628 : //
629 0 : Int_t nbins1=his0->GetAxis(1)->GetNbins();
630 0 : Int_t first1=his0->GetAxis(1)->GetFirst();
631 0 : Int_t last1 =his0->GetAxis(1)->GetLast();
632 : //
633 0 : Double_t bz=AliTrackerBase::GetBz();
634 0 : Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z
635 : //
636 0 : for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta
637 : //
638 0 : Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1);
639 0 : his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1));
640 : //
641 0 : THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1
642 0 : Int_t nbins3 = his1->GetAxis(3)->GetNbins();
643 0 : Int_t first3 = his1->GetAxis(3)->GetFirst();
644 0 : Int_t last3 = his1->GetAxis(3)->GetLast();
645 : //
646 0 : for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex"
647 0 : his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3));
648 0 : Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3);
649 0 : if (ibin3<first3) {
650 0 : his1->GetAxis(3)->SetRangeUser(-1,1);
651 0 : x3=0;
652 0 : }
653 0 : THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3
654 0 : Int_t nbins2 = his3->GetAxis(2)->GetNbins();
655 0 : Int_t first2 = his3->GetAxis(2)->GetFirst();
656 0 : Int_t last2 = his3->GetAxis(2)->GetLast();
657 : //
658 0 : for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){
659 0 : his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2));
660 0 : Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
661 0 : hisDelta = his3->Projection(0);
662 : //
663 0 : Double_t entries = hisDelta->GetEntries();
664 0 : Double_t mean=0, rms=0;
665 0 : if (entries>kMinEntries){
666 0 : mean = hisDelta->GetMean();
667 0 : rms = hisDelta->GetRMS();
668 0 : }
669 0 : Double_t sector = 9.*x2/TMath::Pi();
670 0 : if (sector<0) sector+=18;
671 0 : Double_t dsec = sector-Int_t(sector)-0.5;
672 0 : Double_t snp=0; // dummy snp - equal 0
673 0 : (*pcstream)<<hname<<
674 0 : "run="<<run<<
675 0 : "bz="<<bz<< // magnetic field
676 0 : "theta="<<x1<< // theta
677 0 : "phi="<<x2<< // phi (alpha)
678 0 : "z="<<x3<< // z at "vertex"
679 0 : "snp="<<snp<< // dummy snp
680 0 : "entries="<<entries<< // entries in bin
681 0 : "mean="<<mean<< // mean
682 0 : "rms="<<rms<<
683 0 : "refX="<<refX<< // track matching refernce plane
684 0 : "type="<<type<< // parameter type
685 0 : "sector="<<sector<< // sector
686 0 : "dsec="<<dsec<< // dummy delta sector
687 : "\n";
688 0 : delete hisDelta;
689 0 : printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean);
690 0 : }
691 0 : delete his3;
692 0 : }
693 0 : delete his1;
694 0 : }
695 0 : delete his0;
696 0 : }
697 :
698 :
699 :
700 : void 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 0 : for (Int_t isec0=0; isec0<72; isec0++){
741 0 : 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 0 : hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1);
745 0 : hisSector0=hisInput->Projection(7,index0);
746 : //
747 : //
748 0 : for (Int_t isec1=isec0+1; isec1<72; isec1++){
749 : //if (isec1!=isec0+36) continue;
750 0 : if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue;
751 0 : printf("Sectors %d\t%d\n",isec1,isec0);
752 0 : hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1);
753 0 : TH1 * hisX=hisSector0->Projection(5);
754 0 : Double_t refX= hisX->GetMean();
755 0 : delete hisX;
756 0 : TH1 *hisDelta=hisSector0->Projection(0);
757 0 : Double_t dmean = hisDelta->GetMean();
758 0 : Double_t drms = hisDelta->GetRMS();
759 0 : hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms);
760 0 : delete hisDelta;
761 : //
762 : // 1. make default selections
763 : //
764 0 : Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi }
765 0 : THnBase *hisSector1= hisSector0->ProjectionND(5,idim0);
766 : //
767 : // 2. Get mean in diferent bins
768 : //
769 0 : 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 0 : Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst();
773 0 : Int_t lastPhi =hisSector1->GetAxis(4)->GetLast();
774 : //
775 0 : for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi
776 : //
777 : // Phi loop
778 : //
779 0 : Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi);
780 0 : Double_t psec = (9*xPhi/TMath::Pi());
781 0 : if (psec<0) psec+=18;
782 : Bool_t isOK0=kFALSE;
783 : Bool_t isOK1=kFALSE;
784 0 : if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE;
785 0 : if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE;
786 0 : if (!isOK0) continue;
787 0 : if (!isOK1) continue;
788 : //
789 0 : hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi));
790 0 : if (isec1!=isec0+36) {
791 0 : hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi));
792 0 : }
793 : //
794 0 : htemp = hisSector1->Projection(4);
795 0 : xPhi=htemp->GetMean();
796 0 : delete htemp;
797 0 : THnBase * hisPhi = hisSector1->ProjectionND(4,idim);
798 : //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins();
799 0 : Int_t firstZ = hisPhi->GetAxis(3)->GetFirst();
800 0 : Int_t lastZ = hisPhi->GetAxis(3)->GetLast();
801 : //
802 0 : for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z
803 : //
804 : // Z loop
805 : //
806 0 : hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ));
807 0 : if (isec1!=isec0+36) {
808 0 : hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ));
809 0 : }
810 0 : htemp = hisPhi->Projection(3);
811 0 : Double_t xZ= htemp->GetMean();
812 0 : delete htemp;
813 0 : THnBase * hisZ= hisPhi->ProjectionND(3,idim);
814 : //projected histogram according selection 3 -z
815 : //
816 : //
817 : //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins();
818 0 : Int_t firstSnp = hisZ->GetAxis(2)->GetFirst();
819 0 : Int_t lastSnp = hisZ->GetAxis(2)->GetLast();
820 0 : for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp
821 : //
822 : // Snp loop
823 : //
824 0 : hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp));
825 0 : if (isec1!=isec0+36) {
826 0 : hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp));
827 0 : }
828 0 : htemp = hisZ->Projection(2);
829 0 : Double_t xSnp= htemp->GetMean();
830 0 : delete htemp;
831 0 : THnBase * hisSnp= hisZ->ProjectionND(2,idim);
832 : //projected histogram according selection 2 - snp
833 :
834 : //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins();
835 0 : Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst();
836 0 : Int_t lastTheta = hisSnp->GetAxis(1)->GetLast();
837 : //
838 0 : for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta
839 :
840 :
841 0 : hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta));
842 0 : if (isec1!=isec0+36) {
843 0 : hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta));
844 0 : }
845 0 : htemp = hisSnp->Projection(1);
846 0 : Double_t xTheta=htemp->GetMean();
847 0 : delete htemp;
848 0 : hisDelta = hisSnp->Projection(0);
849 : //
850 0 : Double_t entries = hisDelta->GetEntries();
851 0 : Double_t mean=0, rms=0;
852 0 : if (entries>kMinEntries){
853 0 : mean = hisDelta->GetMean();
854 0 : rms = hisDelta->GetRMS();
855 0 : }
856 0 : Double_t sector = 9.*xPhi/TMath::Pi();
857 0 : if (sector<0) sector+=18;
858 0 : Double_t dsec = sector-Int_t(sector)-0.5;
859 0 : Int_t dtype=1; // TPC alignment type
860 0 : (*pcstream)<<hname<<
861 0 : "run="<<run<<
862 0 : "bz="<<bz<< // magnetic field
863 0 : "ptype="<<type<< // parameter type
864 0 : "dtype="<<dtype<< // parameter type
865 0 : "isec0="<<isec0<< // sector 0
866 0 : "isec1="<<isec1<< // sector 1
867 0 : "sector="<<sector<< // sector as float
868 0 : "dsec="<<dsec<< // delta sector
869 : //
870 0 : "theta="<<xTheta<< // theta
871 0 : "phi="<<xPhi<< // phi (alpha)
872 0 : "z="<<xZ<< // z
873 0 : "snp="<<xSnp<< // snp
874 :
875 : //
876 0 : "entries="<<entries<< // entries in bin
877 0 : "mean="<<mean<< // mean
878 0 : "rms="<<rms<< // rms
879 0 : "refX="<<refX<< // track matching reference plane
880 : "\n";
881 0 : 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 0 : }//ibinTheta
885 0 : delete hisSnp;
886 0 : } //ibinSnp
887 0 : delete hisZ;
888 0 : }//ibinZ
889 0 : delete hisPhi;
890 0 : }//ibinPhi
891 0 : delete hisSector1;
892 0 : }//isec1
893 0 : delete hisSector0;
894 0 : }//isec0
895 0 : }
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.",<r);
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 :
1079 : void 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 0 : const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1104 : // const Double_t kB2C=-0.299792458e-3;
1105 : const Int_t kMinEntries=20;
1106 0 : Double_t phi,theta, snp, mean,rms, entries,sector,dsec;
1107 0 : Float_t refX;
1108 0 : Int_t run;
1109 0 : tinput->SetBranchAddress("run",&run);
1110 0 : tinput->SetBranchAddress("theta",&theta);
1111 0 : tinput->SetBranchAddress("phi", &phi);
1112 0 : tinput->SetBranchAddress("snp",&snp);
1113 0 : tinput->SetBranchAddress("mean",&mean);
1114 0 : tinput->SetBranchAddress("rms",&rms);
1115 0 : tinput->SetBranchAddress("entries",&entries);
1116 0 : tinput->SetBranchAddress("sector",§or);
1117 0 : tinput->SetBranchAddress("dsec",&dsec);
1118 0 : tinput->SetBranchAddress("refX",&refX);
1119 0 : TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset));
1120 : //
1121 0 : Int_t nentries=tinput->GetEntries();
1122 0 : Int_t ncorr=corrArray->GetEntries();
1123 0 : Double_t corrections[100]={0}; //
1124 0 : Double_t tPar[5];
1125 0 : Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1126 : Int_t dir=0;
1127 0 : if (dtype==5 || dtype==6) dtype=4;
1128 0 : if (dtype==0) { dir=-1;}
1129 0 : if (dtype==1) { dir=1;}
1130 0 : if (dtype==2) { dir=-1;}
1131 0 : if (dtype==3) { dir=1;}
1132 0 : if (dtype==4) { dir=-1;}
1133 : //
1134 0 : for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1135 0 : tinput->GetEntry(ientry);
1136 0 : if (TMath::Abs(snp)>kMaxSnp) continue;
1137 0 : tPar[0]=0;
1138 0 : tPar[1]=theta*refX;
1139 0 : if (dtype==2) tPar[1]=theta*kRadiusTPC;
1140 0 : tPar[2]=snp;
1141 0 : tPar[3]=theta;
1142 0 : tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0
1143 0 : if (dtype==4){
1144 : // tracks crossing CE
1145 0 : tPar[1]=0; // track at the CE
1146 : //if (TMath::Abs(theta) <0.05) continue; // deep cross
1147 0 : }
1148 :
1149 0 : if (TMath::Abs(snp) >kcutSnp) continue;
1150 0 : if (TMath::Abs(theta) >kcutTheta) continue;
1151 0 : printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1152 0 : Double_t bz=AliTrackerBase::GetBz();
1153 0 : if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks
1154 0 : if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2);
1155 :
1156 0 : if (dtype==2 && TMath::Abs(bz)>0.1 ) {
1157 0 : tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);//
1158 : // snp at the TPC inner radius in case the vertex match used
1159 0 : }
1160 : }
1161 : //
1162 0 : tPar[4]+=(gRandom->Rndm()-0.5)*0.02;
1163 0 : AliExternalTrackParam track(refX,phi,tPar,cov);
1164 0 : Double_t xyz[3];
1165 0 : track.GetXYZ(xyz);
1166 0 : Int_t id=0;
1167 0 : Double_t pt=1./tPar[4];
1168 0 : 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 0 : Double_t refXD=refX;
1171 0 : (*pcstream)<<"fit"<<
1172 0 : "run="<<run<< // run number
1173 0 : "bz="<<bz<< // magnetic filed used
1174 0 : "dtype="<<dtype<< // detector match type
1175 0 : "ptype="<<ptype<< // parameter type
1176 0 : "theta="<<theta<< // theta
1177 0 : "phi="<<phi<< // phi
1178 0 : "snp="<<snp<< // snp
1179 0 : "mean="<<mean<< // mean dist value
1180 0 : "rms="<<rms<< // rms
1181 0 : "sector="<<sector<<
1182 0 : "dsec="<<dsec<<
1183 0 : "refX="<<refXD<< // referece X as double
1184 0 : "gx="<<xyz[0]<< // global position at reference
1185 0 : "gy="<<xyz[1]<< // global position at reference
1186 0 : "gz="<<xyz[2]<< // global position at reference
1187 0 : "dRrec="<<dRrec<< // delta Radius in reconstruction
1188 0 : "pt="<<pt<< // pt
1189 0 : "id="<<id<< // track id
1190 0 : "entries="<<entries;// number of entries in bin
1191 : //
1192 0 : Bool_t isOK=kTRUE;
1193 0 : if (entries<kMinEntries) isOK=kFALSE;
1194 : //
1195 0 : if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1196 0 : AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1197 0 : corrections[icorr]=0;
1198 0 : if (entries>kMinEntries){
1199 0 : AliExternalTrackParam trackIn(refX,phi,tPar,cov);
1200 : AliExternalTrackParam *trackOut = 0;
1201 0 : if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream);
1202 0 : if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0);
1203 0 : if (dtype==0) {dir= -1;}
1204 0 : if (dtype==1) {dir= 1;}
1205 0 : if (dtype==2) {dir= -1;}
1206 0 : if (dtype==3) {dir= 1;}
1207 : //
1208 0 : if (trackOut){
1209 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1210 0 : if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE;
1211 0 : if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1212 : // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz());
1213 : //
1214 0 : corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype];
1215 0 : delete trackOut;
1216 : }else{
1217 0 : corrections[icorr]=0;
1218 0 : isOK=kFALSE;
1219 : }
1220 : //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out
1221 0 : }
1222 0 : (*pcstream)<<"fit"<<
1223 0 : Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1224 0 : }
1225 :
1226 0 : if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) {
1227 : //
1228 : // special case of the TPC tracks crossing the CE
1229 : //
1230 0 : AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1231 0 : corrections[icorr]=0;
1232 0 : if (entries>kMinEntries){
1233 0 : AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex
1234 0 : AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet
1235 : AliExternalTrackParam *trackOut0 = 0;
1236 : AliExternalTrackParam *trackOut1 = 0;
1237 : //
1238 0 : if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1239 0 : if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1240 0 : if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1241 0 : if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1242 : //
1243 0 : if (trackOut0 && trackOut1){
1244 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1245 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1246 0 : if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1247 0 : if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1248 : //
1249 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE;
1250 0 : if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1251 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1252 0 : if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE;
1253 0 : if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE;
1254 : //
1255 0 : corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1256 0 : corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1257 0 : if (isOK)
1258 0 : if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)||
1259 0 : (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)||
1260 0 : (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)||
1261 0 : (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)||
1262 0 : (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)||
1263 0 : (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001)
1264 : ){
1265 0 : isOK=kFALSE;
1266 0 : }
1267 0 : delete trackOut0;
1268 0 : delete trackOut1;
1269 : }else{
1270 0 : corrections[icorr]=0;
1271 0 : isOK=kFALSE;
1272 : }
1273 : //
1274 : //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup
1275 0 : }
1276 0 : (*pcstream)<<"fit"<<
1277 0 : Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1278 0 : }
1279 : //
1280 0 : (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1281 0 : }
1282 :
1283 :
1284 0 : delete pcstream;
1285 0 : }
1286 :
1287 :
1288 :
1289 : void 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 0 : const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1312 : // const Double_t kB2C=-0.299792458e-3;
1313 0 : Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ;
1314 0 : Int_t isec1, isec0;
1315 0 : Double_t refXD;
1316 : Float_t refX;
1317 0 : Int_t run;
1318 0 : tinput->SetBranchAddress("run",&run);
1319 0 : tinput->SetBranchAddress("theta",&theta);
1320 0 : tinput->SetBranchAddress("phi", &phi);
1321 0 : tinput->SetBranchAddress("snp",&snp);
1322 0 : tinput->SetBranchAddress("mean",&mean);
1323 0 : tinput->SetBranchAddress("rms",&rms);
1324 0 : tinput->SetBranchAddress("entries",&entries);
1325 0 : tinput->SetBranchAddress("sector",§or);
1326 0 : tinput->SetBranchAddress("dsec",&dsec);
1327 0 : tinput->SetBranchAddress("refX",&refXD);
1328 0 : tinput->SetBranchAddress("z",&globalZ);
1329 0 : tinput->SetBranchAddress("isec0",&isec0);
1330 0 : tinput->SetBranchAddress("isec1",&isec1);
1331 0 : TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset));
1332 : //
1333 0 : Int_t nentries=tinput->GetEntries();
1334 0 : Int_t ncorr=corrArray->GetEntries();
1335 0 : Double_t corrections[100]={0}; //
1336 0 : Double_t tPar[5];
1337 0 : 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 0 : for (Int_t ientry=offset; ientry<nentries; ientry+=step){
1341 0 : tinput->GetEntry(ientry);
1342 0 : refX=refXD;
1343 0 : Int_t id=-1;
1344 0 : if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side
1345 0 : if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side
1346 0 : if (dtype==10 && id==-1) continue;
1347 : //
1348 : dir=-1;
1349 0 : tPar[0]=0;
1350 0 : tPar[1]=globalZ;
1351 0 : tPar[2]=snp;
1352 0 : tPar[3]=theta;
1353 0 : tPar[4]=(gRandom->Rndm()-0.1)*0.2; //
1354 0 : Double_t pt=1./tPar[4];
1355 : //
1356 0 : printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms);
1357 0 : Double_t bz=AliTrackerBase::GetBz();
1358 0 : AliExternalTrackParam track(refX,phi,tPar,cov);
1359 0 : Double_t xyz[3],xyzIn[3],xyzOut[3];
1360 0 : track.GetXYZ(xyz);
1361 0 : track.GetXYZAt(85,bz,xyzIn);
1362 0 : track.GetXYZAt(245,bz,xyzOut);
1363 0 : Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]);
1364 0 : Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]);
1365 0 : Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]);
1366 0 : Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5);
1367 0 : Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5);
1368 0 : Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5);
1369 : //
1370 0 : Bool_t isOK=kTRUE;
1371 0 : if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector
1372 0 : if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector
1373 0 : if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin
1374 : // Do downscale
1375 0 : if (TMath::Abs(theta)>1) isOK=kFALSE;
1376 : //
1377 0 : Double_t dRrec=0; // dummy value - needed for points - e.g for laser
1378 : //
1379 0 : (*pcstream)<<"fit"<<
1380 0 : "run="<<run<< //run
1381 0 : "bz="<<bz<< // magnetic filed used
1382 0 : "dtype="<<dtype<< // detector match type
1383 0 : "ptype="<<ptype<< // parameter type
1384 0 : "theta="<<theta<< // theta
1385 0 : "phi="<<phi<< // phi
1386 0 : "snp="<<snp<< // snp
1387 0 : "mean="<<mean<< // mean dist value
1388 0 : "rms="<<rms<< // rms
1389 0 : "sector="<<sector<<
1390 0 : "dsec="<<dsec<<
1391 0 : "refX="<<refXD<< // referece X
1392 0 : "gx="<<xyz[0]<< // global position at reference
1393 0 : "gy="<<xyz[1]<< // global position at reference
1394 0 : "gz="<<xyz[2]<< // global position at reference
1395 0 : "dRrec="<<dRrec<< // delta Radius in reconstruction
1396 0 : "pt="<<pt<< //pt
1397 0 : "id="<<id<< // track id
1398 0 : "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 0 : for (Int_t icorr=0; icorr<ncorr; icorr++) {
1406 : //
1407 : // special case of the TPC tracks crossing the CE
1408 : //
1409 0 : AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr);
1410 0 : corrections[icorr]=0;
1411 0 : if (entries>kMinEntries &&isOK){
1412 0 : AliExternalTrackParam trackIn0(refX,phi,tPar,cov);
1413 0 : AliExternalTrackParam trackIn1(refX,phi,tPar,cov);
1414 : ptrackIn1=&trackIn0;
1415 : ptrackIn0=&trackIn1;
1416 : //
1417 0 : if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream);
1418 0 : if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0);
1419 0 : if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream);
1420 0 : if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0);
1421 : //
1422 0 : if (trackOut0 && trackOut1){
1423 : //
1424 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE;
1425 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1426 : // rotate all tracks to the same frame
1427 0 : if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1428 0 : if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1429 0 : if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE;
1430 : //
1431 0 : if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1432 0 : if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1433 0 : if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE;
1434 : //
1435 0 : corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]);
1436 0 : corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]);
1437 0 : (*pcstream)<<"fitDebug"<< // just to debug the correction
1438 0 : "mean="<<mean<<
1439 0 : "pIn0.="<<ptrackIn0<<
1440 0 : "pIn1.="<<ptrackIn1<<
1441 0 : "pOut0.="<<trackOut0<<
1442 0 : "pOut1.="<<trackOut1<<
1443 0 : "refX="<<refXD<<
1444 : "\n";
1445 0 : delete trackOut0;
1446 0 : delete trackOut1;
1447 : }else{
1448 0 : corrections[icorr]=0;
1449 0 : isOK=kFALSE;
1450 : }
1451 0 : }
1452 0 : (*pcstream)<<"fit"<<
1453 0 : Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value
1454 : }
1455 : //
1456 0 : (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1457 0 : }
1458 0 : delete pcstream;
1459 0 : }
1460 :
1461 :
|
