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 Only Linear and parabolic fit used
19 To be used for laser or for data without field
24 #include "TLinearFitter.h"
26 #include "Riostream.h"
32 #include "THnSparse.h"
40 #include "TGraphErrors.h"
43 #include "AliTPCclusterMI.h"
44 #include "AliTPCseed.h"
45 #include "AliESDVertex.h"
46 #include "AliESDEvent.h"
47 #include "AliESDfriend.h"
48 #include "AliESDInputHandler.h"
49 #include "AliAnalysisManager.h"
50 #include "AliTracker.h"
51 #include "AliMagFMaps.h"
52 #include "AliTPCCalROC.h"
57 #include "TTreeStream.h"
58 #include "AliTPCTracklet.h"
59 #include "TTimeStamp.h"
60 #include "AliTPCcalibDB.h"
61 #include "AliTPCcalibLaser.h"
62 #include "AliDCSSensorArray.h"
63 #include "AliDCSSensor.h"
64 #include "TLinearFitter.h"
67 #include "AliTPCcalibUnlinearity.h"
70 ClassImp(AliTPCcalibUnlinearity)
72 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
80 // Defualt constructor
84 AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
85 AliTPCcalibBase(name,title),
92 // Non default constructor
97 AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
101 if (fDiffHistoLine) delete fDiffHistoLine;
102 if (fDiffHistoPar) delete fDiffHistoPar;
106 void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
110 const Int_t kdrow=10;
111 const Int_t kMinCluster=35;
112 if (!fDiffHistoLine) MakeHisto();
115 for (Int_t i=0;i<72;i++) counter[i]=0;
116 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
117 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
118 if (!cluster) continue;
119 if (cluster->GetType()<0) continue;
120 counter[cluster->GetDetector()]++;
123 for (Int_t is=0; is<72;is++){
124 if (counter[is]>kMinCluster) ProcessDiff(seed, is);
129 void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
131 // process differnce of the cluster in respect with linear and parabolic fit
133 const Double_t kXmean=134;
134 const Int_t kdrow=10;
135 static TLinearFitter fy1(2,"pol1");
136 static TLinearFitter fy2(3,"pol2");
137 static TLinearFitter fz1(2,"pol1");
138 static TLinearFitter fz2(3,"pol2");
140 static TVectorD py1(2);
141 static TVectorD py2(3);
142 static TVectorD pz1(2);
143 static TVectorD pz2(3);
150 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
151 AliTPCclusterMI *c=track->GetClusterPointer(irow);
153 if (c->GetDetector()!=isec) continue;
154 if (c->GetType()<0) continue;
155 Float_t dx = c->GetX()-kXmean;
156 Double_t x[2]={dx, dx*dx};
157 fy1.AddPoint(x,c->GetY());
158 fy2.AddPoint(x,c->GetY());
159 fz1.AddPoint(x,c->GetZ());
160 fz2.AddPoint(x,c->GetZ());
162 fy1.Eval(); fy1.GetParameters(py1);
163 fy2.Eval(); fy2.GetParameters(py2);
164 fz1.Eval(); fz1.GetParameters(pz1);
165 fz2.Eval(); fz2.GetParameters(pz2);
168 for (Int_t irow=0;irow<159;irow++) {
169 AliTPCclusterMI *c=track->GetClusterPointer(irow);
171 if (c->GetDetector()!=isec) continue;
172 if (c->GetType()<0) continue;
173 Float_t dx = c->GetX()-kXmean;
174 Float_t y1 = py1[0]+py1[1]*dx;
175 Float_t y2 = py2[0]+py2[1]*dx+py2[2]*dx*dx;
177 Float_t z1 = pz1[0]+pz1[1]*dx;
178 Float_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
188 x[6]=py2[2]*150*150/4; // sagita 150 cm
192 fDiffHistoLine->Fill(x);
195 fDiffHistoPar->Fill(x);
197 if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(py2[2]*150*150/4)<1){
199 AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
200 py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
204 TTreeSRedirector *cstream = GetDebugStreamer();
224 void AliTPCcalibUnlinearity::MakeHisto(){
228 TString axisName[10];
229 Double_t xmin[10], xmax[10];
233 axisName[0]="sector";
234 xmin[0]=0; xmax[0]=72; nbins[0]=72;
237 xmin[1]=0; xmax[1]=159; nbins[1]=159;
240 xmin[2]=-50; xmax[2]=50; nbins[2]=10;
243 xmin[3]=-250; xmax[3]=250; nbins[3]=50;
246 xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
249 xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
252 xmin[6]=-2; xmax[6]=2; nbins[6]=20;
255 xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
258 xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
260 fDiffHistoLine = new THnSparseF("DistHistoLine","DistHistoLine",9,nbins,xmin,xmax);
261 fDiffHistoPar = new THnSparseF("DistHistoPar","DistHistoPar",9,nbins,xmin,xmax);
262 for (Int_t i=0; i<9;i++){
263 fDiffHistoLine->GetAxis(i)->SetName(axisName[i].Data());
264 fDiffHistoPar->GetAxis(i)->SetName(axisName[i].Data());
265 fDiffHistoLine->GetAxis(i)->SetTitle(axisName[i].Data());
266 fDiffHistoPar->GetAxis(i)->SetTitle(axisName[i].Data());
271 void AliTPCcalibUnlinearity::Terminate(){
273 // Terminate function
274 // call base terminate + Eval of fitters
276 Info("AliTPCcalibUnlinearities","Terminate");
279 AliTPCcalibBase::Terminate();
283 void AliTPCcalibUnlinearity::DumpTree(){
287 if (fStreamLevel==0) return;
288 TTreeSRedirector *cstream = GetDebugStreamer();
289 if (!cstream) return;
292 Double_t position[10];
294 Int_t *bins = new Int_t[10];
296 for (Int_t ihis=0; ihis<2; ihis++){
297 his = (ihis==0)? fDiffHistoLine:fDiffHistoPar;
300 Int_t ndim = his->GetNdimensions();
302 for (Long64_t i = 0; i < his->GetNbins(); ++i) {
303 value = his->GetBinContent(i, bins);
304 for (Int_t idim = 0; idim < ndim; idim++) {
305 position[idim] = his->GetAxis(idim)->GetBinCenter(bins[idim]);
307 (*cstream)<<"Resol"<<
310 for (Int_t idim=0;idim<ndim;idim++){
311 (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
313 (*cstream)<<"Resol"<<
320 void AliTPCcalibUnlinearity::AddPoint(Int_t sector, Int_t row, Float_t dz, Float_t dy, Float_t p2, Float_t p3, Float_t dr, Int_t npoints){
323 // Simple distortion fit in outer sectors
325 // Model - 2 exponential decay toward the center of chamber
326 // - Decay length 10 and 20 cm
327 // - Decay amplitude - Z dependent
330 chainUnlin->SetAlias("side","(-1+((sector%36)<18)*2)");
331 chainUnlin->SetAlias("sa","sin(pi*sector*20/180)");
332 chainUnlin->SetAlias("ca","cos(pi*sector*20/180)");
333 chainUnlin->SetAlias("dzt","dz*side");
334 chainUnlin->SetAlias("dr","(1-abs(lz/250))");
335 chainUnlin->SetAlias("dout","(159-row)*1.5");
336 chainUnlin->SetAlias("din","row*0.75");
337 chainUnlin->SetAlias("eout10","exp(-(dout)/10.)");
338 chainUnlin->SetAlias("eout20","exp(-(dout)/20.)");
340 Float_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1
341 Float_t dzt = dz*side;
342 Float_t dout = (159-row)*1.5; // distance to the last pad row - (valid only for long pads)
343 Float_t eout10 = TMath::Exp(-dout/10.);
344 Float_t eout20 = TMath::Exp(-dout/20.);
346 Double_t xxxR[6], xxxZ[6], xxxRZ[6];
347 //TString fstring="";
348 xxxZ[0]=eout20; //fstring+="eout20++";
349 xxxZ[1]=eout10; //fstring+="eout10++";
350 xxxZ[1]=dr*eout20; //fstring+="dr*eout20++";
351 xxxZ[1]=dr*eout10; //fstring+="dr*eout10++";
352 xxxZ[1]=dr*dr*eout20; //fstring+="dr*dr*eout20++";
353 xxxZ[1]=dr*dr*eout10; //fstring+="dr*dr*eout10++";
355 xxxR[0]=p2*eout20; //fstring+="eout20++";
356 xxxR[1]=p2*eout10; //fstring+="eout10++";
357 xxxR[1]=p2*dr*eout20; //fstring+="dr*eout20++";
358 xxxR[1]=p2*dr*eout10; //fstring+="dr*eout10++";
359 xxxR[1]=p2*dr*dr*eout20; //fstring+="dr*dr*eout20++";
360 xxxR[1]=p2*dr*dr*eout10; //fstring+="dr*dr*eout10++";
362 xxxRZ[0]=p3*eout20; //fstring+="eout20++";
363 xxxRZ[1]=p3*eout10; //fstring+="eout10++";
364 xxxRZ[1]=p3*dr*eout20; //fstring+="dr*eout20++";
365 xxxRZ[1]=p3*dr*eout10; //fstring+="dr*eout10++";
366 xxxRZ[1]=p3*dr*dr*eout20; //fstring+="dr*dr*eout20++";
367 xxxRZ[1]=p3*dr*dr*eout10; //fstring+="dr*dr*eout10++";
369 TLinearFitter * fitter=0;
370 for (Int_t ip=0; ip<npoints; ip++){
372 fitter = (TLinearFitter*)fFittersOutR.At(sector%36);
373 fitter->AddPoint(xxxR,dy);
374 fitter->AddPoint(xxxRZ,dz);
375 fitter = (TLinearFitter*)fFittersOutZ.At(sector%36);
376 fitter->AddPoint(xxxZ,dzt);
379 fitter = (TLinearFitter*)fFittersOutR.At(36);
380 fitter->AddPoint(xxxR,dy);
381 fitter->AddPoint(xxxRZ,dz);
382 fitter = (TLinearFitter*)fFittersOutZ.At(36);
383 fitter->AddPoint(xxxZ,dzt);
386 fitter = (TLinearFitter*)fFittersOutR.At(36);
387 fitter->AddPoint(xxxR,dy);
388 fitter->AddPoint(xxxRZ,dz);
389 fitter = (TLinearFitter*)fFittersOutZ.At(36);
390 fitter->AddPoint(xxxZ,dzt);
396 void AliTPCcalibUnlinearity::MakeFitters(){
400 // Make outer fitters
401 TLinearFitter * fitter = 0;
402 for (Int_t ifit=0; ifit<38; ifit++){
403 fitter = new TLinearFitter(7,"hyp6");
404 fitter->StoreData(kFALSE);
405 fFittersOutR.AddAt(fitter,ifit);
406 fitter = new TLinearFitter(7,"hyp6");
407 fitter->StoreData(kFALSE);
408 fFittersOutZ.AddAt(fitter,ifit);
412 void AliTPCcalibUnlinearity::EvalFitters(){
417 TLinearFitter * fitter = 0;
418 for (Int_t ifit=0; ifit<38; ifit++){
419 fitter=(TLinearFitter*)fFittersOutR.At(ifit);
420 if (fitter) fitter->Eval();
421 fitter=(TLinearFitter*)fFittersOutZ.At(ifit);
422 if (fitter) fitter->Eval();
427 void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Int_t nmaxPoints){
431 // TTree * tree = chainUnlinP;
432 AliTPCcalibUnlinearity *calib = this;
434 AliTPCclusterMI * cl=0;
436 TVectorD *vy=0, *vz=0;
437 Int_t nentries = TMath::Min(Int_t(tree->GetEntries()), nmaxPoints);
440 for (Int_t i=0; i<nentries; i++){
441 tree->SetBranchAddress("Cl.",&cl);
442 tree->SetBranchAddress("y1",&ty);
443 tree->SetBranchAddress("z1",&tz);
444 tree->SetBranchAddress("py1.",&vy);
445 tree->SetBranchAddress("pz1.",&vz);
449 calib->AddPoint(cl->GetDetector(), cl->GetRow(), cl->GetZ()-tz, cl->GetY()-ty,
450 (*vy)(1),(*vz)(1), 1-TMath::Abs(cl->GetZ()/250),1);