1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
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
27 Unlinearities fitting:
28 Unlinearities at the edge aproximated using two exponential decays.
31 dz = dz0(r,z) +dr(r,z)*tan(theta)
32 dy = +dr(r,z)*tan(phi)
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");
46 #include "TLinearFitter.h"
48 #include "Riostream.h"
54 #include "THnSparse.h"
62 #include "TGraphErrors.h"
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"
76 #include "AliTPCCalROC.h"
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"
92 #include "AliTPCcalibUnlinearity.h"
93 #include "AliTPCPointCorrection.h"
96 ClassImp(AliTPCcalibUnlinearity)
98 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
113 fFitterQuadrantY(16*38),
114 fFitterQuadrantPhi(16*38)
117 // Defualt constructor
121 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
122 AliTPCcalibBase(name,title),
136 fFitterQuadrantY(16*38),
137 fFitterQuadrantPhi(16*38)
140 // Non default constructor
145 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
149 if (fDiffHistoLineZ) delete fDiffHistoLineY;
150 if (fDiffHistoParZ) delete fDiffHistoParY;
151 if (fDiffHistoLineZ) delete fDiffHistoLineZ;
152 if (fDiffHistoParZ) delete fDiffHistoParZ;
159 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
164 Init(); // work around for PROOF - initialize firs time used
166 const Int_t kdrow=10;
167 const Int_t kMinCluster=35;
169 if (TMath::Abs(fMagF)>0.01) return; // use only data without field
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()]++;
179 for (Int_t is=0; is<72;is++){
180 if (counter[is]>kMinCluster) ProcessDiff(seed, is);
181 if (counter[is]>kMinCluster) {
188 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
190 // process differnce of the cluster in respect with linear and parabolic fit
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");
200 static TVectorD py1(2);
201 static TVectorD py2(3);
202 static TVectorD pz1(2);
203 static TVectorD pz2(3);
210 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
211 AliTPCclusterMI *c=track->GetClusterPointer(irow);
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());
222 fy1.Eval(); fy1.GetParameters(py1);
223 fy2.Eval(); fy2.GetParameters(py2);
224 fz1.Eval(); fz1.GetParameters(pz1);
225 fz2.Eval(); fz2.GetParameters(pz2);
228 for (Int_t irow=0;irow<159;irow++) {
229 AliTPCclusterMI *c=track->GetClusterPointer(irow);
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;
236 Double_t z1 = pz1[0]+pz1[1]*dx;
237 Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
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;
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);
258 TTreeSRedirector *cstream = GetDebugStreamer();
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
270 "crR="<<corrR<< // radial correction
271 "cclY="<<corrclY<< // edge effect correction cl
272 "ctrY="<<corrtrY<< // edge effect correction using track
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){
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
299 "crR="<<corrR<< // radial correction
300 "cclY="<<corrclY<< // edge effect correction cl
301 "ctrY="<<corrtrY<< // edge effect correction using track
314 if (c->GetType()<0) continue; // don't use edge rphi cluster
318 x[0]=c->GetDetector();
320 x[2]=c->GetY()/c->GetX();
325 if (TMath::Abs(py2[2]*150*150/4)<kSagitaCut &&
326 TMath::Abs(pz2[2]*150*150/4)<kSagitaCut){
330 fDiffHistoLineY->Fill(x);
332 //fDiffHistoParY->Fill(x);
336 fDiffHistoLineZ->Fill(x);
338 //fDiffHistoParZ->Fill(x);
340 AddPoint(isec,c->GetX(),c->GetY(),c->GetZ(),y1,z1,py1[1],pz1[1],1);
347 void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
349 // fit the tracklet in OROC sepatately in Quadrant
353 const Int_t kMinClusterF=35;
354 const Int_t kMinClusterQ=10;
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
364 AliTPCROC * roc = AliTPCROC::Instance();
365 AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
366 const Double_t kXmean=roc->GetPadRowRadii(isec,53);
368 // full track fit parameters
370 TLinearFitter fyf(2,"pol1");
371 TLinearFitter fzf(2,"pol1");
372 TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
375 // make full fit as reference
377 for (Int_t iter=0; iter<2; iter++){
379 for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
380 AliTPCclusterMI *c=track->GetClusterPointer(irow);
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;
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;
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;
397 fyf.AddPoint(x,c->GetY(),0.1);
398 fzf.AddPoint(x,c->GetZ(),0.1);
400 nf = fyf.GetNpoints();
401 if (nf<kMinClusterF) return; // not enough points - skip
403 fyf.GetParameters(pyf);
406 fzf.GetParameters(pzf);
410 // Make Fitters and params for 3 fitters
412 TLinearFitter *fitters[4];
417 for (Int_t i=0;i<4;i++) {
418 fitters[i] = new TLinearFitter(2,"pol1");
420 params[i].ResizeTo(2);
421 errors[i].ResizeTo(2);
426 for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
427 AliTPCclusterMI *c=track->GetClusterPointer(irow);
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;
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;
442 if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
443 if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
446 if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
447 if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
451 for (Int_t i=0;i<4;i++) {
452 npoints[i] = fitters[i]->GetNpoints();
453 if (npoints[i]>=kMinClusterQ){
455 Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
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;
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;
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]);
483 Int_t sector = isec-36;
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);
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);
500 // Dump debug information
503 TTreeSRedirector *cstream = GetDebugStreamer();
504 Int_t sector = isec-36;
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;
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
520 "nf="<<nf<< // total number of points
521 "pyf.="<<&pyf<< // full OROC fit y
522 "pzf.="<<&pzf<< // full OROC fit z
524 "i0="<<i0<< // quadrant number
526 "n0="<<npoints[i0]<< // number of points
528 "p0.="<<¶ms[i0]<< // parameters
529 "p1.="<<¶ms[i1]<<
530 "e0.="<<&errors[i0]<< // errors
531 "e1.="<<&errors[i1]<<
532 "chi0="<<chi2C[i0]<< // chi2s
546 void AliTPCcalibUnlinearity::MakeHisto(){
550 TString axisName[10];
551 Double_t xmin[10], xmax[10];
555 axisName[0]="sector";
556 xmin[0]=0; xmax[0]=72; nbins[0]=72;
559 xmin[1]=0; xmax[1]=96; nbins[1]=96;
561 axisName[2]="tphi"; // tan phi - ly/lx
562 xmin[2]=-0.17; xmax[2]=0.17; nbins[2]=5;
564 axisName[3]="lz"; // global z
565 xmin[3]=-250; xmax[3]=250; nbins[3]=10;
567 axisName[4]="k"; // tangent - in angle
568 xmin[4]=-0.5; xmax[4]=0.5; nbins[4]=10;
571 axisName[5]="delta"; // delta
572 xmin[5]=-0.5; xmax[5]=0.5; nbins[5]=100;
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());
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);
608 void AliTPCcalibUnlinearity::Terminate(){
610 // Terminate function
611 // call base terminate + Eval of fitters
613 Info("AliTPCcalibUnlinearities","Terminate");
616 AliTPCcalibBase::Terminate();
619 Long64_t AliTPCcalibUnlinearity::Merge(TCollection *list) {
623 TIterator* iter = list->MakeIterator();
624 AliTPCcalibUnlinearity* cal = 0;
626 while ((cal = (AliTPCcalibUnlinearity*)iter->Next())) {
627 if (!cal->InheritsFrom(AliTPCcalibUnlinearity::Class())) {
635 void AliTPCcalibUnlinearity::Add(AliTPCcalibUnlinearity * calib){
640 if (fDiffHistoLineY && calib->fDiffHistoLineY){
641 fDiffHistoLineY->Add(calib->fDiffHistoLineY);
643 if (fDiffHistoParY && calib->fDiffHistoParY){
644 fDiffHistoParY->Add(calib->fDiffHistoParY);
646 if (fDiffHistoLineZ && calib->fDiffHistoLineZ){
647 fDiffHistoLineZ->Add(calib->fDiffHistoLineZ);
649 if (fDiffHistoParZ && calib->fDiffHistoParZ){
650 fDiffHistoParZ->Add(calib->fDiffHistoParZ);
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);
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);
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);
685 void AliTPCcalibUnlinearity::DumpTree(const char *fname){
689 TTreeSRedirector *cstream =new TTreeSRedirector(fname);
690 if (!cstream) return;
693 Double_t position[10];
695 Int_t *bins = new Int_t[10];
697 for (Int_t ihis=0; ihis<2; ihis++){
698 his = (ihis==0)? fDiffHistoLineY:fDiffHistoParY;
701 Int_t ndim = his->GetNdimensions();
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]);
708 (*cstream)<<"Resol"<<
711 for (Int_t idim=0;idim<ndim;idim++){
712 (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
714 (*cstream)<<"Resol"<<
722 void 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){
725 // Simple distortion fit in outer sectors
727 // Model - 2 exponential decay toward the center of chamber
728 // - Decay length 10 and 5 cm
729 // - Decay amplitude - Z dependent
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;
740 Double_t dr = 0.5 - TMath::Abs(cz/250.);
742 Double_t eout10 = TMath::Exp(-dout/10.);
743 Double_t eout5 = TMath::Exp(-dout/5.);
745 Double_t eoutp = (eout10+eout5)*0.5;
746 Double_t eoutm = (eout10-eout5)*0.5;
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++";
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++";
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++";
772 TLinearFitter * fitter=0;
774 for (Int_t ip=0; ip<npoints; ip++){
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);
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);
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);
798 void 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){
802 Float_t kMaxDz = 0.5; // cut on z in cm
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.);
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);
822 void AliTPCcalibUnlinearity::Init(){
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);
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);
849 void AliTPCcalibUnlinearity::EvalFitters(){
854 TLinearFitter * fitter = 0;
856 for (Int_t ifit=0; ifit<38; ifit++){
857 fitter=(TLinearFitter*)fFittersOutR.At(ifit);
860 fitter->GetParameters(vec);
861 fParamsOutR.AddAt(vec.Clone(),ifit);
862 fitter->GetErrors(vec);
863 fErrorsOutR.AddAt(vec.Clone(),ifit);
865 fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
868 fitter->GetParameters(vec);
869 fParamsOutZ.AddAt(vec.Clone(),ifit);
870 fitter->GetErrors(vec);
871 fErrorsOutZ.AddAt(vec.Clone(),ifit);
880 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){
884 // TTree * tree = chainUnlinP;
885 AliTPCcalibUnlinearity *calib = this;
887 AliTPCclusterMI * cl=0;
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);
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);
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);
912 if (cl->GetType()<0) continue;
913 Double_t dy = cl->GetY()-ty;
914 Double_t dz = cl->GetZ()-tz;
917 if (TMath::Abs(dy)>0.25) continue;
918 if (TMath::Abs(dz)>0.25) continue;
920 if (TMath::Abs((*vy2)[2]*150*150/4)>1) continue;
921 if (TMath::Abs((*vz2)[2]*150*150/4)>1) continue;
923 if (cl->GetType()>=0) {
924 calib->AddPoint(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(),
925 ty,tz,(*vy)[1],(*vz)[1],1);
932 void AliTPCcalibUnlinearity::MakeQPosNormAll(TTree * chainUnlinD, AliTPCClusterParam * param, Int_t maxPoints){
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
941 // Current implementation - only using debug streamers
946 Int_t maxPoints=100000;
951 gSystem->Load("libANALYSIS");
952 gSystem->Load("libSTAT");
953 gSystem->Load("libTPCcalib");
956 //1. Load Parameters to be modified:
959 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
960 AliCDBManager::Instance()->SetRun(0)
961 AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
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();
972 //3. Do fits and store results
974 AliTPCcalibUnlinearity::MakeQPosNormAll(chainUnlinD,param,200000,0) ;
975 TFile f("paramout.root","recreate");
976 param->Write("clusterParam");
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
987 TStatToolkit toolkit;
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");
1002 chainUnlinD->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
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;
1013 TCut cutA("Cl.fZ>0");
1014 TCut cutC("Cl.fZ<0");
1016 TString fstringY="";
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
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();
1037 TString *strZ1 = toolkit.FitPlane(chainUnlinD,"(Cl.fZ-z2):1",fstringZ.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamZ1,covMatrix,-1,0,maxPoints);
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();
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();
1052 TString *strY1 = toolkit.FitPlane(chainUnlinD,"(Cl.fY-y2):1",fstringY.Data(), "Cl.fDetector>36"+cutAll, chi2,npoints,fitParamY1,covMatrix,-1,0,maxPoints);
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();
1061 chainUnlinD->SetAlias("fitZ0",strZ0->Data());
1062 chainUnlinD->SetAlias("fitZ1",strZ1->Data());
1063 chainUnlinD->SetAlias("fitY0",strY0->Data());
1064 chainUnlinD->SetAlias("fitY1",strY1->Data());