]>
Commit | Line | Data |
---|---|---|
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 | ||
24 | 1. Creation of the distortion maps from the residual histograms | |
25 | 2. Making fit trees | |
26 | 3. Parameters to calculate- dO/dp | |
27 | ||
28 | Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be | |
29 | extracted from the old macros | |
30 | ||
31 | ||
32 | */ | |
33 | ||
34 | #include "Riostream.h" | |
35 | #include <fstream> | |
36 | #include "TMap.h" | |
37 | #include "TGraphErrors.h" | |
38 | #include "AliExternalTrackParam.h" | |
39 | #include "TFile.h" | |
40 | #include "TGraph.h" | |
41 | #include "TMultiGraph.h" | |
42 | #include "TCanvas.h" | |
43 | #include "THnSparse.h" | |
44 | #include "THnBase.h" | |
45 | #include "TProfile.h" | |
46 | #include "TROOT.h" | |
47 | #include "TLegend.h" | |
48 | #include "TPad.h" | |
49 | #include "TH2D.h" | |
50 | #include "TH3D.h" | |
51 | #include "AliTPCROC.h" | |
52 | #include "AliTPCCalROC.h" | |
53 | #include "AliESDfriend.h" | |
54 | #include "AliTPCcalibTime.h" | |
55 | #include "AliSplineFit.h" | |
56 | #include "AliCDBMetaData.h" | |
57 | #include "AliCDBId.h" | |
58 | #include "AliCDBManager.h" | |
59 | #include "AliCDBStorage.h" | |
60 | #include "AliTPCcalibBase.h" | |
61 | #include "AliTPCcalibDB.h" | |
62 | #include "AliTPCcalibDButil.h" | |
63 | #include "AliRelAlignerKalman.h" | |
64 | #include "AliTPCParamSR.h" | |
65 | #include "AliTPCcalibTimeGain.h" | |
66 | #include "AliTPCcalibGainMult.h" | |
67 | #include "AliSplineFit.h" | |
68 | #include "AliTPCComposedCorrection.h" | |
69 | #include "AliTPCExBTwist.h" | |
70 | #include "AliTPCCalibGlobalMisalignment.h" | |
71 | #include "TStatToolkit.h" | |
72 | #include "TChain.h" | |
73 | #include "TCut.h" | |
74 | #include "AliTrackerBase.h" | |
75 | #include "AliTPCCorrectionFit.h" | |
76 | #include "AliTPCLaserTrack.h" | |
77 | #include "TDatabasePDG.h" | |
78 | #include "AliTPCcalibAlign.h" | |
4425a337 | 79 | #include "AliLog.h" |
93096a07 | 80 | |
81 | ClassImp(AliTPCCorrectionFit) | |
82 | ||
83 | AliTPCCorrectionFit::AliTPCCorrectionFit(): | |
84 | TNamed("TPCCorrectionFit","TPCCorrectionFit") | |
85 | { | |
86 | // | |
87 | // default constructor | |
88 | // | |
89 | } | |
90 | ||
91 | AliTPCCorrectionFit::~AliTPCCorrectionFit() { | |
92 | // | |
93 | // Destructor | |
94 | // | |
95 | } | |
96 | ||
97 | ||
98 | Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){ | |
99 | // | |
100 | // | |
101 | // | |
102 | Double_t sector = 9*phi/TMath::Pi(); | |
103 | if (sector<0) sector+=18; | |
104 | Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr); | |
105 | Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr); | |
106 | if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.); | |
107 | if (ptype==2) return (y245-y85)/(245.-85.); | |
108 | return 0; | |
109 | } | |
110 | ||
111 | ||
112 | ||
113 | ||
114 | Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){ | |
115 | // | |
116 | // Fit the distortion along the line with the parabolic model | |
117 | // Parameters: | |
118 | // phi0 - phi at the entrance of the TPC | |
119 | // snp - local inclination angle at the entrance of the TPC | |
120 | // refX - ref X where the distortion is evanluated | |
121 | // theta | |
122 | // | |
123 | static TLinearFitter fitter(3,"pol2"); | |
124 | fitter.ClearPoints(); | |
125 | if (nsteps<3) nsteps=3; | |
126 | Double_t deltaX=(245-85)/(nsteps); | |
127 | for (Int_t istep=0; istep<(nsteps+1); istep++){ | |
128 | // | |
129 | Double_t localX =85.+deltaX*istep; | |
130 | Double_t localPhi=phi0+deltaX*snp*istep; | |
131 | Double_t sector = 9*localPhi/TMath::Pi(); | |
132 | if (sector<0) sector+=18; | |
133 | Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr); | |
134 | Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr); | |
135 | Double_t x[1]={localX-dlocalX}; | |
136 | fitter.AddPoint(x,y); | |
137 | } | |
138 | fitter.Eval(); | |
139 | Double_t par[3]; | |
140 | par[0]=fitter.GetParameter(0); | |
141 | par[1]=fitter.GetParameter(1); | |
142 | par[2]=fitter.GetParameter(2); | |
143 | ||
144 | if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX; | |
145 | if (ptype==2) return par[1]+2*par[2]*refX; | |
146 | if (ptype==4) return par[2]; | |
147 | return 0; | |
148 | } | |
149 | ||
150 | ||
151 | ||
152 | ||
153 | void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){ | |
154 | // | |
155 | // Create cluster distortion map | |
156 | // | |
157 | TFile *falign = TFile::Open("CalibObjects.root"); | |
158 | TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign"); | |
63d14e61 | 159 | if (!arrayAlign) { |
4425a337 | 160 | AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task"); |
63d14e61 | 161 | return; |
162 | } | |
93096a07 | 163 | AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC"); |
63d14e61 | 164 | if (!align) { |
4425a337 | 165 | AliWarningGeneral("AliTPCCorrectionFit::CreateAlignMaps","Alignment was not included in the calibration task"); |
63d14e61 | 166 | return; |
167 | } | |
93096a07 | 168 | TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root"); |
169 | ||
170 | THnBase * hdY = (THnBase*)align->GetClusterDelta(0); | |
171 | //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1); | |
172 | AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz); | |
173 | // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz); | |
174 | ||
175 | const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"}; | |
176 | for (Int_t ihis=0; ihis<4; ihis++){ | |
177 | THnSparse * hisAlign =align->GetTrackletDelta(ihis); | |
178 | AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz); | |
179 | } | |
180 | delete pcstream; | |
181 | delete falign; | |
182 | } | |
183 | ||
184 | ||
185 | ||
186 | ||
187 | ||
188 | void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){ | |
189 | // | |
190 | // Make cluster residual map from the n-dimensional histogram | |
191 | // hisInput supposed to have given format: | |
192 | // 4 Dim: | |
193 | // delta, | |
194 | // sector | |
195 | // localX | |
196 | // kZ | |
197 | // Vertex position assumed to be at (0,0,0) | |
198 | // | |
199 | //TTreeSRedirector *pcstream=new TTreeSRedirector(sname); | |
200 | // | |
201 | Int_t nbins1=hisInput->GetAxis(1)->GetNbins(); | |
202 | Int_t nbins2=hisInput->GetAxis(2)->GetNbins(); | |
203 | Int_t nbins3=hisInput->GetAxis(3)->GetNbins(); | |
204 | TF1 *fgaus=0; | |
205 | TH3F * hisResMap3D = | |
206 | new TH3F("his3D","his3D", | |
207 | nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(), | |
208 | nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(), | |
209 | nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax()); | |
210 | hisResMap3D->GetXaxis()->SetTitle("sector"); | |
211 | hisResMap3D->GetYaxis()->SetTitle("localX"); | |
212 | hisResMap3D->GetZaxis()->SetTitle("kZ"); | |
213 | ||
214 | TH2F * hisResMap2D[4] ={0,0,0,0}; | |
215 | for (Int_t i=0; i<4; i++){ | |
216 | hisResMap2D[i]= | |
217 | new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i), | |
218 | nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(), | |
219 | nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax()); | |
220 | hisResMap2D[i]->GetXaxis()->SetTitle("sector"); | |
221 | hisResMap2D[i]->GetYaxis()->SetTitle("localX"); | |
222 | } | |
223 | // | |
224 | // | |
225 | // | |
226 | TF1 * f1= 0; | |
227 | Int_t axis0[4]={0,1,2,3}; | |
228 | Int_t axis1[4]={0,1,2,3}; | |
229 | Int_t counter=0; | |
230 | for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){ | |
231 | // phi- sector range | |
232 | hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1); | |
233 | THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0); | |
234 | Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1); | |
235 | // | |
236 | for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){ | |
237 | // local x range | |
238 | // kz fits | |
239 | his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1); | |
240 | THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1); | |
241 | Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2); | |
242 | // | |
243 | //A side | |
244 | his2->GetAxis(3)->SetRangeUser(0.01,0.3); | |
245 | TH1 * hisA = his2->Projection(0); | |
246 | Double_t meanA= hisA->GetMean(); | |
247 | Double_t rmsA= hisA->GetRMS(); | |
248 | Double_t entriesA= hisA->GetEntries(); | |
249 | delete hisA; | |
250 | //C side | |
251 | his2->GetAxis(3)->SetRangeUser(0.01,0.3); | |
252 | TH1 * hisC = his2->Projection(0); | |
253 | Double_t meanC= hisC->GetMean(); | |
254 | Double_t rmsC= hisC->GetRMS(); | |
255 | Double_t entriesC= hisC->GetEntries(); | |
256 | delete hisC; | |
257 | his2->GetAxis(3)->SetRangeUser(-1.2,1.2); | |
258 | TH2 * hisAC = his2->Projection(0,3); | |
259 | TProfile *profAC = hisAC->ProfileX(); | |
260 | delete hisAC; | |
261 | profAC->Fit("pol1","QNR","QNR",0.05,1); | |
262 | if (!f1) f1=(TF1*)gROOT->FindObject("pol1"); | |
263 | Double_t offsetA=f1->GetParameter(0); | |
264 | Double_t slopeA=f1->GetParameter(1); | |
265 | Double_t offsetAE=f1->GetParError(0); | |
266 | Double_t slopeAE=f1->GetParError(1); | |
267 | Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters(); | |
268 | profAC->Fit("pol1","QNR","QNR",-1.1,-0.1); | |
269 | f1=(TF1*)gROOT->FindObject("pol1"); | |
270 | Double_t offsetC=f1->GetParameter(0); | |
271 | Double_t slopeC=f1->GetParameter(1); | |
272 | Double_t offsetCE=f1->GetParError(0); | |
273 | Double_t slopeCE=f1->GetParError(1); | |
274 | Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters(); | |
275 | if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C); | |
276 | counter++; | |
277 | (*pcstream)<<"deltaFit"<< | |
278 | "sector="<<sector<< | |
279 | "localX="<<localX<< | |
280 | "meanA="<<meanA<< | |
281 | "rmsA="<<rmsA<< | |
282 | "entriesA="<<entriesA<< | |
283 | "meanC="<<meanC<< | |
284 | "rmsC="<<rmsC<< | |
285 | "entriesC="<<entriesC<< | |
286 | "offsetA="<<offsetA<< | |
287 | "slopeA="<<slopeA<< | |
288 | "offsetAE="<<offsetAE<< | |
289 | "slopeAE="<<slopeAE<< | |
290 | "chi2A="<<chi2A<< | |
291 | "offsetC="<<offsetC<< | |
292 | "slopeC="<<slopeC<< | |
293 | "offsetCE="<<offsetCE<< | |
294 | "slopeCE="<<slopeCE<< | |
295 | "chi2C="<<chi2C<< | |
296 | "\n"; | |
297 | // | |
298 | hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA); | |
299 | hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA); | |
300 | hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC); | |
301 | hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC); | |
302 | ||
303 | for (Int_t ibin3=1; ibin3<nbins3; ibin3++){ | |
304 | Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3); | |
305 | if (TMath::Abs(kZ)<0.05) continue; // crossing | |
306 | his2->GetAxis(3)->SetRange(ibin3,ibin3); | |
307 | if (TMath::Abs(kZ)>0.15){ | |
308 | his2->GetAxis(3)->SetRange(ibin3,ibin3); | |
309 | } | |
310 | TH1 * his = his2->Projection(0); | |
311 | Double_t mean= his->GetMean(); | |
312 | Double_t rms= his->GetRMS(); | |
313 | Double_t entries= his->GetEntries(); | |
314 | //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms); | |
315 | hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean); | |
316 | Double_t phi=TMath::Pi()*sector/9; | |
317 | if (phi>TMath::Pi()) phi+=TMath::Pi(); | |
318 | Double_t meanG=0; | |
319 | Double_t rmsG=0; | |
320 | if (entries>50){ | |
321 | if (!fgaus) { | |
322 | his->Fit("gaus","Q","goff"); | |
323 | fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone()); | |
324 | } | |
325 | if (fgaus) { | |
326 | his->Fit(fgaus,"Q","goff"); | |
327 | meanG=fgaus->GetParameter(1); | |
328 | rmsG=fgaus->GetParameter(2); | |
329 | } | |
330 | } | |
331 | Double_t dsec=sector-Int_t(sector)-0.5; | |
332 | Double_t snp=dsec*TMath::Pi()/9.; | |
333 | (*pcstream)<<"delta"<< | |
334 | "ptype="<<ptype<< | |
335 | "dtype="<<dtype<< | |
336 | "sector="<<sector<< | |
337 | "dsec="<<dsec<< | |
338 | "snp="<<snp<< | |
339 | "phi="<<phi<< | |
340 | "localX="<<localX<< | |
341 | "kZ="<<kZ<< | |
342 | "theta="<<kZ<< | |
343 | "mean="<<mean<< | |
344 | "rms="<<rms<< | |
345 | "meanG="<<meanG<< | |
346 | "rmsG="<<rmsG<< | |
347 | "entries="<<entries<< | |
348 | "meanA="<<meanA<< | |
349 | "rmsA="<<rmsA<< | |
350 | "entriesA="<<entriesA<< | |
351 | "meanC="<<meanC<< | |
352 | "rmsC="<<rmsC<< | |
353 | "entriesC="<<entriesC<< | |
354 | "offsetA="<<offsetA<< | |
355 | "slopeA="<<slopeA<< | |
356 | "chi2A="<<chi2A<< | |
357 | "offsetC="<<offsetC<< | |
358 | "slopeC="<<slopeC<< | |
359 | "chi2C="<<chi2C<< | |
360 | "\n"; | |
361 | delete his; | |
362 | } | |
363 | delete his2; | |
364 | } | |
365 | delete his1; | |
366 | } | |
367 | hisResMap3D->Write(); | |
368 | hisResMap2D[0]->Write(); | |
369 | hisResMap2D[1]->Write(); | |
370 | hisResMap2D[2]->Write(); | |
371 | hisResMap2D[3]->Write(); | |
372 | // delete pcstream; | |
373 | } | |
374 | ||
375 | ||
376 | ||
377 | void AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){ | |
378 | // | |
379 | // make a distortion map out ou fthe residual histogram | |
380 | // Results are written to the debug streamer - pcstream | |
381 | // Parameters: | |
382 | // his0 - input (4D) residual histogram | |
383 | // pcstream - file to write the tree | |
384 | // run - run number | |
385 | // refX - track matching reference X | |
386 | // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt | |
387 | // THnSparse axes: | |
388 | // OBJ: TAxis #Delta #Delta | |
389 | // OBJ: TAxis tanTheta tan(#Theta) | |
390 | // OBJ: TAxis phi #phi | |
391 | // OBJ: TAxis snp snp | |
392 | ||
393 | // marian.ivanov@cern.ch | |
394 | const Int_t kMinEntries=10; | |
395 | Int_t idim[4]={0,1,2,3}; | |
396 | // | |
397 | // | |
398 | // | |
399 | Int_t nbins3=his0->GetAxis(3)->GetNbins(); | |
400 | Int_t first3=his0->GetAxis(3)->GetFirst(); | |
401 | Int_t last3 =his0->GetAxis(3)->GetLast(); | |
402 | // | |
403 | for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle | |
404 | his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3)); | |
405 | Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3); | |
406 | THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3 | |
407 | // | |
408 | Int_t nbins2 = his3->GetAxis(2)->GetNbins(); | |
409 | Int_t first2 = his3->GetAxis(2)->GetFirst(); | |
410 | Int_t last2 = his3->GetAxis(2)->GetLast(); | |
411 | // | |
412 | for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi | |
413 | his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2)); | |
414 | Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); | |
415 | THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2 | |
416 | Int_t nbins1 = his2->GetAxis(1)->GetNbins(); | |
417 | Int_t first1 = his2->GetAxis(1)->GetFirst(); | |
418 | Int_t last1 = his2->GetAxis(1)->GetLast(); | |
419 | for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta | |
420 | // | |
421 | Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1); | |
422 | his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1)); | |
423 | if (TMath::Abs(x1)<0.1){ | |
424 | if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1)); | |
425 | if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1)); | |
426 | } | |
427 | if (TMath::Abs(x1)<0.06){ | |
428 | his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1)); | |
429 | } | |
430 | TH1 * hisDelta = his2->Projection(0); | |
431 | // | |
432 | Double_t entries = hisDelta->GetEntries(); | |
433 | Double_t mean=0, rms=0; | |
434 | if (entries>kMinEntries){ | |
435 | mean = hisDelta->GetMean(); | |
436 | rms = hisDelta->GetRMS(); | |
437 | } | |
438 | Double_t sector = 9.*x2/TMath::Pi(); | |
439 | if (sector<0) sector+=18; | |
440 | Double_t dsec = sector-Int_t(sector)-0.5; | |
441 | Double_t z=refX*x1; | |
442 | (*pcstream)<<hname<< | |
443 | "run="<<run<< | |
444 | "bz="<<bz<< | |
445 | "theta="<<x1<< | |
446 | "phi="<<x2<< | |
447 | "z="<<z<< // dummy z | |
448 | "snp="<<x3<< | |
449 | "entries="<<entries<< | |
450 | "mean="<<mean<< | |
451 | "rms="<<rms<< | |
452 | "refX="<<refX<< // track matching refernce plane | |
453 | "type="<<type<< // | |
454 | "sector="<<sector<< | |
455 | "dsec="<<dsec<< | |
456 | "\n"; | |
457 | delete hisDelta; | |
458 | //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean); | |
459 | } | |
460 | delete his2; | |
461 | } | |
462 | delete his3; | |
463 | } | |
464 | } | |
465 | ||
466 | ||
467 | ||
468 | void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){ | |
469 | // | |
470 | // make a distortion map out ou fthe residual histogram | |
471 | // Results are written to the debug streamer - pcstream | |
472 | // Parameters: | |
473 | // his0 - input (4D) residual histogram | |
474 | // pcstream - file to write the tree | |
475 | // run - run number | |
476 | // refX - track matching reference X | |
477 | // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt | |
478 | // marian.ivanov@cern.ch | |
479 | // | |
480 | // Histo axeses | |
481 | // Collection name='TObjArray', class='TObjArray', size=16 | |
482 | // 0. OBJ: TAxis #Delta #Delta | |
483 | // 1. OBJ: TAxis N_{cl} N_{cl} | |
484 | // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm) | |
485 | // 3. OBJ: TAxis z (cm) z (cm) | |
486 | // 4. OBJ: TAxis sin(#phi) sin(#phi) | |
487 | // 5. OBJ: TAxis tan(#theta) tan(#theta) | |
488 | // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV) | |
489 | // 7. OBJ: TAxis pt (GeV) pt (GeV) | |
490 | // 8. OBJ: TAxis alpha alpha | |
491 | const Int_t kMinEntries=10; | |
492 | // | |
493 | // 1. make default selections | |
494 | // | |
495 | TH1 * hisDelta=0; | |
496 | Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z | |
497 | hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks | |
498 | hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe | |
499 | hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance | |
500 | hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks" | |
501 | hisDelta= hisInput->Projection(0); | |
502 | hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS()); | |
503 | delete hisDelta; | |
504 | THnSparse *his0= hisInput->Projection(4,idim0); | |
505 | // | |
506 | // 2. Get mean in diferent bins | |
507 | // | |
508 | Int_t nbins1=his0->GetAxis(1)->GetNbins(); | |
509 | Int_t first1=his0->GetAxis(1)->GetFirst(); | |
510 | Int_t last1 =his0->GetAxis(1)->GetLast(); | |
511 | // | |
512 | Double_t bz=AliTrackerBase::GetBz(); | |
513 | Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z | |
514 | // | |
515 | for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta | |
516 | // | |
517 | Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1); | |
518 | his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1)); | |
519 | // | |
520 | THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1 | |
521 | Int_t nbins3 = his1->GetAxis(3)->GetNbins(); | |
522 | Int_t first3 = his1->GetAxis(3)->GetFirst(); | |
523 | Int_t last3 = his1->GetAxis(3)->GetLast(); | |
524 | // | |
525 | for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex" | |
526 | his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3)); | |
527 | Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3); | |
528 | if (ibin3<first3) { | |
529 | his1->GetAxis(3)->SetRangeUser(-1,1); | |
530 | x3=0; | |
531 | } | |
532 | THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3 | |
533 | Int_t nbins2 = his3->GetAxis(2)->GetNbins(); | |
534 | Int_t first2 = his3->GetAxis(2)->GetFirst(); | |
535 | Int_t last2 = his3->GetAxis(2)->GetLast(); | |
536 | // | |
537 | for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){ | |
538 | his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2)); | |
539 | Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); | |
540 | hisDelta = his3->Projection(0); | |
541 | // | |
542 | Double_t entries = hisDelta->GetEntries(); | |
543 | Double_t mean=0, rms=0; | |
544 | if (entries>kMinEntries){ | |
545 | mean = hisDelta->GetMean(); | |
546 | rms = hisDelta->GetRMS(); | |
547 | } | |
548 | Double_t sector = 9.*x2/TMath::Pi(); | |
549 | if (sector<0) sector+=18; | |
550 | Double_t dsec = sector-Int_t(sector)-0.5; | |
551 | Double_t snp=0; // dummy snp - equal 0 | |
552 | (*pcstream)<<hname<< | |
553 | "run="<<run<< | |
554 | "bz="<<bz<< // magnetic field | |
555 | "theta="<<x1<< // theta | |
556 | "phi="<<x2<< // phi (alpha) | |
557 | "z="<<x3<< // z at "vertex" | |
558 | "snp="<<snp<< // dummy snp | |
559 | "entries="<<entries<< // entries in bin | |
560 | "mean="<<mean<< // mean | |
561 | "rms="<<rms<< | |
562 | "refX="<<refX<< // track matching refernce plane | |
563 | "type="<<type<< // parameter type | |
564 | "sector="<<sector<< // sector | |
565 | "dsec="<<dsec<< // dummy delta sector | |
566 | "\n"; | |
567 | delete hisDelta; | |
568 | printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean); | |
569 | } | |
570 | delete his3; | |
571 | } | |
572 | delete his1; | |
573 | } | |
574 | delete his0; | |
575 | } | |
576 | ||
577 | ||
578 | ||
579 | void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){ | |
580 | // | |
581 | // make a distortion map out of the residual histogram | |
582 | // Results are written to the debug streamer - pcstream | |
583 | // Parameters: | |
584 | // his0 - input (4D) residual histogram | |
585 | // pcstream - file to write the tree | |
586 | // run - run number | |
587 | // type - 0- y 1-z,2 -snp, 3-theta | |
588 | // bz - magnetic field | |
589 | // marian.ivanov@cern.ch | |
590 | ||
591 | //Collection name='TObjArray', class='TObjArray', size=16 | |
592 | //0 OBJ: TAxis delta delta | |
593 | //1 OBJ: TAxis phi phi | |
594 | //2 OBJ: TAxis localX localX | |
595 | //3 OBJ: TAxis kY kY | |
596 | //4 OBJ: TAxis kZ kZ | |
597 | //5 OBJ: TAxis is1 is1 | |
598 | //6 OBJ: TAxis is0 is0 | |
599 | //7. OBJ: TAxis z z | |
600 | //8. OBJ: TAxis IsPrimary IsPrimary | |
601 | ||
602 | const Int_t kMinEntries=10; | |
603 | THnSparse * hisSector0=0; | |
604 | TH1 * htemp=0; // histogram to calculate mean value of parameter | |
605 | // Double_t bz=AliTrackerBase::GetBz(); | |
606 | ||
607 | // | |
608 | // Loop over pair of sector: | |
609 | // isPrim - 8 ==> 8 | |
610 | // isec0 - 6 ==> 7 | |
611 | // isec1 - 5 ==> 6 | |
612 | // refX - 2 ==> 5 | |
613 | // | |
614 | // phi - 1 ==> 4 | |
615 | // z - 7 ==> 3 | |
616 | // snp - 3 ==> 2 | |
617 | // theta- 4 ==> 1 | |
618 | // 0 ==> 0; | |
619 | for (Int_t isec0=0; isec0<72; isec0++){ | |
620 | Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces | |
621 | // | |
622 | //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ? | |
623 | hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1); | |
624 | hisSector0=hisInput->Projection(7,index0); | |
625 | // | |
626 | // | |
627 | for (Int_t isec1=isec0+1; isec1<72; isec1++){ | |
628 | //if (isec1!=isec0+36) continue; | |
629 | if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue; | |
630 | printf("Sectors %d\t%d\n",isec1,isec0); | |
631 | hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1); | |
632 | TH1 * hisX=hisSector0->Projection(5); | |
633 | Double_t refX= hisX->GetMean(); | |
634 | delete hisX; | |
635 | TH1 *hisDelta=hisSector0->Projection(0); | |
636 | Double_t dmean = hisDelta->GetMean(); | |
637 | Double_t drms = hisDelta->GetRMS(); | |
638 | hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms); | |
639 | delete hisDelta; | |
640 | // | |
641 | // 1. make default selections | |
642 | // | |
643 | Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi } | |
644 | THnBase *hisSector1= hisSector0->ProjectionND(5,idim0); | |
645 | // | |
646 | // 2. Get mean in diferent bins | |
647 | // | |
648 | Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4} | |
649 | // | |
650 | // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins(); | |
651 | Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst(); | |
652 | Int_t lastPhi =hisSector1->GetAxis(4)->GetLast(); | |
653 | // | |
654 | for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi | |
655 | // | |
656 | // Phi loop | |
657 | // | |
658 | Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi); | |
659 | Double_t psec = (9*xPhi/TMath::Pi()); | |
660 | if (psec<0) psec+=18; | |
661 | Bool_t isOK0=kFALSE; | |
662 | Bool_t isOK1=kFALSE; | |
663 | if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE; | |
664 | if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE; | |
665 | if (!isOK0) continue; | |
666 | if (!isOK1) continue; | |
667 | // | |
668 | hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi)); | |
669 | if (isec1!=isec0+36) { | |
670 | hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi)); | |
671 | } | |
672 | // | |
673 | htemp = hisSector1->Projection(4); | |
674 | xPhi=htemp->GetMean(); | |
675 | delete htemp; | |
676 | THnBase * hisPhi = hisSector1->ProjectionND(4,idim); | |
677 | //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins(); | |
678 | Int_t firstZ = hisPhi->GetAxis(3)->GetFirst(); | |
679 | Int_t lastZ = hisPhi->GetAxis(3)->GetLast(); | |
680 | // | |
681 | for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z | |
682 | // | |
683 | // Z loop | |
684 | // | |
685 | hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ)); | |
686 | if (isec1!=isec0+36) { | |
687 | hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ)); | |
688 | } | |
689 | htemp = hisPhi->Projection(3); | |
690 | Double_t xZ= htemp->GetMean(); | |
691 | delete htemp; | |
692 | THnBase * hisZ= hisPhi->ProjectionND(3,idim); | |
693 | //projected histogram according selection 3 -z | |
694 | // | |
695 | // | |
696 | //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins(); | |
697 | Int_t firstSnp = hisZ->GetAxis(2)->GetFirst(); | |
698 | Int_t lastSnp = hisZ->GetAxis(2)->GetLast(); | |
699 | for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp | |
700 | // | |
701 | // Snp loop | |
702 | // | |
703 | hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp)); | |
704 | if (isec1!=isec0+36) { | |
705 | hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp)); | |
706 | } | |
707 | htemp = hisZ->Projection(2); | |
708 | Double_t xSnp= htemp->GetMean(); | |
709 | delete htemp; | |
710 | THnBase * hisSnp= hisZ->ProjectionND(2,idim); | |
711 | //projected histogram according selection 2 - snp | |
712 | ||
713 | //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins(); | |
714 | Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst(); | |
715 | Int_t lastTheta = hisSnp->GetAxis(1)->GetLast(); | |
716 | // | |
717 | for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta | |
718 | ||
719 | ||
720 | hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta)); | |
721 | if (isec1!=isec0+36) { | |
722 | hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta)); | |
723 | } | |
724 | htemp = hisSnp->Projection(1); | |
725 | Double_t xTheta=htemp->GetMean(); | |
726 | delete htemp; | |
727 | hisDelta = hisSnp->Projection(0); | |
728 | // | |
729 | Double_t entries = hisDelta->GetEntries(); | |
730 | Double_t mean=0, rms=0; | |
731 | if (entries>kMinEntries){ | |
732 | mean = hisDelta->GetMean(); | |
733 | rms = hisDelta->GetRMS(); | |
734 | } | |
735 | Double_t sector = 9.*xPhi/TMath::Pi(); | |
736 | if (sector<0) sector+=18; | |
737 | Double_t dsec = sector-Int_t(sector)-0.5; | |
738 | Int_t dtype=1; // TPC alignment type | |
739 | (*pcstream)<<hname<< | |
740 | "run="<<run<< | |
741 | "bz="<<bz<< // magnetic field | |
742 | "ptype="<<type<< // parameter type | |
743 | "dtype="<<dtype<< // parameter type | |
744 | "isec0="<<isec0<< // sector 0 | |
745 | "isec1="<<isec1<< // sector 1 | |
746 | "sector="<<sector<< // sector as float | |
747 | "dsec="<<dsec<< // delta sector | |
748 | // | |
749 | "theta="<<xTheta<< // theta | |
750 | "phi="<<xPhi<< // phi (alpha) | |
751 | "z="<<xZ<< // z | |
752 | "snp="<<xSnp<< // snp | |
753 | ||
754 | // | |
755 | "entries="<<entries<< // entries in bin | |
756 | "mean="<<mean<< // mean | |
757 | "rms="<<rms<< // rms | |
758 | "refX="<<refX<< // track matching reference plane | |
759 | "\n"; | |
760 | delete hisDelta; | |
761 | //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean); | |
762 | // | |
763 | }//ibinTheta | |
764 | delete hisSnp; | |
765 | } //ibinSnp | |
766 | delete hisZ; | |
767 | }//ibinZ | |
768 | delete hisPhi; | |
769 | }//ibinPhi | |
770 | delete hisSector1; | |
771 | }//isec1 | |
772 | delete hisSector0; | |
773 | }//isec0 | |
774 | } | |
775 | ||
776 | ||
777 | ||
778 | ||
779 | ||
780 | // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){ | |
781 | // // | |
782 | // // Make a laser fit tree for global minimization | |
783 | // // | |
784 | // AliTPCcalibDB* calib=AliTPCcalibDB::Instance(); | |
785 | // AliTPCCorrection * correction = calib->GetTPCComposedCorrection(); | |
786 | // if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz()); | |
787 | // correction->AddVisualCorrection(correction,0); //register correction | |
788 | ||
789 | // // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; | |
790 | // //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); | |
791 | // // | |
792 | // const Double_t cutErrY=0.05; | |
793 | // const Double_t kSigmaCut=4; | |
794 | // // const Double_t cutErrZ=0.03; | |
795 | // const Double_t kEpsilon=0.00000001; | |
796 | // // const Double_t kMaxDist=1.; // max distance - space correction | |
797 | // TVectorD *vecdY=0; | |
798 | // TVectorD *vecdZ=0; | |
799 | // TVectorD *veceY=0; | |
800 | // TVectorD *veceZ=0; | |
801 | // AliTPCLaserTrack *ltr=0; | |
802 | // AliTPCLaserTrack::LoadTracks(); | |
803 | // tree->SetBranchAddress("dY.",&vecdY); | |
804 | // tree->SetBranchAddress("dZ.",&vecdZ); | |
805 | // tree->SetBranchAddress("eY.",&veceY); | |
806 | // tree->SetBranchAddress("eZ.",&veceZ); | |
807 | // tree->SetBranchAddress("LTr.",<r); | |
808 | // Int_t entries= tree->GetEntries(); | |
809 | // TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root"); | |
810 | // Double_t bz=AliTrackerBase::GetBz(); | |
811 | // // | |
812 | // // Double_t globalXYZ[3]; | |
813 | // //Double_t globalXYZCorr[3]; | |
814 | // for (Int_t ientry=0; ientry<entries; ientry++){ | |
815 | // tree->GetEntry(ientry); | |
816 | // if (!ltr->GetVecGX()){ | |
817 | // ltr->UpdatePoints(); | |
818 | // } | |
819 | // // | |
820 | // TVectorD fit10(5); | |
821 | // TVectorD fit5(5); | |
822 | // printf("Entry\t%d\n",ientry); | |
823 | // for (Int_t irow0=0; irow0<158; irow0+=1){ | |
824 | // // | |
825 | // TLinearFitter fitter10(4,"hyp3"); | |
826 | // TLinearFitter fitter5(2,"hyp1"); | |
827 | // Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0]; | |
828 | // if (sector<0) continue; | |
829 | // //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue; | |
830 | ||
831 | // Double_t refX= (*ltr->GetVecLX())[irow0]; | |
832 | // Int_t firstRow1 = TMath::Max(irow0-10,0); | |
833 | // Int_t lastRow1 = TMath::Min(irow0+10,158); | |
834 | // Double_t padWidth=(irow0<64)?0.4:0.6; | |
835 | // // make long range fit | |
836 | // for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){ | |
837 | // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue; | |
838 | // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue; | |
839 | // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue; | |
840 | // Double_t idealX= (*ltr->GetVecLX())[irow1]; | |
841 | // Double_t idealY= (*ltr->GetVecLY())[irow1]; | |
842 | // // Double_t idealZ= (*ltr->GetVecLZ())[irow1]; | |
843 | // Double_t gx= (*ltr->GetVecGX())[irow1]; | |
844 | // Double_t gy= (*ltr->GetVecGY())[irow1]; | |
845 | // Double_t gz= (*ltr->GetVecGZ())[irow1]; | |
846 | // Double_t measY=(*vecdY)[irow1]+idealY; | |
847 | // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0); | |
848 | // // deltaR = R distorted -R ideal | |
849 | // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)}; | |
850 | // fitter10.AddPoint(xxx,measY,1); | |
851 | // } | |
852 | // Bool_t isOK=kTRUE; | |
853 | // Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4)); | |
854 | // Double_t mean10 =0;// fitter10.GetParameter(0); | |
855 | // Double_t slope10 =0;// fitter10.GetParameter(0); | |
856 | // Double_t cosPart10 = 0;// fitter10.GetParameter(2); | |
857 | // Double_t sinPart10 =0;// fitter10.GetParameter(3); | |
858 | ||
859 | // if (fitter10.GetNpoints()>10){ | |
860 | // fitter10.Eval(); | |
861 | // rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4)); | |
862 | // mean10 = fitter10.GetParameter(0); | |
863 | // slope10 = fitter10.GetParameter(1); | |
864 | // cosPart10 = fitter10.GetParameter(2); | |
865 | // sinPart10 = fitter10.GetParameter(3); | |
866 | // // | |
867 | // // make short range fit | |
868 | // // | |
869 | // for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){ | |
870 | // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue; | |
871 | // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue; | |
872 | // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue; | |
873 | // Double_t idealX= (*ltr->GetVecLX())[irow1]; | |
874 | // Double_t idealY= (*ltr->GetVecLY())[irow1]; | |
875 | // // Double_t idealZ= (*ltr->GetVecLZ())[irow1]; | |
876 | // Double_t gx= (*ltr->GetVecGX())[irow1]; | |
877 | // Double_t gy= (*ltr->GetVecGY())[irow1]; | |
878 | // Double_t gz= (*ltr->GetVecGZ())[irow1]; | |
879 | // Double_t measY=(*vecdY)[irow1]+idealY; | |
880 | // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0); | |
881 | // // deltaR = R distorted -R ideal | |
882 | // Double_t expY= mean10+slope10*(idealX+deltaR-refX); | |
883 | // if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue; | |
884 | // // | |
885 | // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth); | |
886 | // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)}; | |
887 | // fitter5.AddPoint(xxx,measY-corr,1); | |
888 | // } | |
889 | // }else{ | |
890 | // isOK=kFALSE; | |
891 | // } | |
892 | // if (fitter5.GetNpoints()<8) isOK=kFALSE; | |
893 | ||
894 | // Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4)); | |
895 | // Double_t offset5 =0;// fitter5.GetParameter(0); | |
896 | // Double_t slope5 =0;// fitter5.GetParameter(0); | |
897 | // if (isOK){ | |
898 | // fitter5.Eval(); | |
899 | // rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4)); | |
900 | // offset5 = fitter5.GetParameter(0); | |
901 | // slope5 = fitter5.GetParameter(0); | |
902 | // } | |
903 | // // | |
904 | // Double_t dtype=5; | |
905 | // Double_t ptype=0; | |
906 | // Double_t phi =(*ltr->GetVecPhi())[irow0]; | |
907 | // Double_t theta =ltr->GetTgl(); | |
908 | // Double_t mean=(vecdY)->GetMatrixArray()[irow0]; | |
909 | // Double_t gx=0,gy=0,gz=0; | |
910 | // Double_t snp = (*ltr->GetVecP2())[irow0]; | |
911 | // Int_t bundle= ltr->GetBundle(); | |
912 | // Int_t id= ltr->GetId(); | |
913 | // // Double_t rms = err->GetMatrixArray()[irow]; | |
914 | // // | |
915 | // gx = (*ltr->GetVecGX())[irow0]; | |
916 | // gy = (*ltr->GetVecGY())[irow0]; | |
917 | // gz = (*ltr->GetVecGZ())[irow0]; | |
918 | // Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0); | |
919 | // fitter10.GetParameters(fit10); | |
920 | // fitter5.GetParameters(fit5); | |
921 | // Double_t idealY= (*ltr->GetVecLY())[irow0]; | |
922 | // Double_t measY=(*vecdY)[irow0]+idealY; | |
923 | // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth); | |
924 | // if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE; | |
925 | // // | |
926 | // (*pcstream)<<"fitFull"<< // dumpe also intermediate results | |
927 | // "bz="<<bz<< // magnetic filed used | |
928 | // "dtype="<<dtype<< // detector match type | |
929 | // "ptype="<<ptype<< // parameter type | |
930 | // "theta="<<theta<< // theta | |
931 | // "phi="<<phi<< // phi | |
932 | // "snp="<<snp<< // snp | |
933 | // "sector="<<sector<< | |
934 | // "bundle="<<bundle<< | |
935 | // // // "dsec="<<dsec<< | |
936 | // "refX="<<refX<< // reference radius | |
937 | // "gx="<<gx<< // global position | |
938 | // "gy="<<gy<< // global position | |
939 | // "gz="<<gz<< // global position | |
940 | // "dRrec="<<dRrec<< // delta Radius in reconstruction | |
941 | // "id="<<id<< //bundle | |
942 | // "rms10="<<rms10<< | |
943 | // "rms5="<<rms5<< | |
944 | // "fit10.="<<&fit10<< | |
945 | // "fit5.="<<&fit5<< | |
946 | // "measY="<<measY<< | |
947 | // "mean="<<mean<< | |
948 | // "idealY="<<idealY<< | |
949 | // "corr="<<corr<< | |
950 | // "isOK="<<isOK<< | |
951 | // "\n"; | |
952 | // } | |
953 | // } | |
954 | // delete pcstream; | |
955 | // } | |
956 | ||
957 | ||
958 | void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){ | |
959 | // | |
960 | // Make a fit tree: | |
961 | // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX) | |
962 | // calculates partial distortions | |
963 | // Partial distortion is stored in the resulting tree | |
964 | // Output is storred in the file distortion_<dettype>_<partype>.root | |
965 | // Partial distortion is stored with the name given by correction name | |
966 | // | |
967 | // | |
968 | // Parameters of function: | |
969 | // input - input tree | |
970 | // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing | |
971 | // ppype - parameter type | |
972 | // corrArray - array with partial corrections | |
973 | // step - skipe entries - if 1 all entries processed - it is slow | |
974 | // debug 0 if debug on also space points dumped - it is slow | |
975 | ||
976 | const Double_t kMaxSnp = 0.85; | |
977 | const Double_t kcutSnp=0.25; | |
978 | const Double_t kcutTheta=1.; | |
979 | const Double_t kRadiusTPC=85; | |
980 | // AliTPCROC *tpcRoc =AliTPCROC::Instance(); | |
981 | // | |
982 | const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); | |
983 | // const Double_t kB2C=-0.299792458e-3; | |
984 | const Int_t kMinEntries=20; | |
985 | Double_t phi,theta, snp, mean,rms, entries,sector,dsec; | |
986 | Float_t refX; | |
987 | Int_t run; | |
988 | tinput->SetBranchAddress("run",&run); | |
989 | tinput->SetBranchAddress("theta",&theta); | |
990 | tinput->SetBranchAddress("phi", &phi); | |
991 | tinput->SetBranchAddress("snp",&snp); | |
992 | tinput->SetBranchAddress("mean",&mean); | |
993 | tinput->SetBranchAddress("rms",&rms); | |
994 | tinput->SetBranchAddress("entries",&entries); | |
995 | tinput->SetBranchAddress("sector",§or); | |
996 | tinput->SetBranchAddress("dsec",&dsec); | |
997 | tinput->SetBranchAddress("refX",&refX); | |
998 | TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset)); | |
999 | // | |
1000 | Int_t nentries=tinput->GetEntries(); | |
1001 | Int_t ncorr=corrArray->GetEntries(); | |
1002 | Double_t corrections[100]={0}; // | |
1003 | Double_t tPar[5]; | |
1004 | Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; | |
1005 | Int_t dir=0; | |
1006 | if (dtype==5 || dtype==6) dtype=4; | |
1007 | if (dtype==0) { dir=-1;} | |
1008 | if (dtype==1) { dir=1;} | |
1009 | if (dtype==2) { dir=-1;} | |
1010 | if (dtype==3) { dir=1;} | |
1011 | if (dtype==4) { dir=-1;} | |
1012 | // | |
1013 | for (Int_t ientry=offset; ientry<nentries; ientry+=step){ | |
1014 | tinput->GetEntry(ientry); | |
1015 | if (TMath::Abs(snp)>kMaxSnp) continue; | |
1016 | tPar[0]=0; | |
1017 | tPar[1]=theta*refX; | |
1018 | if (dtype==2) tPar[1]=theta*kRadiusTPC; | |
1019 | tPar[2]=snp; | |
1020 | tPar[3]=theta; | |
1021 | tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0 | |
1022 | if (dtype==4){ | |
1023 | // tracks crossing CE | |
1024 | tPar[1]=0; // track at the CE | |
1025 | //if (TMath::Abs(theta) <0.05) continue; // deep cross | |
1026 | } | |
1027 | ||
1028 | if (TMath::Abs(snp) >kcutSnp) continue; | |
1029 | if (TMath::Abs(theta) >kcutTheta) continue; | |
1030 | printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms); | |
1031 | Double_t bz=AliTrackerBase::GetBz(); | |
1032 | if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks | |
1033 | if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2); | |
1034 | ||
1035 | if (dtype==2 && TMath::Abs(bz)>0.1 ) { | |
1036 | tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);// | |
1037 | // snp at the TPC inner radius in case the vertex match used | |
1038 | } | |
1039 | } | |
1040 | // | |
1041 | tPar[4]+=(gRandom->Rndm()-0.5)*0.02; | |
1042 | AliExternalTrackParam track(refX,phi,tPar,cov); | |
1043 | Double_t xyz[3]; | |
1044 | track.GetXYZ(xyz); | |
1045 | Int_t id=0; | |
1046 | Double_t pt=1./tPar[4]; | |
1047 | Double_t dRrec=0; // dummy value - needed for points - e.g for laser | |
1048 | //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used | |
1049 | Double_t refXD=refX; | |
1050 | (*pcstream)<<"fit"<< | |
1051 | "run="<<run<< // run number | |
1052 | "bz="<<bz<< // magnetic filed used | |
1053 | "dtype="<<dtype<< // detector match type | |
1054 | "ptype="<<ptype<< // parameter type | |
1055 | "theta="<<theta<< // theta | |
1056 | "phi="<<phi<< // phi | |
1057 | "snp="<<snp<< // snp | |
1058 | "mean="<<mean<< // mean dist value | |
1059 | "rms="<<rms<< // rms | |
1060 | "sector="<<sector<< | |
1061 | "dsec="<<dsec<< | |
1062 | "refX="<<refXD<< // referece X as double | |
1063 | "gx="<<xyz[0]<< // global position at reference | |
1064 | "gy="<<xyz[1]<< // global position at reference | |
1065 | "gz="<<xyz[2]<< // global position at reference | |
1066 | "dRrec="<<dRrec<< // delta Radius in reconstruction | |
1067 | "pt="<<pt<< // pt | |
1068 | "id="<<id<< // track id | |
1069 | "entries="<<entries;// number of entries in bin | |
1070 | // | |
1071 | Bool_t isOK=kTRUE; | |
1072 | if (entries<kMinEntries) isOK=kFALSE; | |
1073 | // | |
1074 | if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) { | |
1075 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); | |
1076 | corrections[icorr]=0; | |
1077 | if (entries>kMinEntries){ | |
1078 | AliExternalTrackParam trackIn(refX,phi,tPar,cov); | |
1079 | AliExternalTrackParam *trackOut = 0; | |
1080 | if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream); | |
1081 | if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0); | |
1082 | if (dtype==0) {dir= -1;} | |
1083 | if (dtype==1) {dir= 1;} | |
1084 | if (dtype==2) {dir= -1;} | |
1085 | if (dtype==3) {dir= 1;} | |
1086 | // | |
1087 | if (trackOut){ | |
1088 | if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; | |
1089 | if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE; | |
1090 | if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1091 | // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz()); | |
1092 | // | |
1093 | corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype]; | |
1094 | delete trackOut; | |
1095 | }else{ | |
1096 | corrections[icorr]=0; | |
1097 | isOK=kFALSE; | |
1098 | } | |
1099 | //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out | |
1100 | } | |
1101 | (*pcstream)<<"fit"<< | |
1102 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value | |
1103 | } | |
1104 | ||
1105 | if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) { | |
1106 | // | |
1107 | // special case of the TPC tracks crossing the CE | |
1108 | // | |
1109 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); | |
1110 | corrections[icorr]=0; | |
1111 | if (entries>kMinEntries){ | |
1112 | AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex | |
1113 | AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet | |
1114 | AliExternalTrackParam *trackOut0 = 0; | |
1115 | AliExternalTrackParam *trackOut1 = 0; | |
1116 | // | |
1117 | if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream); | |
1118 | if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0); | |
1119 | if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream); | |
1120 | if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0); | |
1121 | // | |
1122 | if (trackOut0 && trackOut1){ | |
1123 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; | |
1124 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1125 | if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; | |
1126 | if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1127 | // | |
1128 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; | |
1129 | if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE; | |
1130 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1131 | if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE; | |
1132 | if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1133 | // | |
1134 | corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]); | |
1135 | corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]); | |
1136 | if (isOK) | |
1137 | if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)|| | |
1138 | (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)|| | |
1139 | (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)|| | |
1140 | (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)|| | |
1141 | (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)|| | |
1142 | (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001) | |
1143 | ){ | |
1144 | isOK=kFALSE; | |
1145 | } | |
1146 | delete trackOut0; | |
1147 | delete trackOut1; | |
1148 | }else{ | |
1149 | corrections[icorr]=0; | |
1150 | isOK=kFALSE; | |
1151 | } | |
1152 | // | |
1153 | //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup | |
1154 | } | |
1155 | (*pcstream)<<"fit"<< | |
1156 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value | |
1157 | } | |
1158 | // | |
1159 | (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n"; | |
1160 | } | |
1161 | ||
1162 | ||
1163 | delete pcstream; | |
1164 | } | |
1165 | ||
1166 | ||
1167 | ||
1168 | void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){ | |
1169 | // | |
1170 | // Make a fit tree: | |
1171 | // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX) | |
1172 | // calculates partial distortions | |
1173 | // Partial distortion is stored in the resulting tree | |
1174 | // Output is storred in the file distortion_<dettype>_<partype>.root | |
1175 | // Partial distortion is stored with the name given by correction name | |
1176 | // | |
1177 | // | |
1178 | // Parameters of function: | |
1179 | // input - input tree | |
1180 | // dtype - distortion type 10 - IROC-OROC | |
1181 | // ppype - parameter type | |
1182 | // corrArray - array with partial corrections | |
1183 | // step - skipe entries - if 1 all entries processed - it is slow | |
1184 | // debug 0 if debug on also space points dumped - it is slow | |
1185 | ||
1186 | const Double_t kMaxSnp = 0.8; | |
1187 | const Int_t kMinEntries=200; | |
1188 | // AliTPCROC *tpcRoc =AliTPCROC::Instance(); | |
1189 | // | |
1190 | const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); | |
1191 | // const Double_t kB2C=-0.299792458e-3; | |
1192 | Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ; | |
1193 | Int_t isec1, isec0; | |
1194 | Double_t refXD; | |
1195 | Float_t refX; | |
1196 | Int_t run; | |
1197 | tinput->SetBranchAddress("run",&run); | |
1198 | tinput->SetBranchAddress("theta",&theta); | |
1199 | tinput->SetBranchAddress("phi", &phi); | |
1200 | tinput->SetBranchAddress("snp",&snp); | |
1201 | tinput->SetBranchAddress("mean",&mean); | |
1202 | tinput->SetBranchAddress("rms",&rms); | |
1203 | tinput->SetBranchAddress("entries",&entries); | |
1204 | tinput->SetBranchAddress("sector",§or); | |
1205 | tinput->SetBranchAddress("dsec",&dsec); | |
1206 | tinput->SetBranchAddress("refX",&refXD); | |
1207 | tinput->SetBranchAddress("z",&globalZ); | |
1208 | tinput->SetBranchAddress("isec0",&isec0); | |
1209 | tinput->SetBranchAddress("isec1",&isec1); | |
1210 | TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset)); | |
1211 | // | |
1212 | Int_t nentries=tinput->GetEntries(); | |
1213 | Int_t ncorr=corrArray->GetEntries(); | |
1214 | Double_t corrections[100]={0}; // | |
1215 | Double_t tPar[5]; | |
1216 | Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; | |
1217 | Int_t dir=0; | |
1218 | // | |
1219 | for (Int_t ientry=offset; ientry<nentries; ientry+=step){ | |
1220 | tinput->GetEntry(ientry); | |
1221 | refX=refXD; | |
1222 | Int_t id=-1; | |
1223 | if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side | |
1224 | if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side | |
1225 | if (dtype==10 && id==-1) continue; | |
1226 | // | |
1227 | dir=-1; | |
1228 | tPar[0]=0; | |
1229 | tPar[1]=globalZ; | |
1230 | tPar[2]=snp; | |
1231 | tPar[3]=theta; | |
1232 | tPar[4]=(gRandom->Rndm()-0.1)*0.2; // | |
1233 | Double_t pt=1./tPar[4]; | |
1234 | // | |
1235 | printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms); | |
1236 | Double_t bz=AliTrackerBase::GetBz(); | |
1237 | AliExternalTrackParam track(refX,phi,tPar,cov); | |
1238 | Double_t xyz[3],xyzIn[3],xyzOut[3]; | |
1239 | track.GetXYZ(xyz); | |
1240 | track.GetXYZAt(85,bz,xyzIn); | |
1241 | track.GetXYZAt(245,bz,xyzOut); | |
1242 | Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]); | |
1243 | Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]); | |
1244 | Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]); | |
1245 | Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5); | |
1246 | Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5); | |
1247 | Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5); | |
1248 | // | |
1249 | Bool_t isOK=kTRUE; | |
1250 | if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector | |
1251 | if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector | |
1252 | if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin | |
1253 | // Do downscale | |
1254 | if (TMath::Abs(theta)>1) isOK=kFALSE; | |
1255 | // | |
1256 | Double_t dRrec=0; // dummy value - needed for points - e.g for laser | |
1257 | // | |
1258 | (*pcstream)<<"fit"<< | |
1259 | "run="<<run<< //run | |
1260 | "bz="<<bz<< // magnetic filed used | |
1261 | "dtype="<<dtype<< // detector match type | |
1262 | "ptype="<<ptype<< // parameter type | |
1263 | "theta="<<theta<< // theta | |
1264 | "phi="<<phi<< // phi | |
1265 | "snp="<<snp<< // snp | |
1266 | "mean="<<mean<< // mean dist value | |
1267 | "rms="<<rms<< // rms | |
1268 | "sector="<<sector<< | |
1269 | "dsec="<<dsec<< | |
1270 | "refX="<<refXD<< // referece X | |
1271 | "gx="<<xyz[0]<< // global position at reference | |
1272 | "gy="<<xyz[1]<< // global position at reference | |
1273 | "gz="<<xyz[2]<< // global position at reference | |
1274 | "dRrec="<<dRrec<< // delta Radius in reconstruction | |
1275 | "pt="<<pt<< //pt | |
1276 | "id="<<id<< // track id | |
1277 | "entries="<<entries;// number of entries in bin | |
1278 | // | |
1279 | AliExternalTrackParam *trackOut0 = 0; | |
1280 | AliExternalTrackParam *trackOut1 = 0; | |
1281 | AliExternalTrackParam *ptrackIn0 = 0; | |
1282 | AliExternalTrackParam *ptrackIn1 = 0; | |
1283 | ||
1284 | for (Int_t icorr=0; icorr<ncorr; icorr++) { | |
1285 | // | |
1286 | // special case of the TPC tracks crossing the CE | |
1287 | // | |
1288 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); | |
1289 | corrections[icorr]=0; | |
1290 | if (entries>kMinEntries &&isOK){ | |
1291 | AliExternalTrackParam trackIn0(refX,phi,tPar,cov); | |
1292 | AliExternalTrackParam trackIn1(refX,phi,tPar,cov); | |
1293 | ptrackIn1=&trackIn0; | |
1294 | ptrackIn0=&trackIn1; | |
1295 | // | |
1296 | if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream); | |
1297 | if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0); | |
1298 | if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream); | |
1299 | if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0); | |
1300 | // | |
1301 | if (trackOut0 && trackOut1){ | |
1302 | // | |
1303 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE; | |
1304 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1305 | // rotate all tracks to the same frame | |
1306 | if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; | |
1307 | if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE; | |
1308 | if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; | |
1309 | // | |
1310 | if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1311 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1312 | if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; | |
1313 | // | |
1314 | corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]); | |
1315 | corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]); | |
1316 | (*pcstream)<<"fitDebug"<< // just to debug the correction | |
1317 | "mean="<<mean<< | |
1318 | "pIn0.="<<ptrackIn0<< | |
1319 | "pIn1.="<<ptrackIn1<< | |
1320 | "pOut0.="<<trackOut0<< | |
1321 | "pOut1.="<<trackOut1<< | |
1322 | "refX="<<refXD<< | |
1323 | "\n"; | |
1324 | delete trackOut0; | |
1325 | delete trackOut1; | |
1326 | }else{ | |
1327 | corrections[icorr]=0; | |
1328 | isOK=kFALSE; | |
1329 | } | |
1330 | } | |
1331 | (*pcstream)<<"fit"<< | |
1332 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value | |
1333 | } | |
1334 | // | |
1335 | (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n"; | |
1336 | } | |
1337 | delete pcstream; | |
1338 | } | |
1339 | ||
1340 |