]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCcalibUnlinearity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Mac users should be happy now. Thanks to Laurent for providing the fix
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibUnlinearity.cxx
... / ...
CommitLineData
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 Class for histogramming of cluster residuals
18 and fitting of the unlinearities. To be used only for data without
19 magnetic field. The laser tracks should be also rejected.
20 //
21
22 Track fitting:
23 The track is fitted using linear and parabolic model
24 The edge clusters are removed from the fit.
25 Edge pad-row - first and last 15 padrows removed from the fit
26
27 Unlinearities fitting:
28 Unlinearities at the edge aproximated using two exponential decays.
29
30 Model:
31 dz = dz0(r,z) +dr(r,z)*tan(theta)
32 dy = +dr(r,z)*tan(phi)
33
34
35 .x ~/NimStyle.C
36 gSystem->Load("libANALYSIS");
37 gSystem->Load("libTPCcalib");
38 TFile fcalib("CalibObjects.root");
39 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
40 AliTPCcalibUnlinearity * calibUnlin = ( AliTPCcalibUnlinearity *)array->FindObject("calibUnlinearity");
41 //
42
43*/
44
45
46#include "TLinearFitter.h"
47
48#include "Riostream.h"
49#include "TChain.h"
50#include "TTree.h"
51#include "TH1F.h"
52#include "TH2F.h"
53#include "TH3F.h"
54#include "THnSparse.h"
55#include "TList.h"
56#include "TMath.h"
57#include "TCanvas.h"
58#include "TFile.h"
59#include "TF1.h"
60#include "TVectorD.h"
61#include "TProfile.h"
62#include "TGraphErrors.h"
63#include "TCanvas.h"
64#include "TCut.h"
65
66
67#include "AliTPCclusterMI.h"
68#include "AliTPCseed.h"
69#include "AliESDVertex.h"
70#include "AliESDEvent.h"
71#include "AliESDfriend.h"
72#include "AliESDInputHandler.h"
73#include "AliAnalysisManager.h"
74#include "AliTracker.h"
75#include "AliMagF.h"
76#include "AliTPCCalROC.h"
77
78#include "AliLog.h"
79
80
81#include "TTreeStream.h"
82#include "AliTPCTracklet.h"
83#include "TTimeStamp.h"
84#include "AliTPCcalibDB.h"
85#include "AliDCSSensorArray.h"
86#include "AliDCSSensor.h"
87#include "TLinearFitter.h"
88#include "AliTPCClusterParam.h"
89#include "TStatToolkit.h"
90
91
92#include "AliTPCcalibUnlinearity.h"
93#include "AliTPCPointCorrection.h"
94
95
96ClassImp(AliTPCcalibUnlinearity)
97
98AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
99 AliTPCcalibBase(),
100 fInit(kFALSE),
101 fDiffHistoLineY(0),
102 fDiffHistoParY(0),
103 fDiffHistoLineZ(0),
104 fDiffHistoParZ(0),
105 fFittersOutR(38),
106 fFittersOutZ(38),
107 fParamsOutR(38),
108 fParamsOutZ(38),
109 fErrorsOutR(38),
110 fErrorsOutZ(38),
111 fDistRPHIPlus(74),
112 fDistRPHIMinus(74),
113 fFitterQuadrantY(16*38),
114 fFitterQuadrantPhi(16*38)
115{
116 //
117 // Defualt constructor
118 //
119}
120
121AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
122 AliTPCcalibBase(name,title),
123 fInit(kFALSE),
124 fDiffHistoLineY(0),
125 fDiffHistoParY(0),
126 fDiffHistoLineZ(0),
127 fDiffHistoParZ(0),
128 fFittersOutR(38),
129 fFittersOutZ(38),
130 fParamsOutR(38),
131 fParamsOutZ(38),
132 fErrorsOutR(38),
133 fErrorsOutZ(38),
134 fDistRPHIPlus(74),
135 fDistRPHIMinus(74),
136 fFitterQuadrantY(16*38),
137 fFitterQuadrantPhi(16*38)
138{
139 //
140 // Non default constructor
141 //
142 MakeHisto();
143}
144
145AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
146 //
147 //
148 //
149 if (fDiffHistoLineZ) delete fDiffHistoLineY;
150 if (fDiffHistoParZ) delete fDiffHistoParY;
151 if (fDiffHistoLineZ) delete fDiffHistoLineZ;
152 if (fDiffHistoParZ) delete fDiffHistoParZ;
153}
154
155
156
157
158
159void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
160 //
161 //
162 //
163 if (!fInit) {
164 Init(); // work around for PROOF - initialize firs time used
165 }
166 const Int_t kdrow=10;
167 const Int_t kMinCluster=35;
168 if (!seed) return;
169 if (TMath::Abs(fMagF)>0.01) return; // use only data without field
170 Int_t counter[72];
171 for (Int_t i=0;i<72;i++) counter[i]=0;
172 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
173 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
174 if (!cluster) continue;
175 if (cluster->GetType()<0) continue;
176 counter[cluster->GetDetector()]++;
177 }
178 //
179 for (Int_t is=0; is<72;is++){
180 if (counter[is]>kMinCluster) ProcessDiff(seed, is);
181 if (counter[is]>kMinCluster) {
182 AlignOROC(seed, is);
183 }
184 }
185}
186
187
188void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
189 //
190 // process differnce of the cluster in respect with linear and parabolic fit
191 //
192 const Double_t kXmean=134;
193 const Int_t kdrow=10;
194 const Float_t kSagitaCut = 1;
195 static TLinearFitter fy1(2,"pol1");
196 static TLinearFitter fy2(3,"pol2");
197 static TLinearFitter fz1(2,"pol1");
198 static TLinearFitter fz2(3,"pol2");
199 //
200 static TVectorD py1(2);
201 static TVectorD py2(3);
202 static TVectorD pz1(2);
203 static TVectorD pz2(3);
204 //
205 fy1.ClearPoints();
206 fy2.ClearPoints();
207 fz1.ClearPoints();
208 fz2.ClearPoints();
209 //
210 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
211 AliTPCclusterMI *c=track->GetClusterPointer(irow);
212 if (!c) continue;
213 if (c->GetDetector()!=isec) continue;
214 if (c->GetType()<0) continue;
215 Double_t dx = c->GetX()-kXmean;
216 Double_t x[2]={dx, dx*dx};
217 fy1.AddPoint(x,c->GetY());
218 fy2.AddPoint(x,c->GetY());
219 fz1.AddPoint(x,c->GetZ());
220 fz2.AddPoint(x,c->GetZ());
221 }
222 fy1.Eval(); fy1.GetParameters(py1);
223 fy2.Eval(); fy2.GetParameters(py2);
224 fz1.Eval(); fz1.GetParameters(pz1);
225 fz2.Eval(); fz2.GetParameters(pz2);
226 //
227
228 for (Int_t irow=0;irow<159;irow++) {
229 AliTPCclusterMI *c=track->GetClusterPointer(irow);
230 if (!c) continue;
231 if (c->GetDetector()!=isec) continue;
232 Double_t dx = c->GetX()-kXmean;
233 Double_t y1 = py1[0]+py1[1]*dx;
234 Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
235 //
236 Double_t z1 = pz1[0]+pz1[1]*dx;
237 Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
238 //
239 Double_t edgeMinus= -y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
240 Double_t edgePlus = y1+c->GetX()*TMath::Tan(TMath::Pi()/18.);
241 Int_t npads = AliTPCROC::Instance()->GetNPads(isec,c->GetRow());
242 Float_t dpad = TMath::Min(c->GetPad(), npads-1- c->GetPad());
243 Float_t dy =c->GetY()-y1;
244 //
245 // Corrections
246 //
247 AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
248 Double_t corrclY=0, corrtrY=0, corrR=0;
249 corrclY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
250 c->GetY(),c->GetY(), c->GetZ(), py1[1],c->GetMax(),2.5);
251 corrtrY = corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
252 c->GetY(),y1, c->GetZ(), py1[1], c->GetMax(),2.5);
253 corrR = corr->CorrectionOutR0(kFALSE,kFALSE,c->GetX(),c->GetY(),c->GetZ(),isec);
254 //
255 //
256 //
257 if (fStreamLevel>1){
258 TTreeSRedirector *cstream = GetDebugStreamer();
259 if (cstream){
260 (*cstream)<<"Diff"<<
261 "run="<<fRun<< // run number
262 "event="<<fEvent<< // event number
263 "time="<<fTime<< // time stamp of event
264 "trigger="<<fTrigger<< // trigger
265 "mag="<<fMagF<< // magnetic field
266 "isec="<<isec<< // current sector
267 "npads="<<npads<< // number of pads at given sector
268 "dpad="<<dpad<< // distance to edge pad
269 //
270 "crR="<<corrR<< // radial correction
271 "cclY="<<corrclY<< // edge effect correction cl
272 "ctrY="<<corrtrY<< // edge effect correction using track
273 //
274 "Cl.="<<c<<
275 "y1="<<y1<<
276 "y2="<<y2<<
277 "z1="<<z1<<
278 "z2="<<z2<<
279 "py1.="<<&py1<<
280 "py2.="<<&py2<<
281 "pz1.="<<&pz1<<
282 "pz2.="<<&pz2<<
283 "eP="<<edgePlus<<
284 "eM="<<edgeMinus<<
285 "\n";
286 }
287 }
288 if (TMath::Min(edgeMinus,edgePlus)<6){
289 AddPointRPHI(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
290 TTreeSRedirector *cstream = GetDebugStreamer();
291 if (TMath::Abs(dy)<2 && cstream){
292 (*cstream)<<"rphi"<<
293 "isec="<<isec<< // current sector
294 "npads="<<npads<< // number of pads at given sector
295 "dpad="<<dpad<< // distance to edge pad
296 "eP="<<edgePlus<< // distance to edge
297 "eM="<<edgeMinus<< // distance to edge minus
298 //
299 "crR="<<corrR<< // radial correction
300 "cclY="<<corrclY<< // edge effect correction cl
301 "ctrY="<<corrtrY<< // edge effect correction using track
302 //
303 "dy="<<dy<< // dy
304 "Cl.="<<c<<
305 "y1="<<y1<<
306 "y2="<<y2<<
307 "z1="<<z1<<
308 "z2="<<z2<<
309 "py1.="<<&py1<<
310 "pz1.="<<&pz1<<
311 "\n";
312 }
313 }
314 if (c->GetType()<0) continue; // don't use edge rphi cluster
315 //
316 //
317 Double_t x[10];
318 x[0]=c->GetDetector();
319 x[1]=c->GetRow();
320 x[2]=c->GetY()/c->GetX();
321 x[3]=c->GetZ();
322 //
323 // apply sagita cut
324 //
325 if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut &&
326 TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
327 //
328 x[4]=py1[1];
329 x[5]=c->GetY()-y1;
330 fDiffHistoLineY->Fill(x);
331 x[5]=c->GetY()-y2;
332 //fDiffHistoParY->Fill(x);
333 //
334 x[4]=pz1[1];
335 x[5]=c->GetY()-z1;
336 fDiffHistoLineZ->Fill(x);
337 x[5]=c->GetY()-z2;
338 //fDiffHistoParZ->Fill(x);
339 // Apply sagita cut
340 AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
341 }
342
343 }
344}
345
346
347void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
348 //
349 // fit the tracklet in OROC sepatately in Quadrant
350 //
351 //
352 if (isec<36) return;
353 const Int_t kMinClusterF=35;
354 const Int_t kMinClusterQ=10;
355 //
356 const Int_t kdrow1 =8; // rows to skip at the end
357 const Int_t kdrow0 =2; // rows to skip at beginning
358 const Float_t kedgey=3;
359 const Float_t kMaxDist=0.5;
360 const Float_t kMaxCorrY=0.1;
361 const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
362 //
363 //
364 AliTPCROC * roc = AliTPCROC::Instance();
365 AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
366 const Double_t kXmean=roc->GetPadRowRadii(isec,53);
367 //
368 // full track fit parameters
369 //
370 TLinearFitter fyf(2,"pol1");
371 TLinearFitter fzf(2,"pol1");
372 TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
373 Int_t nf=0;
374 //
375 // make full fit as reference
376 //
377 for (Int_t iter=0; iter<2; iter++){
378 fyf.ClearPoints();
379 for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
380 AliTPCclusterMI *c=track->GetClusterPointer(irow);
381 if (!c) continue;
382 if (c->GetDetector()!=isec) continue;
383 if (c->GetRow()<kdrow0) continue;
384 Double_t dx = c->GetX()-kXmean;
385 Double_t x[2]={dx, dx*dx};
386 if (iter==0 &&c->GetType()<0) continue;
387 if (iter==1){
388 Double_t yfit = pyf[0]+pyf[1]*dx;
389 Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
390 if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
391 if (dedge<kedgey) continue;
392 Double_t corrtrY =
393 corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
394 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
395 if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
396 }
397 fyf.AddPoint(x,c->GetY(),0.1);
398 fzf.AddPoint(x,c->GetZ(),0.1);
399 }
400 nf = fyf.GetNpoints();
401 if (nf<kMinClusterF) return; // not enough points - skip
402 fyf.Eval();
403 fyf.GetParameters(pyf);
404 fyf.GetErrors(peyf);
405 fzf.Eval();
406 fzf.GetParameters(pzf);
407 fzf.GetErrors(pezf);
408 }
409 //
410 // Make Fitters and params for 3 fitters
411 //
412 TLinearFitter *fitters[4];
413 Int_t npoints[4];
414 TVectorD params[4];
415 TVectorD errors[4];
416 Double_t chi2C[4];
417 for (Int_t i=0;i<4;i++) {
418 fitters[i] = new TLinearFitter(2,"pol1");
419 npoints[i]=0;
420 params[i].ResizeTo(2);
421 errors[i].ResizeTo(2);
422 }
423 //
424 // Update fitters
425 //
426 for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
427 AliTPCclusterMI *c=track->GetClusterPointer(irow);
428 if (!c) continue;
429 if (c->GetDetector()!=isec) continue;
430 if (c->GetRow()<kdrow0) continue;
431 Double_t dx = c->GetX()-kXmean;
432 Double_t x[2]={dx, dx*dx};
433 Double_t yfit = pyf[0]+pyf[1]*dx;
434 Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
435 if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
436 if (dedge<kedgey) continue;
437 Double_t corrtrY =
438 corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
439 c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
440 if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
441 if (dx<0){
442 if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
443 if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
444 }
445 if (dx>0){
446 if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
447 if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
448 }
449 }
450
451 for (Int_t i=0;i<4;i++) {
452 npoints[i] = fitters[i]->GetNpoints();
453 if (npoints[i]>=kMinClusterQ){
454 fitters[i]->Eval();
455 Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
456 chi2C[i]=chi2Fac;
457 fitters[i]->GetParameters(params[i]);
458 fitters[i]->GetErrors(errors[i]);
459 // renormalize errors
460 errors[i][0]*=chi2Fac;
461 errors[i][1]*=chi2Fac;
462 }
463 }
464 //
465 // Fill fitters
466 //
467 for (Int_t i0=0;i0<4;i0++){
468 for (Int_t i1=0;i1<4;i1++){
469 if (i0==i1) continue;
470 if(npoints[i0]<kMinClusterQ) continue;
471 if(npoints[i1]<kMinClusterQ) continue;
472 Int_t index0=i0*4+i1;
473 Double_t xy[1];
474 Double_t xfi[1];
475 xy[0] = pyf[1];
476 xfi[0] = pyf[1];
477 //
478 Double_t dy = params[i1][0]-params[i0][0];
479 Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
480 Double_t dphi = params[i1][1]-params[i0][1];
481 Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
482 //
483 Int_t sector = isec-36;
484 //
485 if (TMath::Abs(dy/(sy+0.1))>5.) continue; // 5 sigma cut
486 if (TMath::Abs(dphi/(sphi+0.004))>5.) continue; // 5 sigma cut
487 TLinearFitter * fitterY,*fitterPhi;
488 fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
489 if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
490 fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
491 if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
492 //
493 fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
494 if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
495 fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
496 if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
497 }
498 }
499 //
500 // Dump debug information
501 //
502 if (fStreamLevel>0){
503 TTreeSRedirector *cstream = GetDebugStreamer();
504 Int_t sector = isec-36;
505 if (cstream){
506 for (Int_t i0=0;i0<4;i0++){
507 for (Int_t i1=i0;i1<4;i1++){
508 if (i0==i1) continue;
509 if(npoints[i0]<kMinClusterQ) continue;
510 if(npoints[i1]<kMinClusterQ) continue;
511 (*cstream)<<"fitQ"<<
512 "run="<<fRun<< // run number
513 "event="<<fEvent<< // event number
514 "time="<<fTime<< // time stamp of event
515 "trigger="<<fTrigger<< // trigger
516 "mag="<<fMagF<< // magnetic field
517 "sec="<<sector<< // current sector from 0 to 36
518 "isec="<<isec<< // current sector
519 // full fit
520 "nf="<<nf<< // total number of points
521 "pyf.="<<&pyf<< // full OROC fit y
522 "pzf.="<<&pzf<< // full OROC fit z
523 // quadrant fit
524 "i0="<<i0<< // quadrant number
525 "i1="<<i1<<
526 "n0="<<npoints[i0]<< // number of points
527 "n1="<<npoints[i1]<<
528 "p0.="<<&params[i0]<< // parameters
529 "p1.="<<&params[i1]<<
530 "e0.="<<&errors[i0]<< // errors
531 "e1.="<<&errors[i1]<<
532 "chi0="<<chi2C[i0]<< // chi2s
533 "chi1="<<chi2C[i1]<<
534
535 "\n";
536 }
537 }
538 }
539 }
540}
541
542
543
544
545
546void AliTPCcalibUnlinearity::MakeHisto(){
547 //
548 //
549 //
550 TString axisName[10];
551 Double_t xmin[10], xmax[10];
552 Int_t nbins[10];
553 //
554 //
555 axisName[0]="sector";
556 xmin[0]=0; xmax[0]=72; nbins[0]=72;
557 //
558 axisName[1]="row";
559 xmin[1]=0; xmax[1]=96; nbins[1]=96;
560 //
561 axisName[2]="tphi"; // tan phi - ly/lx
562 xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
563 //
564 axisName[3]="lz"; // global z
565 xmin[3]=-250; xmax[3]=250; nbins[3]=10;
566 //
567 axisName[4]="k"; // tangent - in angle
568 xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
569 //
570 //
571 axisName[5]="delta"; // delta
572 xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
573 //
574 //
575 fDiffHistoLineY = new THnSparseF("hnDistHistoLineY","hnDistHistoLineY",6,nbins,xmin,xmax);
576 fDiffHistoParY = new THnSparseF("hnDistHistoParY","hnDistHistoParY",6,nbins,xmin,xmax);
577 fDiffHistoLineZ = new THnSparseF("hnDistHistoLineZ","hnDistHistoLineZ",6,nbins,xmin,xmax);
578 fDiffHistoParZ = new THnSparseF("hnDistHistoParZ","hnDistHistoParZ",6,nbins,xmin,xmax);
579 // for (Int_t i=0; i<8;i++){
580// fDiffHistoLineY->GetAxis(i)->SetName(axisName[i].Data());
581// fDiffHistoParY->GetAxis(i)->SetName(axisName[i].Data());
582// fDiffHistoLineY->GetAxis(i)->SetTitle(axisName[i].Data());
583// fDiffHistoParY->GetAxis(i)->SetTitle(axisName[i].Data());
584// fDiffHistoLineZ->GetAxis(i)->SetName(axisName[i].Data());
585// fDiffHistoParZ->GetAxis(i)->SetName(axisName[i].Data());
586// fDiffHistoLineZ->GetAxis(i)->SetTitle(axisName[i].Data());
587// fDiffHistoParZ->GetAxis(i)->SetTitle(axisName[i].Data());
588// }
589
590 //
591 //
592 //
593 char hname[1000];
594 TH2F * his=0;
595 for (Int_t isec=0;isec<74;isec++){
596 sprintf(hname,"DeltaRPhiPlus_Sector%d",isec);
597 his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
598 his->SetDirectory(0);
599 fDistRPHIPlus.AddAt(his,isec);
600 sprintf(hname,"DeltaRPhiMinus_Sector%d",isec);
601 his = new TH2F(hname,hname,100,1,4,100,-0.5,0.5);
602 his->SetDirectory(0);
603 fDistRPHIMinus.AddAt(his,isec);
604 }
605}
606
607
608void AliTPCcalibUnlinearity::Terminate(){
609 //
610 // Terminate function
611 // call base terminate + Eval of fitters
612 //
613 Info("AliTPCcalibUnlinearities","Terminate");
614 EvalFitters();
615 DumpTree();
616 AliTPCcalibBase::Terminate();
617}
618
619Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
620 //
621 // merge results
622 //
623 TIterator* iter = list->MakeIterator();
624 AliTPCcalibUnlinearity* cal = 0;
625 //
626 while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
627 if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
628 return -1;
629 }
630 Add(cal);
631 }
632 return 0;
633}
634
635void AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
636 //
637 //
638 //
639 if (!fInit) Init();
640 if (fDiffHistoLineY && calib->fDiffHistoLineY){
641 fDiffHistoLineY->Add(calib->fDiffHistoLineY);
642 }
643 if (fDiffHistoParY && calib->fDiffHistoParY){
644 fDiffHistoParY->Add(calib->fDiffHistoParY);
645 }
646 if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
647 fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
648 }
649 if (fDiffHistoParZ && calib->fDiffHistoParZ){
650 fDiffHistoParZ->Add(calib->fDiffHistoParZ);
651 }
652 for (Int_t isec=0;isec<38;isec++){
653 TLinearFitter *f0r = (TLinearFitter*)fFittersOutR.At(isec);
654 TLinearFitter *f1r = (TLinearFitter*)(calib->fFittersOutR.At(isec));
655 if (f0r&&f1r) f0r->Add(f1r);
656 TLinearFitter *f0z = (TLinearFitter*)fFittersOutZ.At(isec);
657 TLinearFitter *f1z = (TLinearFitter*)(calib->fFittersOutZ.At(isec));
658 if (f0z&&f1z) f0z->Add(f1z);
659 }
660
661 for (Int_t isec=0;isec<16*38;isec++){
662 TLinearFitter *f0y = (TLinearFitter*)fFitterQuadrantY.At(isec);
663 TLinearFitter *f1y = (TLinearFitter*)(calib->fFitterQuadrantY.At(isec));
664 if (f0y&&f1y) f0y->Add(f1y);
665 TLinearFitter *f0phi = (TLinearFitter*)fFitterQuadrantPhi.At(isec);
666 TLinearFitter *f1phi = (TLinearFitter*)(calib->fFitterQuadrantPhi.At(isec));
667 if (f0phi&&f1phi) f0phi->Add(f1phi);
668 }
669
670
671 for (Int_t isec=0;isec<74;isec++){
672 TH2F * h0p= (TH2F*)(fDistRPHIPlus.At(isec));
673 TH2F * h1p= (TH2F*)(calib->fDistRPHIPlus.At(isec));
674 if (h0p&&h1p) h0p->Add(h1p);
675 TH2F * h0m= (TH2F*)(fDistRPHIMinus.At(isec));
676 TH2F * h1m= (TH2F*)(calib->fDistRPHIMinus.At(isec));
677 if (h0m&&h1m) h0m->Add(h1m);
678 }
679
680
681}
682
683
684
685void AliTPCcalibUnlinearity::DumpTree(const char *fname){
686 //
687 //
688 //
689 TTreeSRedirector *cstream =new TTreeSRedirector(fname);
690 if (!cstream) return;
691 //
692 THnSparse *his=0;
693 Double_t position[10];
694 Double_t value;
695 Int_t *bins = new Int_t[10];
696 //
697 for (Int_t ihis=0; ihis<2; ihis++){
698 his = (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
699 if (!his) continue;
700 //
701 Int_t ndim = his->GetNdimensions();
702 //
703 for (Long64_t i = 0; i < his->GetNbins(); ++i) {
704 value = his->GetBinContent(i, bins);
705 for (Int_t idim = 0; idim < ndim; idim++) {
706 position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
707 }
708 (*cstream)<<"Resol"<<
709 "ihis="<<ihis<<
710 "bincont="<<value;
711 for (Int_t idim=0;idim<ndim;idim++){
712 (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
713 }
714 (*cstream)<<"Resol"<<
715 "\n";
716 }
717 }
718 delete cstream;
719}
720
721
722void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t ky, Double_t kz, Int_t npoints){
723 //
724 //
725 // Simple distortion fit in outer sectors
726 //
727 // Model - 2 exponential decay toward the center of chamber
728 // - Decay length 10 and 5 cm
729 // - Decay amplitude - Z dependent
730 //
731
732 Double_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1
733 Double_t dzt = (cz-tz)*side;
734 Double_t radius = TMath::Sqrt(cx*cx+cy*cy);
735 AliTPCROC * roc = AliTPCROC::Instance();
736 Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
737 Double_t dout = xout-radius;
738 if (dout>30) return;
739 //drift
740 Double_t dr = 0.5 - TMath::Abs(cz/250.);
741 Double_t dy = cy-ty;
742 Double_t eout10 = TMath::Exp(-dout/10.);
743 Double_t eout5 = TMath::Exp(-dout/5.);
744 //
745 Double_t eoutp = (eout10+eout5)*0.5;
746 Double_t eoutm = (eout10-eout5)*0.5;
747
748 //
749 Double_t xxxR[6], xxxZ[6], xxxRZ[6];
750 //TString fstring="";
751 xxxZ[0]=eoutp; //fstring+="eoutp++";
752 xxxZ[1]=eoutm; //fstring+="eoutm++";
753 xxxZ[2]=dr*eoutp; //fstring+="dr*eoutp++";
754 xxxZ[3]=dr*eoutm; //fstring+="dr*eoutm++";
755 xxxZ[4]=dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
756 xxxZ[5]=dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
757 //
758 xxxR[0]=ky*eoutp; //fstring+="eoutp++";
759 xxxR[1]=ky*eoutm; //fstring+="eoutm++";
760 xxxR[2]=ky*dr*eoutp; //fstring+="dr*eoutp++";
761 xxxR[3]=ky*dr*eoutm; //fstring+="dr*eoutm++";
762 xxxR[4]=ky*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
763 xxxR[5]=ky*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
764 //
765 xxxRZ[0]=kz*eoutp; //fstring+="eoutp++";
766 xxxRZ[1]=kz*eoutm; //fstring+="eoutm++";
767 xxxRZ[2]=kz*dr*eoutp; //fstring+="dr*eoutp++";
768 xxxRZ[3]=kz*dr*eoutm; //fstring+="dr*eoutm++";
769 xxxRZ[4]=kz*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
770 xxxRZ[5]=kz*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
771 //
772 TLinearFitter * fitter=0;
773 Double_t err=0.1;
774 for (Int_t ip=0; ip<npoints; ip++){
775 //
776 fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
777 fitter->AddPoint(xxxR,dy,err);
778 //fitter->AddPoint(xxxRZ,dz);
779 fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
780 fitter->AddPoint(xxxZ,dzt,err);
781 //
782 if (side>0){
783 fitter = (TLinearFitter*)fFittersOutR.At(36);
784 fitter->AddPoint(xxxR,dy,err);
785 //fitter->AddPoint(xxxRZ,dz);
786 fitter = (TLinearFitter*)fFittersOutZ.At(36);
787 fitter->AddPoint(xxxZ,dzt,err);
788 }
789 if (side<0){
790 fitter = (TLinearFitter*)fFittersOutR.At(37);
791 fitter->AddPoint(xxxR,dy,err);
792 //fitter->AddPoint(xxxRZ,dz);
793 fitter = (TLinearFitter*)fFittersOutZ.At(37);
794 fitter->AddPoint(xxxZ,dzt,err);
795 }
796 }
797}
798void AliTPCcalibUnlinearity::AddPointRPHI(Int_t sector, Double_t cx, Double_t cy, Double_t cz, Double_t ty, Double_t tz, Double_t /*ky*/, Double_t /*kz*/, Int_t npoints){
799 //
800 //
801 //
802 Float_t kMaxDz = 0.5; // cut on z in cm
803 Double_t dy = cy-ty;
804 Double_t dz = cz-tz;
805 if (TMath::Abs(dz)>kMaxDz) return;
806 Double_t edgeMinus= -ty+cx*TMath::Tan(TMath::Pi()/18.);
807 Double_t edgePlus = ty+cx*TMath::Tan(TMath::Pi()/18.);
808 //
809 TH2F * his =0;
810 his = (TH2F*)fDistRPHIPlus.At(sector);
811 his->Fill(edgePlus,dy,npoints);
812 his = (TH2F*)((sector<36)? fDistRPHIPlus.At(72):fDistRPHIPlus.At(73));
813 his->Fill(edgePlus,dy,npoints);
814 his = (TH2F*)fDistRPHIMinus.At(sector);
815 his->Fill(edgeMinus,dy,npoints);
816 his = (TH2F*)((sector<36)? fDistRPHIMinus.At(72):fDistRPHIMinus.At(73));
817 his->Fill(edgeMinus,dy,npoints);
818}
819
820
821
822void AliTPCcalibUnlinearity::Init(){
823 //
824 //
825 //
826 //
827 MakeHisto();
828 // Make outer fitters
829 TLinearFitter * fitter = 0;
830 for (Int_t ifit=0; ifit<38; ifit++){
831 fitter = new TLinearFitter(7,"hyp6");
832 fitter->StoreData(kFALSE);
833 fFittersOutR.AddAt(fitter,ifit);
834 fitter = new TLinearFitter(7,"hyp6");
835 fitter->StoreData(kFALSE);
836 fFittersOutZ.AddAt(fitter,ifit);
837 }
838 for (Int_t ifit=0; ifit<16*38;ifit++){
839 fitter = new TLinearFitter(2,"hyp1");
840 fitter->StoreData(kFALSE);
841 fFitterQuadrantY.AddAt(fitter,ifit);
842 fitter = new TLinearFitter(2,"hyp1");
843 fitter->StoreData(kFALSE);
844 fFitterQuadrantPhi.AddAt(fitter,ifit);
845 }
846 fInit= kTRUE;
847}
848
849void AliTPCcalibUnlinearity::EvalFitters(){
850 //
851 //
852 // Evaluate fitters
853 //
854 TLinearFitter * fitter = 0;
855 TVectorD vec;
856 for (Int_t ifit=0; ifit<38; ifit++){
857 fitter=(TLinearFitter*)fFittersOutR.At(ifit);
858 if (fitter) {
859 fitter->Eval();
860 fitter->GetParameters(vec);
861 fParamsOutR.AddAt(vec.Clone(),ifit);
862 fitter->GetErrors(vec);
863 fErrorsOutR.AddAt(vec.Clone(),ifit);
864 }
865 fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
866 if (fitter) {
867 fitter->Eval();
868 fitter->GetParameters(vec);
869 fParamsOutZ.AddAt(vec.Clone(),ifit);
870 fitter->GetErrors(vec);
871 fErrorsOutZ.AddAt(vec.Clone(),ifit);
872 }
873 }
874}
875
876
877
878
879
880void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
881 //
882 //
883 //
884 // TTree * tree = chainUnlinP;
885 AliTPCcalibUnlinearity *calib = this;
886 //
887 AliTPCclusterMI * cl=0;
888 Double_t ty,tz;
889 TVectorD *vy=0, *vz=0;
890 TVectorD *vy2=0, *vz2=0;
891 Long64_t nentries = tree->GetEntries();
892 if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
893 //
894 {
895 for (Long64_t i=0; i<nentries; i++){
896 tree->SetBranchAddress("Cl.",&cl);
897 tree->SetBranchAddress("y1",&ty);
898 tree->SetBranchAddress("z1",&tz);
899 tree->SetBranchAddress("py1.",&vy);
900 tree->SetBranchAddress("pz1.",&vz);
901 tree->SetBranchAddress("py2.",&vy2);
902 tree->SetBranchAddress("pz2.",&vz2);
903 //
904 tree->GetEntry(i);
905 if (!cl) continue;
906 if (i%10000==0) printf("%d\n",(Int_t)i);
907 Int_t row= cl->GetRow();
908 if (cl->GetDetector()>36) row+=64;
909 calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
910 ty,tz,(*vy)[1],(*vz)[1],1);
911
912 if (cl->GetType()<0) continue;
913 Double_t dy = cl->GetY()-ty;
914 Double_t dz = cl->GetZ()-tz;
915 //
916 //
917 if (TMath::Abs(dy)>0.25) continue;
918 if (TMath::Abs(dz)>0.25) continue;
919
920 if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
921 if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
922 // Apply sagita cut
923 if (cl->GetType()>=0) {
924 calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
925 ty,tz,(*vy)[1],(*vz)[1],1);
926 }
927 }
928 }
929}
930
931
932void AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
933 //
934 // Make position corrections
935 // for the moment Only using debug streamer
936 // chainUnlinD - debug tree
937 // param - parameters to be updated
938 // maxPoints - maximal number of points using for fit
939 // verbose - print info flag
940 //
941 // Current implementation - only using debug streamers
942 //
943
944 /*
945 //Defaults
946 Int_t maxPoints=100000;
947 */
948 /*
949 Usage:
950 //0. Load libraries
951 gSystem->Load("libANALYSIS");
952 gSystem->Load("libSTAT");
953 gSystem->Load("libTPCcalib");
954
955
956 //1. Load Parameters to be modified:
957 //e.g:
958
959 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
960 AliCDBManager::Instance()->SetRun(0)
961 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
962 //
963 //2. Load the debug streamer
964 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
965 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
966 AliXRDPROOFtoolkit tool;
967 TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
968 chainUnlin->Lookup();
969 TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
970 chainUnlinD->Lookup();
971
972 //3. Do fits and store results
973 //
974 AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
975 TFile f("paramout.root","recreate");
976 param->Write("clusterParam");
977 f.Close();
978 //4. Verification
979 TFile f2("paramout.root");
980 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
981 param2->SetInstance(param2);
982 chainUnlinD->Draw("fitZ0:AliTPCClusterParam::SPosCorrection(1,0,Cl.fPad,Cl.fTimeBin,Cl.fZ,Cl.fSigmaY2,Cl.fSigmaZ2,Cl.fMax)","Cl.fDetector<36","",10000); // should be line
983
984 */
985
986
987 TStatToolkit toolkit;
988 Double_t chi2;
989 TVectorD fitParamY0;
990 TVectorD fitParamY1;
991 TVectorD fitParamZ0;
992 TVectorD fitParamZ1;
993 TMatrixD covMatrix;
994 Int_t npoints;
995
996 chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
997 chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
998 chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
999 chainUnlinD->SetAlias("st","(sin(dt)-dt)");
1000 chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
1001 //
1002 chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
1003 //
1004 //
1005 //
1006 TCut cutE("Cl.fType>=0");
1007 TCut cutDy("abs(Cl.fY-y1)<0.4");
1008 TCut cutDz("abs(Cl.fZ-z1)<0.4");
1009 TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
1010 TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
1011 TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
1012 //
1013 TCut cutA("Cl.fZ>0");
1014 TCut cutC("Cl.fZ<0");
1015
1016 TString fstringY="";
1017 //
1018 fstringY+="(dp)++"; //1
1019 fstringY+="(dp)*di++"; //2
1020 fstringY+="(sp)++"; //3
1021 fstringY+="(sp)*di++"; //4
1022 fstringY+="(dq)++"; //5
1023 TString fstringZ="";
1024 fstringZ+="(dt)++"; //1
1025 fstringZ+="(dt)*di++"; //2
1026 fstringZ+="(st)++"; //3
1027 fstringZ+="(st)*di++"; //4
1028 fstringZ+="(dq)++"; //5
1029 //
1030 // Z corrections
1031 //
1032 TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
1033 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1034 strZ0->Tokenize("++")->Print();
1035 param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
1036 //
1037 TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
1038 //
1039 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1040 strZ1->Tokenize("++")->Print();
1041 param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
1042 param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
1043 //
1044 // Y corrections
1045 //
1046 TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
1047 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1048 strY0->Tokenize("++")->Print();
1049 param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
1050
1051
1052 TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
1053 //
1054 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
1055 strY1->Tokenize("++")->Print();
1056 param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
1057 param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
1058 //
1059 //
1060 //
1061 chainUnlinD->SetAlias("fitZ0",strZ0->Data());
1062 chainUnlinD->SetAlias("fitZ1",strZ1->Data());
1063 chainUnlinD->SetAlias("fitY0",strY0->Data());
1064 chainUnlinD->SetAlias("fitY1",strY1->Data());
1065}
1066
1067
1068
1069
1070
1071
1072