]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibUnlinearity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moved old AliMagFCheb to AliMagF, small fixes/optimizations in field classes
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibUnlinearity.cxx
CommitLineData
c11fe047 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
28b5b7c8 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
c11fe047 36 //
37*/
38
39
40#include "TLinearFitter.h"
41
42#include "Riostream.h"
43#include "TChain.h"
44#include "TTree.h"
45#include "TH1F.h"
46#include "TH2F.h"
47#include "TH3F.h"
48#include "THnSparse.h"
49#include "TList.h"
50#include "TMath.h"
51#include "TCanvas.h"
52#include "TFile.h"
53#include "TF1.h"
54#include "TVectorD.h"
55#include "TProfile.h"
56#include "TGraphErrors.h"
57#include "TCanvas.h"
28b5b7c8 58#include "TCut.h"
59
c11fe047 60
61#include "AliTPCclusterMI.h"
62#include "AliTPCseed.h"
63#include "AliESDVertex.h"
64#include "AliESDEvent.h"
65#include "AliESDfriend.h"
66#include "AliESDInputHandler.h"
67#include "AliAnalysisManager.h"
68#include "AliTracker.h"
69#include "AliMagFMaps.h"
70#include "AliTPCCalROC.h"
71
72#include "AliLog.h"
73
74
75#include "TTreeStream.h"
76#include "AliTPCTracklet.h"
77#include "TTimeStamp.h"
78#include "AliTPCcalibDB.h"
79#include "AliTPCcalibLaser.h"
80#include "AliDCSSensorArray.h"
81#include "AliDCSSensor.h"
82#include "TLinearFitter.h"
28b5b7c8 83#include "AliTPCClusterParam.h"
84#include "TStatToolkit.h"
c11fe047 85
86
87#include "AliTPCcalibUnlinearity.h"
88
28b5b7c8 89AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::fgInstance = 0;
c11fe047 90
91ClassImp(AliTPCcalibUnlinearity)
92
93AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
94 AliTPCcalibBase(),
95 fDiffHistoLine(0),
96 fDiffHistoPar(0),
97 fFittersOutR(38),
28b5b7c8 98 fFittersOutZ(38),
99 fParamsOutR(38),
100 fParamsOutZ(38),
101 fErrorsOutR(38),
102 fErrorsOutZ(38)
c11fe047 103{
104 //
105 // Defualt constructor
106 //
107}
108
109AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
110 AliTPCcalibBase(name,title),
111 fDiffHistoLine(0),
112 fDiffHistoPar(0),
113 fFittersOutR(38),
28b5b7c8 114 fFittersOutZ(38),
115 fParamsOutR(38),
116 fParamsOutZ(38),
117 fErrorsOutR(38),
118 fErrorsOutZ(38)
c11fe047 119{
120 //
121 // Non default constructor
122 //
123 MakeFitters();
124}
125
126AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
127 //
128 //
129 //
130 if (fDiffHistoLine) delete fDiffHistoLine;
131 if (fDiffHistoPar) delete fDiffHistoPar;
132}
133
134
28b5b7c8 135AliTPCcalibUnlinearity* AliTPCcalibUnlinearity::Instance()
136{
137 //
138 // Singleton implementation
139 // Returns an instance of this class, it is created if neccessary
140 //
141 if (fgInstance == 0){
142 fgInstance = new AliTPCcalibUnlinearity();
143 }
144 return fgInstance;
145}
146
147
148
149
c11fe047 150void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
151 //
152 //
153 //
154 const Int_t kdrow=10;
155 const Int_t kMinCluster=35;
156 if (!fDiffHistoLine) MakeHisto();
157 if (!seed) return;
28b5b7c8 158 if (TMath::Abs(fMagF)>0.01) return; // use only data without field
c11fe047 159 Int_t counter[72];
160 for (Int_t i=0;i<72;i++) counter[i]=0;
161 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
162 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
163 if (!cluster) continue;
164 if (cluster->GetType()<0) continue;
165 counter[cluster->GetDetector()]++;
166 }
167 //
168 for (Int_t is=0; is<72;is++){
169 if (counter[is]>kMinCluster) ProcessDiff(seed, is);
170 }
171}
172
173
174void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
175 //
176 // process differnce of the cluster in respect with linear and parabolic fit
177 //
178 const Double_t kXmean=134;
179 const Int_t kdrow=10;
180 static TLinearFitter fy1(2,"pol1");
181 static TLinearFitter fy2(3,"pol2");
182 static TLinearFitter fz1(2,"pol1");
183 static TLinearFitter fz2(3,"pol2");
184 //
185 static TVectorD py1(2);
186 static TVectorD py2(3);
187 static TVectorD pz1(2);
188 static TVectorD pz2(3);
189 //
190 fy1.ClearPoints();
191 fy2.ClearPoints();
192 fz1.ClearPoints();
193 fz2.ClearPoints();
194 //
195 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
196 AliTPCclusterMI *c=track->GetClusterPointer(irow);
197 if (!c) continue;
198 if (c->GetDetector()!=isec) continue;
199 if (c->GetType()<0) continue;
28b5b7c8 200 Double_t dx = c->GetX()-kXmean;
c11fe047 201 Double_t x[2]={dx, dx*dx};
202 fy1.AddPoint(x,c->GetY());
203 fy2.AddPoint(x,c->GetY());
204 fz1.AddPoint(x,c->GetZ());
205 fz2.AddPoint(x,c->GetZ());
206 }
207 fy1.Eval(); fy1.GetParameters(py1);
208 fy2.Eval(); fy2.GetParameters(py2);
209 fz1.Eval(); fz1.GetParameters(pz1);
210 fz2.Eval(); fz2.GetParameters(pz2);
211 //
212
213 for (Int_t irow=0;irow<159;irow++) {
214 AliTPCclusterMI *c=track->GetClusterPointer(irow);
215 if (!c) continue;
216 if (c->GetDetector()!=isec) continue;
217 if (c->GetType()<0) continue;
28b5b7c8 218 Double_t dx = c->GetX()-kXmean;
219 Double_t y1 = py1[0]+py1[1]*dx;
220 Double_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
c11fe047 221 //
28b5b7c8 222 Double_t z1 = pz1[0]+pz1[1]*dx;
223 Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
c11fe047 224 //
225 //
226 Double_t x[10];
227 x[0]=isec;
228 x[1]=irow;
229 x[2]=c->GetY();
230 x[3]=c->GetZ();
231 x[4]=py1[1];
232 x[5]=pz1[1];
233 x[6]=py2[2]*150*150/4; // sagita 150 cm
234 //
235 x[7]=c->GetY()-y1;
236 x[8]=c->GetZ()-z1;
237 fDiffHistoLine->Fill(x);
238 x[7]=c->GetY()-y2;
239 x[8]=c->GetZ()-z2;
240 fDiffHistoPar->Fill(x);
241
28b5b7c8 242 if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(pz2[2]*150*150/4)<1){
c11fe047 243 // Apply sagita cut
244 AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
245 py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
246 }
247
248 if (fStreamLevel>1){
249 TTreeSRedirector *cstream = GetDebugStreamer();
250 if (cstream){
251 (*cstream)<<"Diff"<<
28b5b7c8 252 "run="<<fRun<< // run number
253 "event="<<fEvent<< // event number
254 "time="<<fTime<< // time stamp of event
255 "trigger="<<fTrigger<< // trigger
256 "mag="<<fMagF<< // magnetic field
c11fe047 257 "isec="<<isec<<
258 "Cl.="<<c<<
259 "y1="<<y1<<
260 "y2="<<y2<<
261 "z1="<<z1<<
262 "z2="<<z2<<
263 "py1.="<<&py1<<
264 "py2.="<<&py2<<
265 "pz1.="<<&pz1<<
266 "pz2.="<<&pz2<<
267 "\n";
268 }
269 }
270 }
271}
272
273
274void AliTPCcalibUnlinearity::MakeHisto(){
275 //
276 //
277 //
278 TString axisName[10];
279 Double_t xmin[10], xmax[10];
280 Int_t nbins[10];
281 //
282 //
283 axisName[0]="sector";
284 xmin[0]=0; xmax[0]=72; nbins[0]=72;
285 //
286 axisName[1]="row";
287 xmin[1]=0; xmax[1]=159; nbins[1]=159;
288 //
289 axisName[2]="ly";
290 xmin[2]=-50; xmax[2]=50; nbins[2]=10;
291 //
292 axisName[3]="lz";
293 xmin[3]=-250; xmax[3]=250; nbins[3]=50;
294 //
295 axisName[4]="p2";
296 xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
297 //
298 axisName[5]="p3";
299 xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
300 //
301 axisName[6]="p4";
302 xmin[6]=-2; xmax[6]=2; nbins[6]=20;
303 //
304 axisName[7]="dy";
305 xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
306 //
307 axisName[8]="dz";
308 xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
309 //
310 fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
311 fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
312 for (Int_t i=0; i<9;i++){
313 fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
314 fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
315 fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
316 fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
317 }
318}
319
320
321void AliTPCcalibUnlinearity::Terminate(){
322 //
323 // Terminate function
324 // call base terminate + Eval of fitters
325 //
326 Info("AliTPCcalibUnlinearities","Terminate");
327 EvalFitters();
328 DumpTree();
329 AliTPCcalibBase::Terminate();
330}
331
332
333void AliTPCcalibUnlinearity::DumpTree(){
334 //
335 //
336 //
337 if (fStreamLevel==0) return;
338 TTreeSRedirector *cstream = GetDebugStreamer();
339 if (!cstream) return;
340 //
341 THnSparse *his=0;
342 Double_t position[10];
343 Double_t value;
344 Int_t *bins = new Int_t[10];
345 //
346 for (Int_t ihis=0; ihis<2; ihis++){
347 his = (ihis==0)? fDiffHistoLine:fDiffHistoPar;
348 if (!his) continue;
349 //
350 Int_t ndim = his->GetNdimensions();
351 //
352 for (Long64_t i = 0; i < his->GetNbins(); ++i) {
353 value = his->GetBinContent(i, bins);
354 for (Int_t idim = 0; idim < ndim; idim++) {
355 position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
356 }
357 (*cstream)<<"Resol"<<
358 "ihis="<<ihis<<
359 "bincont="<<value;
360 for (Int_t idim=0;idim<ndim;idim++){
361 (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
362 }
363 (*cstream)<<"Resol"<<
364 "\n";
365 }
366 }
367}
368
369
28b5b7c8 370void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Double_t dz, Double_t dy, Double_t p2, Double_t p3, Double_t dr, Int_t npoints){
c11fe047 371 //
372 //
373 // Simple distortion fit in outer sectors
374 //
375 // Model - 2 exponential decay toward the center of chamber
28b5b7c8 376 // - Decay length 10 and 5 cm
c11fe047 377 // - Decay amplitude - Z dependent
378 //
379 /*
380 chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
381 chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
382 chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
383 chainUnlin->SetAlias("dzt","dz*side");
384 chainUnlin->SetAlias("dr","(1-abs(lz/250))");
385 chainUnlin->SetAlias("dout","(159-row)*1.5");
386 chainUnlin->SetAlias("din","row*0.75");
387 chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
28b5b7c8 388 chainUnlin->SetAlias("eout5","exp(-(dout)/5.)");
c11fe047 389 */
28b5b7c8 390
391 Double_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1
392 Double_t dzt = dz*side;
393 Double_t dout = (159-row)*1.5; // distance to the last pad row - (valid only for long pads)
394 if (dout>30) return; // use only edge pad rows
395 dr-=0.5; // dr shifted to the middle - reduce numerical instabilities
396
397 Double_t eout10 = TMath::Exp(-dout/10.);
398 Double_t eout5 = TMath::Exp(-dout/5.);
399 //
400 Double_t eoutp = (eout10+eout5)*0.5;
401 Double_t eoutm = (eout10-eout5)*0.5;
402
c11fe047 403 //
404 Double_t xxxR[6], xxxZ[6], xxxRZ[6];
405 //TString fstring="";
28b5b7c8 406 xxxZ[0]=eoutp; //fstring+="eoutp++";
407 xxxZ[1]=eoutm; //fstring+="eoutm++";
408 xxxZ[2]=dr*eoutp; //fstring+="dr*eoutp++";
409 xxxZ[3]=dr*eoutm; //fstring+="dr*eoutm++";
410 xxxZ[4]=dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
411 xxxZ[5]=dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
412 //
413 xxxR[0]=p2*eoutp; //fstring+="eoutp++";
414 xxxR[1]=p2*eoutm; //fstring+="eoutm++";
415 xxxR[2]=p2*dr*eoutp; //fstring+="dr*eoutp++";
416 xxxR[3]=p2*dr*eoutm; //fstring+="dr*eoutm++";
417 xxxR[4]=p2*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
418 xxxR[5]=p2*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
419 //
420 xxxRZ[0]=p3*eoutp; //fstring+="eoutp++";
421 xxxRZ[1]=p3*eoutm; //fstring+="eoutm++";
422 xxxRZ[2]=p3*dr*eoutp; //fstring+="dr*eoutp++";
423 xxxRZ[3]=p3*dr*eoutm; //fstring+="dr*eoutm++";
424 xxxRZ[4]=p3*dr*dr*eoutp; //fstring+="dr*dr*eoutp++";
425 xxxRZ[5]=p3*dr*dr*eoutm; //fstring+="dr*dr*eoutm++";
c11fe047 426 //
427 TLinearFitter * fitter=0;
28b5b7c8 428 Double_t err=0.1;
c11fe047 429 for (Int_t ip=0; ip<npoints; ip++){
430 //
431 fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
28b5b7c8 432 fitter->AddPoint(xxxR,dy,err);
433 //fitter->AddPoint(xxxRZ,dz);
c11fe047 434 fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
28b5b7c8 435 fitter->AddPoint(xxxZ,dzt,err);
c11fe047 436 //
437 if (side>0){
438 fitter = (TLinearFitter*)fFittersOutR.At(36);
28b5b7c8 439 fitter->AddPoint(xxxR,dy,err);
440 //fitter->AddPoint(xxxRZ,dz);
c11fe047 441 fitter = (TLinearFitter*)fFittersOutZ.At(36);
28b5b7c8 442 fitter->AddPoint(xxxZ,dzt,err);
c11fe047 443 }
444 if (side<0){
28b5b7c8 445 fitter = (TLinearFitter*)fFittersOutR.At(37);
446 fitter->AddPoint(xxxR,dy,err);
447 //fitter->AddPoint(xxxRZ,dz);
448 fitter = (TLinearFitter*)fFittersOutZ.At(37);
449 fitter->AddPoint(xxxZ,dzt,err);
c11fe047 450 }
451 }
452}
453
454
455void AliTPCcalibUnlinearity::MakeFitters(){
456 //
457 //
458 //
459 // Make outer fitters
460 TLinearFitter * fitter = 0;
461 for (Int_t ifit=0; ifit<38; ifit++){
462 fitter = new TLinearFitter(7,"hyp6");
463 fitter->StoreData(kFALSE);
464 fFittersOutR.AddAt(fitter,ifit);
465 fitter = new TLinearFitter(7,"hyp6");
466 fitter->StoreData(kFALSE);
467 fFittersOutZ.AddAt(fitter,ifit);
468 }
469}
470
471void AliTPCcalibUnlinearity::EvalFitters(){
472 //
473 //
474 // Evaluate fitters
475 //
476 TLinearFitter * fitter = 0;
28b5b7c8 477 TVectorD vec;
c11fe047 478 for (Int_t ifit=0; ifit<38; ifit++){
479 fitter=(TLinearFitter*)fFittersOutR.At(ifit);
28b5b7c8 480 if (fitter) {
481 fitter->Eval();
482 fitter->GetParameters(vec);
483 fParamsOutR.AddAt(vec.Clone(),ifit);
484 fitter->GetErrors(vec);
485 fErrorsOutR.AddAt(vec.Clone(),ifit);
486 }
c11fe047 487 fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
28b5b7c8 488 if (fitter) {
489 fitter->Eval();
490 fitter->GetParameters(vec);
491 fParamsOutZ.AddAt(vec.Clone(),ifit);
492 fitter->GetErrors(vec);
493 fErrorsOutZ.AddAt(vec.Clone(),ifit);
494 }
c11fe047 495 }
496}
497
28b5b7c8 498Double_t AliTPCcalibUnlinearity::GetDr(Int_t sector, Double_t dout, Double_t dr){
499 //
500 //
501 //
502 TVectorD * vec = GetParamOutR(sector);
503 if (!vec) return 0;
504 dr-=0.5; // dr=0 at the middle drift length
505 Double_t eout10 = TMath::Exp(-dout/10.);
506 Double_t eout5 = TMath::Exp(-dout/5.);
507 Double_t eoutp = (eout10+eout5)*0.5;
508 Double_t eoutm = (eout10-eout5)*0.5;
509
510 Double_t result=0;
511 result+=(*vec)[1]*eoutp;
512 result+=(*vec)[2]*eoutm;
513 result+=(*vec)[3]*eoutp*dr;
514 result+=(*vec)[4]*eoutm*dr;
515 result+=(*vec)[5]*eoutp*dr*dr;
516 result+=(*vec)[6]*eoutm*dr*dr;
517 return result;
518}
519
520
521Double_t AliTPCcalibUnlinearity::GetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
522 //
523 //
524 //
525 Double_t alpha = TMath::ATan2(gy,gx);
526 if (alpha<0) alpha+=TMath::Pi()*2;
527 Int_t lsector = Int_t(18*alpha/(TMath::Pi()*2));
528 alpha = (lsector+0.5)*TMath::Pi()/9.;
529 //
530 Double_t r[3];
531 r[0]=gx; r[1]=gy;
532 Double_t cs=TMath::Cos(-alpha), sn=TMath::Sin(-alpha), x=r[0];
533 r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs;
534 //
535 AliTPCROC * roc = AliTPCROC::Instance();
536 Double_t xout = roc->GetPadRowRadiiUp(roc->GetNRows(37)-1);
537 Double_t dout = xout-r[0];
538 if (dout<0) return 0;
539 Double_t dr = 1-TMath::Abs(gz/250.);
540 if (gz<0) lsector+=18;
541 if (stype==0) lsector = (gz>0) ? 36:37;
542 if (stype<0) return lsector; //
543 return GetDr(lsector,dout,dr);
544}
545
c11fe047 546
28b5b7c8 547Double_t AliTPCcalibUnlinearity::GetDz(Int_t sector, Double_t dout, Double_t dr){
548 //
549 //
550 //
551 TVectorD * vec = GetParamOutZ(sector);
552 if (!vec) return 0;
553 dr-=0.5; // 0 at the middle
554 Double_t eout10 = TMath::Exp(-dout/10.);
555 Double_t eout5 = TMath::Exp(-dout/5.);
556 Double_t eoutp = (eout10+eout5)*0.5;
557 Double_t eoutm = (eout10-eout5)*0.5;
558
559
560 Double_t result=0;
561 result+=(*vec)[1]*eoutp;
562 result+=(*vec)[2]*eoutm;
563 result+=(*vec)[3]*eoutp*dr;
564 result+=(*vec)[4]*eoutm*dr;
565 result+=(*vec)[5]*eoutp*dr*dr;
566 result+=(*vec)[6]*eoutm*dr*dr;
567 return result;
568}
569
570
571Double_t AliTPCcalibUnlinearity::SGetDr(Int_t sector, Double_t dout, Double_t dr){
572 //
573 //
574 //
575 return fgInstance->GetDr(sector,dout,dr);
576}
577Double_t AliTPCcalibUnlinearity::SGetGDr(Int_t stype,Float_t gx, Float_t gy,Float_t gz){
578 //
579 //
580 //
581 return fgInstance->GetGDr(stype,gx,gy,gz);
582}
583Double_t AliTPCcalibUnlinearity::SGetDz(Int_t sector, Double_t dout, Double_t dr){
584 //
585 //
586 //
587 return fgInstance->GetDz(sector,dout,dr);
588}
589
590
591
592
593void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
c11fe047 594 //
595 //
596 //
597 // TTree * tree = chainUnlinP;
598 AliTPCcalibUnlinearity *calib = this;
599 //
600 AliTPCclusterMI * cl=0;
28b5b7c8 601 Double_t ty,tz;
c11fe047 602 TVectorD *vy=0, *vz=0;
28b5b7c8 603 TVectorD *vy2=0, *vz2=0;
604 Long64_t nentries = tree->GetEntries();
605 if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints);
c11fe047 606 //
607 {
28b5b7c8 608 for (Long64_t i=0; i<nentries; i++){
c11fe047 609 tree->SetBranchAddress("Cl.",&cl);
610 tree->SetBranchAddress("y1",&ty);
611 tree->SetBranchAddress("z1",&tz);
612 tree->SetBranchAddress("py1.",&vy);
613 tree->SetBranchAddress("pz1.",&vz);
28b5b7c8 614 tree->SetBranchAddress("py2.",&vy2);
615 tree->SetBranchAddress("pz2.",&vz2);
c11fe047 616 //
617 tree->GetEntry(i);
618 if (!cl) continue;
28b5b7c8 619 if (i%10000==0) printf("%d\n",(Int_t)i);
620 Int_t row= cl->GetRow();
621 if (cl->GetDetector()>36) row+=64;
622 if (cl->GetType()<0) continue;
623 Double_t dy = cl->GetY()-ty;
624 Double_t dz = cl->GetZ()-tz;
625 Double_t dr = 1-TMath::Abs(cl->GetZ()/250);
626 //
627 //
628 if (TMath::Abs(dy)>0.25) continue;
629 if (TMath::Abs(dz)>0.25) continue;
630
631 if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
632 if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
633 // Apply sagita cut
634 calib->AddPoint(cl->GetDetector(), row, dz, dy,
635 (*vy)[1],(*vz)[1], dr, 1);
c11fe047 636 }
637 }
638}
28b5b7c8 639
640
641void AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
642 //
643 // Make position corrections
644 // for the moment Only using debug streamer
645 // chainUnlinD - debug tree
646 // param - parameters to be updated
647 // maxPoints - maximal number of points using for fit
648 // verbose - print info flag
649 //
650 // Current implementation - only using debug streamers
651 //
652
653 /*
654 //Defaults
655 Int_t maxPoints=100000;
656 */
657 /*
658 Usage:
659 //0. Load libraries
660 gSystem->Load("libANALYSIS");
661 gSystem->Load("libSTAT");
662 gSystem->Load("libTPCcalib");
663
664
665 //1. Load Parameters to be modified:
666 //e.g:
667
668 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
669 AliCDBManager::Instance()->SetRun(0)
670 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
671 //
672 //2. Load the debug streamer
673 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
674 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
675 AliXRDPROOFtoolkit tool;
676 TChain * chainUnlin = tool.MakeChain("unlin.txt","Resol",0,10200);
677 chainUnlin->Lookup();
678 TChain * chainUnlinD = tool.MakeChain("unlin.txt","Diff",0,10200);
679 chainUnlinD->Lookup();
680
681 //3. Do fits and store results
682 //
683 AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
684 TFile f("paramout.root","recreate");
685 param->Write("clusterParam");
686 f.Close();
687 //4. Verification
688 TFile f2("paramout.root");
689 AliTPCClusterParam *param2 = (AliTPCClusterParam*)f2.Get("clusterParam");
690 param2->SetInstance(param2);
691 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
692
693 */
694
695
696 TStatToolkit toolkit;
697 Double_t chi2;
698 TVectorD fitParamY0;
699 TVectorD fitParamY1;
700 TVectorD fitParamZ0;
701 TVectorD fitParamZ1;
702 TMatrixD covMatrix;
703 Int_t npoints;
704
705 chainUnlinD->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
706 chainUnlinD->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
707 chainUnlinD->SetAlias("sp","(sin(dp*pi)-dp*pi)");
708 chainUnlinD->SetAlias("st","(sin(dt)-dt)");
709 chainUnlinD->SetAlias("dq","(-1+(Cl.fZ>0)*2)/Cl.fMax");
710 //
711 chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
712 //
713 //
714 //
715 TCut cutE("Cl.fType>=0");
716 TCut cutDy("abs(Cl.fY-y1)<0.4");
717 TCut cutDz("abs(Cl.fZ-z1)<0.4");
718 TCut cutY2("abs(py2.fElements[2])*150^2/4<1");
719 TCut cutZ2("abs(pz2.fElements[2])*150^2/4<1");
720 TCut cutAll = cutE+cutDy+cutDz+cutY2+cutZ2;
721 //
722 TCut cutA("Cl.fZ>0");
723 TCut cutC("Cl.fZ<0");
724
725 TString fstringY="";
726 //
727 fstringY+="(dp)++"; //1
728 fstringY+="(dp)*di++"; //2
729 fstringY+="(sp)++"; //3
730 fstringY+="(sp)*di++"; //4
731 fstringY+="(dq)++"; //5
732 TString fstringZ="";
733 fstringZ+="(dt)++"; //1
734 fstringZ+="(dt)*di++"; //2
735 fstringZ+="(st)++"; //3
736 fstringZ+="(st)*di++"; //4
737 fstringZ+="(dq)++"; //5
738 //
739 // Z corrections
740 //
741 TString *strZ0 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamZ0,covMatrix,-1,0,maxPoints);
742 printf("Z0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
743 strZ0->Tokenize("++")->Print();
744 param->fPosZcor[0] = (TVectorD*) fitParamZ0.Clone();
745 //
746 TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
747 //
748 printf("Z1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
749 strZ1->Tokenize("++")->Print();
750 param->fPosZcor[1] = (TVectorD*) fitParamZ1.Clone();
751 param->fPosZcor[2] = (TVectorD*) fitParamZ1.Clone();
752 //
753 // Y corrections
754 //
755 TString *strY0 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector<36"+cutAll, chi2,npoints,fitParamY0,covMatrix,-1,0,maxPoints);
756 printf("Y0 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
757 strY0->Tokenize("++")->Print();
758 param->fPosYcor[0] = (TVectorD*) fitParamY0.Clone();
759
760
761 TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
762 //
763 printf("Y1 - chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
764 strY1->Tokenize("++")->Print();
765 param->fPosYcor[1] = (TVectorD*) fitParamY1.Clone();
766 param->fPosYcor[2] = (TVectorD*) fitParamY1.Clone();
767 //
768 //
769 //
770 chainUnlinD->SetAlias("fitZ0",strZ0->Data());
771 chainUnlinD->SetAlias("fitZ1",strZ1->Data());
772 chainUnlinD->SetAlias("fitY0",strY0->Data());
773 chainUnlinD->SetAlias("fitY1",strY1->Data());
774}
775
776
777
778
779
780
781