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