]> 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
Adding example hot to draw an example primitive transformation
[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
1cadb1b2 25
93096a07 262. Making fit trees
1cadb1b2 27
283. Space point distortion not directly observable. Instead a derived variables
29 like DCA at vertex, local y distortion in the TPC-*TOF,TRD,ITS) matching
30 in all 5 tracking parameters are obsereved.
31 In the AliTPCcorrection fir code we calculate the derivative of given variables
32 dO_{i}/dp_{i}
33
344. Global fit - later
35 d0 = sum{ki*dO_{i}/dp_{i}} - linear fitting of the amplitudes ki
93096a07 36
37 Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be
38 extracted from the old macros
39
40
41*/
42
43#include "Riostream.h"
44#include <fstream>
45#include "TMap.h"
46#include "TGraphErrors.h"
47#include "AliExternalTrackParam.h"
48#include "TFile.h"
49#include "TGraph.h"
50#include "TMultiGraph.h"
51#include "TCanvas.h"
52#include "THnSparse.h"
53#include "THnBase.h"
54#include "TProfile.h"
55#include "TROOT.h"
56#include "TLegend.h"
57#include "TPad.h"
58#include "TH2D.h"
59#include "TH3D.h"
60#include "AliTPCROC.h"
61#include "AliTPCCalROC.h"
62#include "AliESDfriend.h"
63#include "AliTPCcalibTime.h"
64#include "AliSplineFit.h"
65#include "AliCDBMetaData.h"
66#include "AliCDBId.h"
67#include "AliCDBManager.h"
68#include "AliCDBStorage.h"
69#include "AliTPCcalibBase.h"
70#include "AliTPCcalibDB.h"
71#include "AliTPCcalibDButil.h"
72#include "AliRelAlignerKalman.h"
73#include "AliTPCParamSR.h"
74#include "AliTPCcalibTimeGain.h"
75#include "AliTPCcalibGainMult.h"
76#include "AliSplineFit.h"
77#include "AliTPCComposedCorrection.h"
78#include "AliTPCExBTwist.h"
79#include "AliTPCCalibGlobalMisalignment.h"
80#include "TStatToolkit.h"
81#include "TChain.h"
82#include "TCut.h"
83#include "AliTrackerBase.h"
84#include "AliTPCCorrectionFit.h"
85#include "AliTPCLaserTrack.h"
86#include "TDatabasePDG.h"
87#include "AliTPCcalibAlign.h"
4425a337 88#include "AliLog.h"
1cadb1b2 89#include "AliRieman.h"
93096a07 90
91ClassImp(AliTPCCorrectionFit)
92
93AliTPCCorrectionFit::AliTPCCorrectionFit():
94 TNamed("TPCCorrectionFit","TPCCorrectionFit")
95{
96 //
97 // default constructor
98 //
99}
100
101AliTPCCorrectionFit::~AliTPCCorrectionFit() {
102 //
103 // Destructor
104 //
105}
106
107
108Double_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
124Double_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 172Double_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
221void 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
256void 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
445void 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
536void 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
647void 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.",&ltr);
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
1026void 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",&sector);
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
1236void 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",&sector);
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