]>
Commit | Line | Data |
---|---|---|
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 | ||
37fb53c1 | 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 | ||
c11fe047 | 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" | |
28b5b7c8 | 64 | #include "TCut.h" |
65 | ||
c11fe047 | 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" | |
f7a1cc68 | 75 | #include "AliMagF.h" |
c11fe047 | 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" | |
c11fe047 | 85 | #include "AliDCSSensorArray.h" |
86 | #include "AliDCSSensor.h" | |
87 | #include "TLinearFitter.h" | |
28b5b7c8 | 88 | #include "AliTPCClusterParam.h" |
89 | #include "TStatToolkit.h" | |
c11fe047 | 90 | |
91 | ||
92 | #include "AliTPCcalibUnlinearity.h" | |
37fb53c1 | 93 | #include "AliTPCPointCorrection.h" |
c11fe047 | 94 | |
95 | ||
96 | ClassImp(AliTPCcalibUnlinearity) | |
97 | ||
98 | AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(): | |
99 | AliTPCcalibBase(), | |
37fb53c1 | 100 | fInit(kFALSE), |
101 | fDiffHistoLineY(0), | |
102 | fDiffHistoParY(0), | |
103 | fDiffHistoLineZ(0), | |
104 | fDiffHistoParZ(0), | |
c11fe047 | 105 | fFittersOutR(38), |
28b5b7c8 | 106 | fFittersOutZ(38), |
107 | fParamsOutR(38), | |
108 | fParamsOutZ(38), | |
109 | fErrorsOutR(38), | |
37fb53c1 | 110 | fErrorsOutZ(38), |
111 | fDistRPHIPlus(74), | |
112 | fDistRPHIMinus(74), | |
113 | fFitterQuadrantY(16*38), | |
114 | fFitterQuadrantPhi(16*38) | |
c11fe047 | 115 | { |
116 | // | |
117 | // Defualt constructor | |
118 | // | |
119 | } | |
120 | ||
121 | AliTPCcalibUnlinearity::AliTPCcalibUnlinearity(const Text_t *name, const Text_t *title): | |
122 | AliTPCcalibBase(name,title), | |
37fb53c1 | 123 | fInit(kFALSE), |
124 | fDiffHistoLineY(0), | |
125 | fDiffHistoParY(0), | |
126 | fDiffHistoLineZ(0), | |
127 | fDiffHistoParZ(0), | |
c11fe047 | 128 | fFittersOutR(38), |
28b5b7c8 | 129 | fFittersOutZ(38), |
130 | fParamsOutR(38), | |
131 | fParamsOutZ(38), | |
132 | fErrorsOutR(38), | |
37fb53c1 | 133 | fErrorsOutZ(38), |
134 | fDistRPHIPlus(74), | |
135 | fDistRPHIMinus(74), | |
136 | fFitterQuadrantY(16*38), | |
137 | fFitterQuadrantPhi(16*38) | |
c11fe047 | 138 | { |
139 | // | |
140 | // Non default constructor | |
141 | // | |
37fb53c1 | 142 | MakeHisto(); |
c11fe047 | 143 | } |
144 | ||
145 | AliTPCcalibUnlinearity::~AliTPCcalibUnlinearity(){ | |
146 | // | |
147 | // | |
148 | // | |
37fb53c1 | 149 | if (fDiffHistoLineZ) delete fDiffHistoLineY; |
150 | if (fDiffHistoParZ) delete fDiffHistoParY; | |
151 | if (fDiffHistoLineZ) delete fDiffHistoLineZ; | |
152 | if (fDiffHistoParZ) delete fDiffHistoParZ; | |
c11fe047 | 153 | } |
154 | ||
155 | ||
28b5b7c8 | 156 | |
157 | ||
158 | ||
c11fe047 | 159 | void AliTPCcalibUnlinearity::Process(AliTPCseed *seed) { |
160 | // | |
161 | // | |
162 | // | |
37fb53c1 | 163 | if (!fInit) { |
164 | Init(); // work around for PROOF - initialize firs time used | |
165 | } | |
c11fe047 | 166 | const Int_t kdrow=10; |
167 | const Int_t kMinCluster=35; | |
c11fe047 | 168 | if (!seed) return; |
28b5b7c8 | 169 | if (TMath::Abs(fMagF)>0.01) return; // use only data without field |
c11fe047 | 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); | |
37fb53c1 | 181 | if (counter[is]>kMinCluster) { |
182 | AlignOROC(seed, is); | |
183 | } | |
c11fe047 | 184 | } |
185 | } | |
186 | ||
187 | ||
188 | void 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; | |
37fb53c1 | 194 | const Float_t kSagitaCut = 1; |
c11fe047 | 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; | |
28b5b7c8 | 215 | Double_t dx = c->GetX()-kXmean; |
c11fe047 | 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; | |
28b5b7c8 | 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; | |
c11fe047 | 235 | // |
28b5b7c8 | 236 | Double_t z1 = pz1[0]+pz1[1]*dx; |
237 | Double_t z2 = pz2[0]+pz2[1]*dx+pz2[2]*dx*dx; | |
c11fe047 | 238 | // |
37fb53c1 | 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 | // | |
c11fe047 | 255 | // |
c11fe047 | 256 | // |
c11fe047 | 257 | if (fStreamLevel>1){ |
258 | TTreeSRedirector *cstream = GetDebugStreamer(); | |
259 | if (cstream){ | |
260 | (*cstream)<<"Diff"<< | |
28b5b7c8 | 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 | |
37fb53c1 | 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 | // | |
c11fe047 | 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<< | |
37fb53c1 | 283 | "eP="<<edgePlus<< |
284 | "eM="<<edgeMinus<< | |
c11fe047 | 285 | "\n"; |
286 | } | |
287 | } | |
37fb53c1 | 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 | ||
c11fe047 | 343 | } |
344 | } | |
345 | ||
346 | ||
37fb53c1 | 347 | void 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.="<<¶ms[i0]<< // parameters | |
529 | "p1.="<<¶ms[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 | ||
c11fe047 | 546 | void 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; | |
37fb53c1 | 557 | // |
c11fe047 | 558 | axisName[1]="row"; |
37fb53c1 | 559 | xmin[1]=0; xmax[1]=96; nbins[1]=96; |
c11fe047 | 560 | // |
37fb53c1 | 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 | ||
c11fe047 | 590 | // |
37fb53c1 | 591 | // |
c11fe047 | 592 | // |
37fb53c1 | 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); | |
c11fe047 | 604 | } |
605 | } | |
606 | ||
607 | ||
608 | void 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 | ||
37fb53c1 | 619 | Long64_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 | ||
635 | void 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 | ||
c11fe047 | 684 | |
37fb53c1 | 685 | void AliTPCcalibUnlinearity::DumpTree(const char *fname){ |
c11fe047 | 686 | // |
687 | // | |
688 | // | |
37fb53c1 | 689 | TTreeSRedirector *cstream =new TTreeSRedirector(fname); |
c11fe047 | 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++){ | |
37fb53c1 | 698 | his = (ihis==0)? fDiffHistoLineY:fDiffHistoParY; |
c11fe047 | 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 | } | |
37fb53c1 | 718 | delete cstream; |
c11fe047 | 719 | } |
720 | ||
721 | ||
37fb53c1 | 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){ |
c11fe047 | 723 | // |
724 | // | |
725 | // Simple distortion fit in outer sectors | |
726 | // | |
727 | // Model - 2 exponential decay toward the center of chamber | |
28b5b7c8 | 728 | // - Decay length 10 and 5 cm |
c11fe047 | 729 | // - Decay amplitude - Z dependent |
730 | // | |
37fb53c1 | 731 | |
28b5b7c8 | 732 | Double_t side = (-1+((sector%36)<18)*2); // A side +1 C side -1 |
37fb53c1 | 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; | |
28b5b7c8 | 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 | ||
c11fe047 | 748 | // |
749 | Double_t xxxR[6], xxxZ[6], xxxRZ[6]; | |
750 | //TString fstring=""; | |
28b5b7c8 | 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 | // | |
37fb53c1 | 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++"; | |
28b5b7c8 | 764 | // |
37fb53c1 | 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++"; | |
c11fe047 | 771 | // |
772 | TLinearFitter * fitter=0; | |
28b5b7c8 | 773 | Double_t err=0.1; |
c11fe047 | 774 | for (Int_t ip=0; ip<npoints; ip++){ |
775 | // | |
776 | fitter = (TLinearFitter*)fFittersOutR.At(sector%36); | |
28b5b7c8 | 777 | fitter->AddPoint(xxxR,dy,err); |
778 | //fitter->AddPoint(xxxRZ,dz); | |
c11fe047 | 779 | fitter = (TLinearFitter*)fFittersOutZ.At(sector%36); |
28b5b7c8 | 780 | fitter->AddPoint(xxxZ,dzt,err); |
c11fe047 | 781 | // |
782 | if (side>0){ | |
783 | fitter = (TLinearFitter*)fFittersOutR.At(36); | |
28b5b7c8 | 784 | fitter->AddPoint(xxxR,dy,err); |
785 | //fitter->AddPoint(xxxRZ,dz); | |
c11fe047 | 786 | fitter = (TLinearFitter*)fFittersOutZ.At(36); |
28b5b7c8 | 787 | fitter->AddPoint(xxxZ,dzt,err); |
c11fe047 | 788 | } |
789 | if (side<0){ | |
28b5b7c8 | 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); | |
c11fe047 | 795 | } |
796 | } | |
797 | } | |
37fb53c1 | 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){ |
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 | ||
c11fe047 | 820 | |
821 | ||
37fb53c1 | 822 | void AliTPCcalibUnlinearity::Init(){ |
c11fe047 | 823 | // |
824 | // | |
825 | // | |
37fb53c1 | 826 | // |
827 | MakeHisto(); | |
c11fe047 | 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 | } | |
37fb53c1 | 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; | |
c11fe047 | 847 | } |
848 | ||
849 | void AliTPCcalibUnlinearity::EvalFitters(){ | |
850 | // | |
851 | // | |
852 | // Evaluate fitters | |
853 | // | |
854 | TLinearFitter * fitter = 0; | |
28b5b7c8 | 855 | TVectorD vec; |
c11fe047 | 856 | for (Int_t ifit=0; ifit<38; ifit++){ |
857 | fitter=(TLinearFitter*)fFittersOutR.At(ifit); | |
28b5b7c8 | 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 | } | |
c11fe047 | 865 | fitter=(TLinearFitter*)fFittersOutZ.At(ifit); |
28b5b7c8 | 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 | } | |
c11fe047 | 873 | } |
874 | } | |
875 | ||
28b5b7c8 | 876 | |
877 | ||
878 | ||
879 | ||
880 | void AliTPCcalibUnlinearity::ProcessTree(TTree * tree, Long64_t nmaxPoints){ | |
c11fe047 | 881 | // |
882 | // | |
883 | // | |
884 | // TTree * tree = chainUnlinP; | |
885 | AliTPCcalibUnlinearity *calib = this; | |
886 | // | |
887 | AliTPCclusterMI * cl=0; | |
28b5b7c8 | 888 | Double_t ty,tz; |
c11fe047 | 889 | TVectorD *vy=0, *vz=0; |
28b5b7c8 | 890 | TVectorD *vy2=0, *vz2=0; |
891 | Long64_t nentries = tree->GetEntries(); | |
892 | if (nmaxPoints>0) nentries = TMath::Min(nentries,nmaxPoints); | |
c11fe047 | 893 | // |
894 | { | |
28b5b7c8 | 895 | for (Long64_t i=0; i<nentries; i++){ |
c11fe047 | 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); | |
28b5b7c8 | 901 | tree->SetBranchAddress("py2.",&vy2); |
902 | tree->SetBranchAddress("pz2.",&vz2); | |
c11fe047 | 903 | // |
904 | tree->GetEntry(i); | |
905 | if (!cl) continue; | |
28b5b7c8 | 906 | if (i%10000==0) printf("%d\n",(Int_t)i); |
907 | Int_t row= cl->GetRow(); | |
908 | if (cl->GetDetector()>36) row+=64; | |
37fb53c1 | 909 | calib->AddPointRPHI(cl->GetDetector(),cl->GetX(),cl->GetY(),cl->GetZ(), |
910 | ty,tz,(*vy)[1],(*vz)[1],1); | |
911 | ||
28b5b7c8 | 912 | if (cl->GetType()<0) continue; |
913 | Double_t dy = cl->GetY()-ty; | |
914 | Double_t dz = cl->GetZ()-tz; | |
28b5b7c8 | 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 | |
37fb53c1 | 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 | } | |
c11fe047 | 927 | } |
928 | } | |
929 | } | |
28b5b7c8 | 930 | |
931 | ||
932 | void 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 | ||
162637e4 | 959 | AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); |
28b5b7c8 | 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(); | |
798017c7 | 1035 | param->PosZcor(0) = (TVectorD*) fitParamZ0.Clone(); |
28b5b7c8 | 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(); | |
798017c7 | 1041 | param->PosZcor(1) = (TVectorD*) fitParamZ1.Clone(); |
1042 | param->PosZcor(2) = (TVectorD*) fitParamZ1.Clone(); | |
28b5b7c8 | 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(); | |
798017c7 | 1049 | param->PosYcor(0) = (TVectorD*) fitParamY0.Clone(); |
28b5b7c8 | 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(); | |
798017c7 | 1056 | param->PosYcor(1) = (TVectorD*) fitParamY1.Clone(); |
1057 | param->PosYcor(2) = (TVectorD*) fitParamY1.Clone(); | |
28b5b7c8 | 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 |