]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCorrectionFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updates from Saehanseul.
[u/mrichter/AliRoot.git] / TPC / AliTPCCorrectionFit.cxx
CommitLineData
93096a07 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17
18/*
19 Responsible: marian.ivanov@cern.ch
20 Tools for fitting of the space point distortion parameters.
21 Functionality
22
23
241. Creation of the distortion maps from the residual histograms
252. Making fit trees
263. 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
81ClassImp(AliTPCCorrectionFit)
82
83AliTPCCorrectionFit::AliTPCCorrectionFit():
84 TNamed("TPCCorrectionFit","TPCCorrectionFit")
85{
86 //
87 // default constructor
88 //
89}
90
91AliTPCCorrectionFit::~AliTPCCorrectionFit() {
92 //
93 // Destructor
94 //
95}
96
97
98Double_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
114Double_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
153void 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
188void 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
377void 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
468void 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
579void 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.",&ltr);
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
958void 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",&sector);
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
1168void 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",&sector);
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