]> 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
New classes for finding multiple vertices (in case of pile-up). They will be used...
[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
18 Only Linear and parabolic fit used
19 To be used for laser or for data without field
20 //
21*/
22
23
24#include "TLinearFitter.h"
25
26#include "Riostream.h"
27#include "TChain.h"
28#include "TTree.h"
29#include "TH1F.h"
30#include "TH2F.h"
31#include "TH3F.h"
32#include "THnSparse.h"
33#include "TList.h"
34#include "TMath.h"
35#include "TCanvas.h"
36#include "TFile.h"
37#include "TF1.h"
38#include "TVectorD.h"
39#include "TProfile.h"
40#include "TGraphErrors.h"
41#include "TCanvas.h"
42
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"
53
54#include "AliLog.h"
55
56
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"
65
66
67#include "AliTPCcalibUnlinearity.h"
68
69
70ClassImp(AliTPCcalibUnlinearity)
71
72AliTPCcalibUnlinearity::AliTPCcalibUnlinearity():
73 AliTPCcalibBase(),
74 fDiffHistoLine(0),
75 fDiffHistoPar(0),
76 fFittersOutR(38),
77 fFittersOutZ(38)
78{
79 //
80 // Defualt constructor
81 //
82}
83
84AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title):
85 AliTPCcalibBase(name,title),
86 fDiffHistoLine(0),
87 fDiffHistoPar(0),
88 fFittersOutR(38),
89 fFittersOutZ(38)
90{
91 //
92 // Non default constructor
93 //
94 MakeFitters();
95}
96
97AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){
98 //
99 //
100 //
101 if (fDiffHistoLine) delete fDiffHistoLine;
102 if (fDiffHistoPar) delete fDiffHistoPar;
103}
104
105
106void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) {
107 //
108 //
109 //
110 const Int_t kdrow=10;
111 const Int_t kMinCluster=35;
112 if (!fDiffHistoLine) MakeHisto();
113 if (!seed) return;
114 Int_t counter[72];
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()]++;
121 }
122 //
123 for (Int_t is=0; is<72;is++){
124 if (counter[is]>kMinCluster) ProcessDiff(seed, is);
125 }
126}
127
128
129void AliTPCcalibUnlinearity::ProcessDiff(AliTPCseed *track, Int_t isec){
130 //
131 // process differnce of the cluster in respect with linear and parabolic fit
132 //
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");
139 //
140 static TVectorD py1(2);
141 static TVectorD py2(3);
142 static TVectorD pz1(2);
143 static TVectorD pz2(3);
144 //
145 fy1.ClearPoints();
146 fy2.ClearPoints();
147 fz1.ClearPoints();
148 fz2.ClearPoints();
149 //
150 for (Int_t irow=kdrow;irow<159-kdrow;irow++) {
151 AliTPCclusterMI *c=track->GetClusterPointer(irow);
152 if (!c) continue;
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());
161 }
162 fy1.Eval(); fy1.GetParameters(py1);
163 fy2.Eval(); fy2.GetParameters(py2);
164 fz1.Eval(); fz1.GetParameters(pz1);
165 fz2.Eval(); fz2.GetParameters(pz2);
166 //
167
168 for (Int_t irow=0;irow<159;irow++) {
169 AliTPCclusterMI *c=track->GetClusterPointer(irow);
170 if (!c) continue;
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;
176 //
177 Float_t z1 = pz1[0]+pz1[1]*dx;
178 Float_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx;
179 //
180 //
181 Double_t x[10];
182 x[0]=isec;
183 x[1]=irow;
184 x[2]=c->GetY();
185 x[3]=c->GetZ();
186 x[4]=py1[1];
187 x[5]=pz1[1];
188 x[6]=py2[2]*150*150/4; // sagita 150 cm
189 //
190 x[7]=c->GetY()-y1;
191 x[8]=c->GetZ()-z1;
192 fDiffHistoLine->Fill(x);
193 x[7]=c->GetY()-y2;
194 x[8]=c->GetZ()-z2;
195 fDiffHistoPar->Fill(x);
196
197 if (TMath::Abs(py2[2]*150*150/4)<1 && TMath::Abs(py2[2]*150*150/4)<1){
198 // Apply sagita cut
199 AddPoint(isec,irow,c->GetZ()-z1, c->GetY()-y1,
200 py1[1],pz1[1],1-TMath::Abs(c->GetZ()/250.),1);
201 }
202
203 if (fStreamLevel>1){
204 TTreeSRedirector *cstream = GetDebugStreamer();
205 if (cstream){
206 (*cstream)<<"Diff"<<
207 "isec="<<isec<<
208 "Cl.="<<c<<
209 "y1="<<y1<<
210 "y2="<<y2<<
211 "z1="<<z1<<
212 "z2="<<z2<<
213 "py1.="<<&py1<<
214 "py2.="<<&py2<<
215 "pz1.="<<&pz1<<
216 "pz2.="<<&pz2<<
217 "\n";
218 }
219 }
220 }
221}
222
223
224void AliTPCcalibUnlinearity::MakeHisto(){
225 //
226 //
227 //
228 TString axisName[10];
229 Double_t xmin[10], xmax[10];
230 Int_t nbins[10];
231 //
232 //
233 axisName[0]="sector";
234 xmin[0]=0; xmax[0]=72; nbins[0]=72;
235 //
236 axisName[1]="row";
237 xmin[1]=0; xmax[1]=159; nbins[1]=159;
238 //
239 axisName[2]="ly";
240 xmin[2]=-50; xmax[2]=50; nbins[2]=10;
241 //
242 axisName[3]="lz";
243 xmin[3]=-250; xmax[3]=250; nbins[3]=50;
244 //
245 axisName[4]="p2";
246 xmin[4]=-0.6; xmax[4]=0.6; nbins[4]=12;
247 //
248 axisName[5]="p3";
249 xmin[5]=-0.6; xmax[5]=0.6; nbins[5]=12;
250 //
251 axisName[6]="p4";
252 xmin[6]=-2; xmax[6]=2; nbins[6]=20;
253 //
254 axisName[7]="dy";
255 xmin[7]=-0.5; xmax[7]=0.5; nbins[7]=100;
256 //
257 axisName[8]="dz";
258 xmin[8]=-0.5; xmax[8]=0.5; nbins[8]=100;
259 //
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());
267 }
268}
269
270
271void AliTPCcalibUnlinearity::Terminate(){
272 //
273 // Terminate function
274 // call base terminate + Eval of fitters
275 //
276 Info("AliTPCcalibUnlinearities","Terminate");
277 EvalFitters();
278 DumpTree();
279 AliTPCcalibBase::Terminate();
280}
281
282
283void AliTPCcalibUnlinearity::DumpTree(){
284 //
285 //
286 //
287 if (fStreamLevel==0) return;
288 TTreeSRedirector *cstream = GetDebugStreamer();
289 if (!cstream) return;
290 //
291 THnSparse *his=0;
292 Double_t position[10];
293 Double_t value;
294 Int_t *bins = new Int_t[10];
295 //
296 for (Int_t ihis=0; ihis<2; ihis++){
297 his = (ihis==0)? fDiffHistoLine:fDiffHistoPar;
298 if (!his) continue;
299 //
300 Int_t ndim = his->GetNdimensions();
301 //
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]);
306 }
307 (*cstream)<<"Resol"<<
308 "ihis="<<ihis<<
309 "bincont="<<value;
310 for (Int_t idim=0;idim<ndim;idim++){
311 (*cstream)<<"Resol"<<Form("%s=", his->GetAxis(idim)->GetName())<<position[idim];
312 }
313 (*cstream)<<"Resol"<<
314 "\n";
315 }
316 }
317}
318
319
320void 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){
321 //
322 //
323 // Simple distortion fit in outer sectors
324 //
325 // Model - 2 exponential decay toward the center of chamber
326 // - Decay length 10 and 20 cm
327 // - Decay amplitude - Z dependent
328 //
329 /*
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.)");
339 */
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.);
345 //
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++";
354 //
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++";
361 //
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++";
368 //
369 TLinearFitter * fitter=0;
370 for (Int_t ip=0; ip<npoints; ip++){
371 //
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);
377 //
378 if (side>0){
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);
384 }
385 if (side<0){
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);
391 }
392 }
393}
394
395
396void AliTPCcalibUnlinearity::MakeFitters(){
397 //
398 //
399 //
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);
409 }
410}
411
412void AliTPCcalibUnlinearity::EvalFitters(){
413 //
414 //
415 // Evaluate fitters
416 //
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();
423 }
424}
425
426
427void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Int_t nmaxPoints){
428 //
429 //
430 //
431 // TTree * tree = chainUnlinP;
432 AliTPCcalibUnlinearity *calib = this;
433 //
434 AliTPCclusterMI * cl=0;
435 Float_t ty,tz;
436 TVectorD *vy=0, *vz=0;
437 Int_t nentries = TMath::Min(Int_t(tree->GetEntries()), nmaxPoints);
438 //
439 {
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);
446 //
447 tree->GetEntry(i);
448 if (!cl) continue;
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);
451 }
452 }
453}