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 **************************************************************************/
20 // 1. The laser track is associated with the mirror
21 // see function FindMirror
23 // 2. The laser track is accepted for the analysis under certain condition
24 // (see function Accpet laser)
26 // 3. The drift velocity and jitter is calculated event by event
27 // (see function drift velocity)
29 // 4. The tracks are refitted at different sectors
33 // 4.c) parabola with common P2 for inner and outer
35 // To make laser scan the user interaction neccessary
38 gSystem->Load("libANALYSIS");
39 gSystem->Load("libTPCcalib");
40 TFile fcalib("CalibObjectsTrain2.root");
41 AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
42 laser->DumpMeanInfo(run)
43 TFile fmean("laserMean.root")
45 // laser track clasification;
47 TCut cutT("cutT","abs(Tr.fP[3])<0.06");
48 TCut cutN("cutN","fTPCncls>70");
49 TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
50 TCut cutA = cutT+cutPt+cutP;
51 TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
58 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
59 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
60 AliXRDPROOFtoolkit tool;
61 TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
63 AliTPCcalibLaser::DumpScanInfo(chain)
64 TFile fscan("laserScan.root")
65 TTree * treeT = (TTree*)fscan.Get("Mean")
69 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
70 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
71 AliXRDPROOFtoolkit tool;
72 AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1)
73 TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
75 TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
76 chainDriftN->Lookup();
79 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
81 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
83 TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
90 #include "TLinearFitter.h"
91 #include "AliTPCcalibLaser.h"
92 #include "AliExternalTrackParam.h"
93 #include "AliESDEvent.h"
94 #include "AliESDfriend.h"
95 #include "AliESDtrack.h"
96 #include "AliTPCTracklet.h"
100 #include "TVectorD.h"
101 #include "TTreeStream.h"
104 #include "TGraphErrors.h"
105 #include "AliTPCclusterMI.h"
106 #include "AliTPCseed.h"
107 #include "AliTracker.h"
109 #include "TClonesArray.h"
115 #include "TStatToolkit.h"
117 #include "TDatabasePDG.h"
120 #include "TTreeStream.h"
123 #include "AliTPCLaserTrack.h"
124 #include "AliTPCcalibDB.h"
125 #include "AliTPCParam.h"
126 #include "AliTPCParamSR.h"
127 #include "TTimeStamp.h"
128 #include "AliDCSSensorArray.h"
129 #include "AliDCSSensor.h"
130 #include "AliGRPObject.h"
131 #include "AliTPCROC.h"
135 ClassImp(AliTPCcalibLaser)
137 AliTPCcalibLaser::AliTPCcalibLaser():
144 fTracksEsdParam(336),
154 fHisLaser(0), // N dim histogram of laser
155 fHisLaserPad(0), // N dim histogram of laser
156 fHisLaserTime(0), // N dim histogram of laser
157 fHisNclIn(0), //->Number of clusters inner
158 fHisNclOut(0), //->Number of clusters outer
159 fHisNclIO(0), //->Number of cluster inner outer
160 fHisLclIn(0), //->Level arm inner
161 fHisLclOut(0), //->Level arm outer
162 fHisLclIO(0), //->Number of cluster inner outer
163 fHisdEdx(0), //->dEdx histo
164 fHisdZfit(0), //->distance to the mirror after linear fit
167 fHisChi2YIn1(0), //->chi2 y inner - line
168 fHisChi2YOut1(0), //->chi2 y inner - line
169 fHisChi2YIn2(0), //->chi2 y inner - parabola
170 fHisChi2YOut2(0), //->chi2 y inner - parabola
171 fHisChi2YIO1(0), //->chi2 y IO - common
172 fHisChi2ZIn1(0), //->chi2 z inner - line
173 fHisChi2ZOut1(0), //->chi2 z inner - line
174 fHisChi2ZIn2(0), //->chi2 z inner - parabola
175 fHisChi2ZOut2(0), //->chi2 z inner - parabola
176 fHisChi2ZIO1(0), //->chi2 z IO - common
179 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
180 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
181 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
182 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
183 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
184 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
185 fHisPy2vP2In(0), //-> Curv P2inner - parabola
186 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
187 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
190 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
191 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
192 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
193 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
194 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
195 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
196 fHisPz2vP2In(0), //-> Curv P2inner - parabola
197 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
198 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
200 fDeltaYres(336), //->2D histo of residuals
201 fDeltaZres(336), //->2D histo fo residuals
202 fDeltaYres2(336), //->2D histo of residuals
203 fDeltaZres2(336), //->2D histo fo residuals
204 fDeltaYresAbs(336), //->2D histo of residuals
205 fHisYAbsErrors(0), //-> total errors per beam in the abs res y analysis
206 fDeltaZresAbs(336), //->2D histo of residuals
207 fHisZAbsErrors(0), //-> total errors per beam in the abs res z analysis
208 //fDeltaYres3(336), //->2D histo of residuals
209 //fDeltaZres3(336), //->2D histo fo residuals
210 fFitAside(new TVectorD(5)),
211 fFitCside(new TVectorD(5)),
212 fFitACside(new TVectorD(6)),
217 fBeamSectorOuter(336),
218 fBeamSectorInner(336),
219 fBeamOffsetYOuter(336),
220 fBeamSlopeYOuter(336),
221 fBeamOffsetYInner(336),
222 fBeamSlopeYInner(336),
223 fBeamOffsetZOuter(336),
224 fBeamSlopeZOuter(336),
225 fBeamOffsetZInner(336),
226 fBeamSlopeZInner(336),
227 fInverseSlopeZ(kTRUE),
229 fFixedFitAside0(0.0),
230 fFixedFitAside1(1.0),
231 fFixedFitCside0(0.0),
237 fTracksEsdParam.SetOwner(kTRUE);
238 for (Int_t i=0; i<336; i++) {
240 fCounter[i]=0; //! counter of usage
241 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
242 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
246 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
253 fTracksEsdParam(336),
257 fDeltaZ(336), // array of histograms of delta z for each track
258 fDeltaP3(336), // array of histograms of delta z for each track
259 fDeltaP4(336), // array of histograms of P3 for each track
260 fDeltaPhi(336), // array of histograms of P4 for each track
261 fDeltaPhiP(336), // array of histograms of delta z for each track
262 fSignals(336), // array of dedx signals
265 fHisLaser(0), // N dim histogram of laser
266 fHisLaserPad(0), // N dim histogram of laser
267 fHisLaserTime(0), // N dim histogram of laser
269 fHisNclIn(0), //->Number of clusters inner
270 fHisNclOut(0), //->Number of clusters outer
271 fHisNclIO(0), //->Number of cluster inner outer
272 fHisLclIn(0), //->Level arm inner
273 fHisLclOut(0), //->Level arm outer
274 fHisLclIO(0), //->Number of cluster inner outer
275 fHisdEdx(0), //->dEdx histo
276 fHisdZfit(0), //->distance to the mirror after linear fit
279 fHisChi2YIn1(0), //->chi2 y inner - line
280 fHisChi2YOut1(0), //->chi2 y inner - line
281 fHisChi2YIn2(0), //->chi2 y inner - parabola
282 fHisChi2YOut2(0), //->chi2 y inner - parabola
283 fHisChi2YIO1(0), //->chi2 y IO - common
284 fHisChi2ZIn1(0), //->chi2 z inner - line
285 fHisChi2ZOut1(0), //->chi2 z inner - line
286 fHisChi2ZIn2(0), //->chi2 z inner - parabola
287 fHisChi2ZOut2(0), //->chi2 z inner - parabola
288 fHisChi2ZIO1(0), //->chi2 z IO - common
291 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
292 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
293 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
294 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
295 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
296 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
297 fHisPy2vP2In(0), //-> Curv P2inner - parabola
298 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
299 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
302 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
303 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
304 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
305 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
306 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
307 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
308 fHisPz2vP2In(0), //-> Curv P2inner - parabola
309 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
310 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
324 fFitAside(new TVectorD(5)), // drift fit - A side
325 fFitCside(new TVectorD(5)), // drift fit - C- side
326 fFitACside(new TVectorD(6)), // drift fit - AC- side
327 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
328 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
329 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
330 fNcuts(0), // number of cuts
331 fBeamSectorOuter(336),
332 fBeamSectorInner(336),
333 fBeamOffsetYOuter(336),
334 fBeamSlopeYOuter(336),
335 fBeamOffsetYInner(336),
336 fBeamSlopeYInner(336),
337 fBeamOffsetZOuter(336),
338 fBeamSlopeZOuter(336),
339 fBeamOffsetZInner(336),
340 fBeamSlopeZInner(336),
341 fInverseSlopeZ(kTRUE),
343 fFixedFitAside0(0.0),
344 fFixedFitAside1(1.0),
345 fFixedFitCside0(0.0),
353 fTracksEsdParam.SetOwner(kTRUE);
354 for (Int_t i=0; i<336; i++) {
356 fCounter[i]=0; //! counter of usage
357 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
358 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
362 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
363 AliTPCcalibBase(calibLaser),
369 fTracksEsdParam(336),
371 fFullCalib(calibLaser.fFullCalib),
373 fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
374 fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
375 fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
376 fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
377 fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
378 fSignals(((calibLaser.fSignals))), // array of dedx signals
381 fHisLaser(0), // N dim histogram of laser
382 fHisLaserPad(0), // N dim histogram of laser
383 fHisLaserTime(0), // N dim histogram of laser
385 fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
386 fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
387 fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
388 fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
389 fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
390 fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
391 fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
392 fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
395 fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
396 fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
397 fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
398 fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
399 fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
400 fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
401 fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
402 fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
403 fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
404 fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
407 fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
408 fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
409 fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
410 fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
411 fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
412 fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
413 fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
414 fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
415 fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
418 fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
419 fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
420 fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
421 fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
422 fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
423 fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
424 fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
425 fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
426 fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
429 fDeltaYres(((calibLaser.fDeltaYres))),
430 fDeltaZres(((calibLaser.fDeltaZres))),
431 fDeltaYres2(((calibLaser.fDeltaYres))),
432 fDeltaZres2(((calibLaser.fDeltaZres))),
433 fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
434 fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
435 fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
436 fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
437 // fDeltaYres3(((calibLaser.fDeltaYres))),
438 //fDeltaZres3(((calibLaser.fDeltaZres))),
439 fFitAside(new TVectorD(5)), // drift fit - A side
440 fFitCside(new TVectorD(5)), // drift fit - C- side
441 fFitACside(new TVectorD(6)), // drift fit - C- side
442 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
443 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
444 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
445 fNcuts(0), // number of cuts
446 fBeamSectorOuter(336),
447 fBeamSectorInner(336),
448 fBeamOffsetYOuter(336),
449 fBeamSlopeYOuter(336),
450 fBeamOffsetYInner(336),
451 fBeamSlopeYInner(336),
452 fBeamOffsetZOuter(336),
453 fBeamSlopeZOuter(336),
454 fBeamOffsetZInner(336),
455 fBeamSlopeZInner(336),
456 fInverseSlopeZ(calibLaser.fInverseSlopeZ),
457 fUseFixedDriftV(calibLaser.fUseFixedDriftV),
458 fFixedFitAside0(calibLaser.fFixedFitAside0),
459 fFixedFitAside1(calibLaser.fFixedFitAside1),
460 fFixedFitCside0(calibLaser.fFixedFitCside0),
461 fFixedFitCside1(calibLaser.fFixedFitCside1)
466 for (Int_t i=0; i<336; i++) {
468 fCounter[i]=0; //! counter of usage
469 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
470 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
476 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
478 // assgnment operator
480 if (this != &calibLaser) {
481 new (this) AliTPCcalibLaser(calibLaser);
490 AliTPCcalibLaser::~AliTPCcalibLaser() {
495 delete fHisLaser; //->
496 delete fHisLaserPad; //->
497 delete fHisLaserTime; //->
499 delete fHisNclIn; //->Number of clusters inner
500 delete fHisNclOut; //->Number of clusters outer
501 delete fHisNclIO; //->Number of cluster inner outer
502 delete fHisLclIn; //->Level arm inner
503 delete fHisLclOut; //->Level arm outer
504 delete fHisLclIO; //->Number of cluster inner outer
509 delete fHisChi2YIn1; //->chi2 y inner - line
510 delete fHisChi2YOut1; //->chi2 y inner - line
511 delete fHisChi2YIn2; //->chi2 y inner - parabola
512 delete fHisChi2YOut2; //->chi2 y inner - parabola
513 delete fHisChi2YIO1; //->chi2 y IO - common
514 delete fHisChi2ZIn1; //->chi2 z inner - line
515 delete fHisChi2ZOut1; //->chi2 z inner - line
516 delete fHisChi2ZIn2; //->chi2 z inner - parabola
517 delete fHisChi2ZOut2; //->chi2 z inner - parabola
518 delete fHisChi2ZIO1; //->chi2 z IO - common
521 delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
522 delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
523 delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
524 delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
525 delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
526 delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
527 delete fHisPy2vP2In; //-> Curv P2inner - parabola
528 delete fHisPy2vP2Out; //-> Curv P2outer - parabola
529 delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
532 delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
533 delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
534 delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
535 delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
536 delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
537 delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
538 delete fHisPz2vP2In; //-> Curv P2inner - parabola
539 delete fHisPz2vP2Out; //-> Curv P2outer - parabola
540 delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
545 fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
546 fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
547 fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
548 fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
549 fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
550 fSignals.SetOwner(); //->Array of dedx signals
552 fDeltaZ.Delete(); //-> array of histograms of delta z for each track
553 fDeltaP3.Delete(); //-> array of histograms of P3 for each track
554 fDeltaP4.Delete(); //-> array of histograms of P4 for each track
555 fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
556 fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
557 fSignals.Delete(); //->Array of dedx signals
559 fDeltaYres.SetOwner();
561 delete fHisYAbsErrors;
562 fDeltaZres.SetOwner();
564 delete fHisZAbsErrors;
565 fDeltaYres2.SetOwner();
566 fDeltaYres2.Delete();
567 fDeltaZres2.SetOwner();
568 fDeltaZres2.Delete();
570 fDeltaYresAbs.SetOwner();
571 fDeltaYresAbs.Delete();
572 fDeltaZresAbs.SetOwner();
573 fDeltaZresAbs.Delete();
578 void AliTPCcalibLaser::Process(AliESDEvent * event) {
581 // Loop over tracks and call Process function
583 const Int_t kMinTracks=20;
584 const Int_t kMinClusters=40;
590 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
594 if (fESDfriend->TestSkipBit()) return;
595 if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
596 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
598 // find CE background if present
600 if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
601 TH1D hisCE("hhisCE","hhisCE",100,-100,100);
602 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
603 AliESDtrack *track=fESD->GetTrack(i);
604 if (!track) continue;
605 hisCE.Fill(track->GetZ());
606 hisCE.Fill(track->GetZ()+2);
607 hisCE.Fill(track->GetZ()-2);
615 fTracksEsdParam.Delete();
616 for (Int_t id=0; id<336;id++) {
618 fClusterCounter[id]=0;
622 Int_t n=fESD->GetNumberOfTracks();
624 for (Int_t i=0;i<n;++i) {
625 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
626 if (!friendTrack) continue;
627 AliESDtrack *track=fESD->GetTrack(i);
628 if (!track) continue;
629 Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
630 if (binC>336) continue; //remove CE background
631 TObject *calibObject=0;
633 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
634 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
636 if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
638 Int_t id = FindMirror(track,seed);
639 if (id>=0) counter++;
644 if (counter<kMinTracks) return;
648 if (!fFullCalib) return;
649 static Bool_t init=kFALSE;
651 init = kTRUE; // way around for PROOF - to be investigated
655 for (Int_t id=0; id<336; id++){
657 if (!fTracksEsdParam.At(id)) continue;
658 if (fClusterSatur[id]>0.3) continue; // tracks in saturation
660 if ( AcceptLaser(id) && fFitZ[id]<0.5){
668 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
672 // for (Int_t id=0; id<336; id++){
675 if (!fTracksEsdParam.At(id)) return;
676 if (!AcceptLaser(id)) return;
677 Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
680 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
681 if (!hisdz) MakeFitHistos();
682 hisdz = (TH1F*)fDeltaZ.At(id);
683 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
684 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
685 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
686 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
687 TH1F * hisSignal = (TH1F*)fSignals.At(id);
690 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
691 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
692 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
701 param->GetPxPyPz(pxyz);
703 ltrp->GetPxPyPz(lpxyz);
705 Float_t dz = param->GetZ()-ltrp->GetZ();
706 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
707 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
708 if (hisdz) hisdz->Fill(dz);
709 if (hisP3) hisP3->Fill(param->GetParameter()[3]);
710 if (hisP4) hisP4->Fill(param->GetParameter()[4]);
711 if (hisdphi) hisdphi->Fill(dphi);
712 if (hisdphiP) hisdphiP->Fill(dphiP);
713 if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
715 xhis[0] = ltrp->GetId();
716 xhis[1] = ltrp->GetSide();
717 xhis[2] = ltrp->GetRod();
718 xhis[3] = ltrp->GetBundle();
719 xhis[4] = ltrp->GetBeam();
722 xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
723 xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
724 xhis[9] = param->GetParameter()[4];
725 xhis[10]= track->GetTPCNcls();
726 xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
728 fHisLaser->Fill(xhis);
733 void AliTPCcalibLaser::FitDriftV(){
735 // Fit corrections to the drift velocity - linear approximation in the z and global y
736 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
741 z = s* (z0 - vd*(t-t0))
745 vd - nominal drift velocity
746 zs - miscalibrated position
748 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
749 vr - relative change of the drift velocity
754 zs ~ z - s*vr*(z0-s*z)+s*dzt
755 --------------------------------
756 1. Correction function vr constant:
759 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
760 dzs/dl = dz/dl +s*s*vr*dz/dl
764 const Float_t kZCut = 240; // remove the closest laser beam
765 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
766 const Float_t kDistCut = 3; // distance sigma cut
767 const Float_t kDistCutAbs = 0.25;
768 const Float_t kMinClusters = 60; // minimal amount of the clusters
769 const Float_t kMinSignal = 16; // minimal mean height of the signal
770 const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
771 static TLinearFitter fdriftA(3,"hyp2");
772 static TLinearFitter fdriftC(3,"hyp2");
773 static TLinearFitter fdriftAC(4,"hyp3");
774 TVectorD fitA(3),fitC(3),fitAC(4);
776 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
777 AliTPCParam * tpcparam = calib->GetParameters();
779 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
789 for (Int_t iter=0; iter<3; iter++){
790 fdriftA.ClearPoints();
791 fdriftC.ClearPoints();
792 fdriftAC.ClearPoints();
794 for (Int_t id=0; id<336; id++){
795 if (!fTracksEsdParam.At(id)) continue;
796 if (!AcceptLaser(id)) continue;
797 if ( fClusterSatur[id]>kSaturCut) continue;
798 if ( fClusterCounter[id]<kMinClusters) continue;
799 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
800 if (track->GetTPCsignal()<kMinSignal) continue;
801 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
802 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
809 param->GetPxPyPz(pxyz);
811 ltrp->GetPxPyPz(lpxyz);
812 if (TMath::Abs(lxyz[2])>kZCut) continue;
813 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
814 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
815 if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
816 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
819 Double_t zlength = tpcparam->GetZLength(0);
820 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
821 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
822 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
823 if (ltrp->GetSide()==0){
824 fdriftA.AddPoint(xxx,mdrift,1);
826 fdriftC.AddPoint(xxx,mdrift,1);
828 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
829 fdriftAC.AddPoint(xxx2,mdrift,1);
832 if (fdriftA.GetNpoints()>10){
835 if (iter==0) fdriftA.EvalRobust(0.7);
836 else fdriftA.EvalRobust(0.8);
837 fdriftA.GetParameters(fitA);
838 npointsA= fdriftA.GetNpoints();
839 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
840 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
841 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
842 (*fFitAside)[0] = fitA[0];
843 (*fFitAside)[1] = fitA[1];
844 (*fFitAside)[2] = fitA[2];
845 (*fFitAside)[3] = fdriftA.GetNpoints();
846 (*fFitAside)[4] = chi2A;
849 if (fdriftC.GetNpoints()>10){
851 if (iter==0) fdriftC.EvalRobust(0.7);
852 else fdriftC.EvalRobust(0.8);
854 fdriftC.GetParameters(fitC);
855 npointsC= fdriftC.GetNpoints();
856 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
857 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
858 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
859 (*fFitCside)[0] = fitC[0];
860 (*fFitCside)[1] = fitC[1];
861 (*fFitCside)[2] = fitC[2];
862 (*fFitCside)[3] = fdriftC.GetNpoints();
863 (*fFitCside)[4] = chi2C;
867 if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
869 if (iter==0) fdriftAC.EvalRobust(0.7);
870 else fdriftAC.EvalRobust(0.8);
872 fdriftAC.GetParameters(fitAC);
873 npointsAC= fdriftAC.GetNpoints();
874 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
875 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
878 for (Int_t id=0; id<336; id++){
879 if (!fTracksEsdParam.At(id)) continue;
881 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
882 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
888 param->GetPxPyPz(pxyz);
890 ltrp->GetPxPyPz(lpxyz);
891 Double_t zlength = tpcparam->GetZLength(0);
892 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
893 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
896 if (ltrp->GetSide()==0){
897 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
899 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
902 fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
907 TTreeSRedirector *cstream = GetDebugStreamer();
908 TTimeStamp tstamp(fTime);
909 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
910 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
911 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
912 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
913 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
914 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
915 TVectorD vecGoofie(20);
916 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
918 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
919 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
920 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
924 (*cstream)<<"driftv"<<
925 "run="<<fRun<< // run number
926 "event="<<fEvent<< // event number
927 "time="<<fTime<< // time stamp of event
928 "trigger="<<fTrigger<< // trigger
929 "mag="<<fMagF<< // magnetic field
930 // Environment values
931 "press0="<<valuePressure0<<
932 "press1="<<valuePressure1<<
933 "pt0="<<ptrelative0<<
934 "pt1="<<ptrelative1<<
937 "vecGoofie.="<<&vecGoofie<<
941 "driftA.="<<fFitAside<<
942 "driftC.="<<fFitCside<<
943 "driftAC.="<<fFitACside<<
955 Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
957 // Fit corrections to the drift velocity - linear approximation in the z and global y
958 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
960 // Source of outlyers :
961 // 0. Track in the saturation - postpeak
962 // 1. gating grid close the part of the signal for first bundle
964 // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
965 // 1. Robust fit is used in the itteration number 0
966 // only fraction of laser used
967 // 2. Only the tracks close to the fit used in the second itteration
971 z = s* (z0 - vd*(t-t0))
975 vd - nominal drift velocity
976 zs - miscalibrated position
978 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
979 vr - relative change of the drift velocity
984 zs ~ z - s*vr*(z0-s*z)+s*dzt
985 --------------------------------
986 1. Correction function vr constant:
989 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
990 dzs/dl = dz/dl +s*s*vr*dz/dl
993 const Int_t knLaser = 336; //n laser tracks
994 const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
996 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
997 const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
998 const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
999 const Float_t kMinClusters = 40.; // minimal amount of the clusters
1000 const Float_t kMinSignal = 2.5; // minimal mean height of the signal
1001 const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
1003 static TLinearFitter fdriftA(3,"hyp2");
1004 static TLinearFitter fdriftC(3,"hyp2");
1005 static TLinearFitter fdriftAC(4,"hyp3");
1006 TVectorD fitA(3),fitC(3),fitAC(4);
1008 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1009 AliTPCParam * tpcparam = calib->GetParameters();
1013 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
1014 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1015 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1016 for (Int_t i=0;i<5; i++){
1024 Float_t chi2AC = 10;
1028 Int_t nbA[4]={0,0,0,0};
1029 Int_t nbC[4]={0,0,0,0};
1030 TVectorD vecZM(336); // measured z potion of laser
1031 TVectorD vecA(336); // accepted laser
1032 TVectorD vecZF(336); // fitted position
1033 TVectorD vecDz(336); // deltaZ
1034 TVectorD vecZS(336); // surveyed position of laser
1035 // additional variable to cut
1036 TVectorD vecdEdx(336); // dEdx
1037 TVectorD vecSy(336); // shape y
1038 TVectorD vecSz(336); // shape z
1041 for (Int_t id=0; id<336; id++){
1043 AliTPCLaserTrack *ltrp =
1044 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1045 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1046 AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1047 vecZM(id)= (param==0) ? 0:param->GetZ();
1048 vecZS(id)= ltrp->GetZ();
1051 vecdEdx(id)=(seed)?seed->GetdEdx():0;
1052 vecSy(id) =(seed)?seed->CookShape(1):0;
1053 vecSz(id) =(seed)?seed->CookShape(2):0;
1056 if (!fTracksEsdParam.At(id)) reject|=1;
1057 if (!param) continue;
1058 if (!AcceptLaser(id)) reject|=2;
1059 if ( fClusterSatur[id]>kSaturCut) reject|=4;
1060 if ( fClusterCounter[id]<kMinClusters) reject|=8;
1062 if (reject>0) continue;
1063 if (ltrp->GetSide()==0){
1065 nbA[ltrp->GetBundle()]++;
1067 if (ltrp->GetSide()==1){
1069 nbC[ltrp->GetBundle()]++;
1073 // reject "bad events"
1078 if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1080 for (Int_t i=0;i<4;i++){
1081 //count accepted for all layers
1082 if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1083 if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1085 if (naA<3 &&naC<3) isOK=kFALSE;
1086 if (isOK==kFALSE) return kFALSE;
1091 for (Int_t iter=0; iter<2; iter++){
1092 fdriftA.ClearPoints();
1093 fdriftC.ClearPoints();
1094 fdriftAC.ClearPoints();
1095 npointsA=0; npointsC=0; npointsAC=0;
1097 for (Int_t id=0; id<336; id++){
1098 if (!fTracksEsdParam.At(id)) continue;
1099 if (!AcceptLaser(id)) continue;
1100 if ( fClusterSatur[id]>kSaturCut) continue;
1101 if ( fClusterCounter[id]<kMinClusters) continue;
1102 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1103 if (track->GetTPCsignal()<kMinSignal) continue;
1104 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1105 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1111 param->GetPxPyPz(pxyz);
1113 ltrp->GetPxPyPz(lpxyz);
1114 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1115 //if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1116 if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1117 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1119 // // drift distance
1120 // Double_t zlength = tpcparam->GetZLength(0);
1121 // Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1122 // Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1124 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1125 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1126 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1129 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1130 if (iter==0 &<rp->GetBundle()==0) continue;
1131 // skip bundle 0 - can be wrong in the 0 iteration
1132 if (ltrp->GetSide()==0){
1133 fdriftA.AddPoint(xxx,mdrift,1);
1135 fdriftC.AddPoint(xxx,mdrift,1);
1137 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1138 fdriftAC.AddPoint(xxx2,mdrift,1);
1141 if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1144 //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
1145 npointsA= fdriftA.GetNpoints();
1146 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1147 fdriftA.EvalRobust(kFraction[iter]);
1148 fdriftA.GetParameters(fitA);
1149 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1150 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1151 (*fFitAside)[0] = fitA[0];
1152 (*fFitAside)[1] = fitA[1];
1153 (*fFitAside)[2] = fitA[2];
1154 (*fFitAside)[3] = fdriftA.GetNpoints();
1155 (*fFitAside)[4] = chi2A;
1158 if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1160 //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
1161 npointsC= fdriftC.GetNpoints();
1162 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1163 fdriftC.EvalRobust(kFraction[iter]);
1164 fdriftC.GetParameters(fitC);
1165 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1166 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1167 (*fFitCside)[0] = fitC[0];
1168 (*fFitCside)[1] = fitC[1];
1169 (*fFitCside)[2] = fitC[2];
1170 (*fFitCside)[3] = fdriftC.GetNpoints();
1171 (*fFitCside)[4] = chi2C;
1175 if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1177 //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
1178 npointsAC= fdriftAC.GetNpoints();
1179 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1180 fdriftAC.EvalRobust(kFraction[iter]);
1181 fdriftAC.GetParameters(fitAC);
1182 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1183 (*fFitACside)[0] = fitAC[0];
1184 (*fFitACside)[1] = fitAC[1];
1185 (*fFitACside)[2] = fitAC[2];
1186 (*fFitACside)[3] = fdriftAC.GetNpoints();
1187 (*fFitACside)[4] = chi2AC;
1190 for (Int_t id=0; id<336; id++){
1191 if (!fTracksEsdParam.At(id)) continue;
1193 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1194 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1200 param->GetPxPyPz(pxyz);
1202 ltrp->GetPxPyPz(lpxyz);
1203 //Double_t zlength = tpcparam->GetZLength(0);
1204 //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1205 //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1206 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1207 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1208 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1212 if (ltrp->GetSide()==0){
1213 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1215 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1218 //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1220 fFitZ[id]=mdrift-fz;
1222 vecDz[id]=mdrift-fz;
1224 if (fStreamLevel>0){
1225 TTreeSRedirector *cstream = GetDebugStreamer();
1226 TTimeStamp tstamp(fTime);
1227 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1228 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1229 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1230 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1231 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1232 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1233 TVectorD vecGoofie(20);
1234 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1236 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1237 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1238 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1242 (*cstream)<<"driftvN"<<
1243 "run="<<fRun<< // run number
1244 "event="<<fEvent<< // event number
1245 "time="<<fTime<< // time stamp of event
1246 "trigger="<<fTrigger<< // trigger
1247 "mag="<<fMagF<< // magnetic field
1248 // Environment values
1249 "press0="<<valuePressure0<<
1250 "press1="<<valuePressure1<<
1251 "pt0="<<ptrelative0<<
1252 "pt1="<<ptrelative1<<
1255 "vecGoofie.="<<&vecGoofie<<
1258 "vecZM.="<<&vecZM<< // measured z position
1259 "vecZS.="<<&vecZS<< // surveyed z position
1260 "vecZF.="<<&vecZF<< // fitted z position
1261 "vecDz.="<<&vecDz<< // fitted z position
1262 "vecA.="<<&vecA<< // accept laser flag
1263 "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
1264 "vecSy.="<<&vecSy<< // shape y - to cut on
1265 "vecSz.="<<&vecSz<< // shape z - to cut on
1268 "driftA.="<<fFitAside<<
1269 "driftC.="<<fFitCside<<
1270 "driftAC.="<<fFitACside<<
1280 variables to check in debug mode:
1282 chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1283 chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1284 chainDriftN->SetAlias("driftF","vecZF.fElements");
1285 chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
1286 TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1287 TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1296 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1298 // get distance between mirror and laser track
1309 dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1312 // apply drift correction if already exist
1315 if (ltrp->GetSide()==0){
1316 if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1318 if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1320 if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1321 dist+=TMath::Abs(dz);
1323 // phi dist - divergence on 50 cm
1325 dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1330 Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
1335 TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1336 TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1337 TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1338 TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1340 TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1342 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1343 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1344 //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1349 if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1350 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
1351 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
1352 if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
1359 Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1361 // Find corresponding mirror
1362 // add the corresponding tracks
1364 if (!track->GetOuterParam()) return -1;
1366 Float_t kRadius0 = 252;
1367 Float_t kRadius = 254.2;
1369 Float_t counterSatur=0;
1372 for (Int_t irow=158;irow>-1;--irow) {
1373 AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1375 Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1376 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1377 if (pedgeY<3) continue;
1378 if (pedgeX<3) continue;
1380 if (c->GetDetector()%36<18) csideA++;
1381 if (c->GetDetector()%36>=18) csideC++;
1382 if (c->GetMax()>900) counterSatur++;
1384 counterSatur/=(countercl+1);
1388 if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
1391 if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1394 AliExternalTrackParam param(*(track->GetOuterParam()));
1395 AliTracker::PropagateTrackTo(¶m,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
1396 AliTracker::PropagateTrackTo(¶m,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
1397 AliTPCLaserTrack ltr;
1398 AliTPCLaserTrack *ltrp=0x0;
1399 // AliTPCLaserTrack *ltrpjw=0x0;
1401 Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
1402 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
1403 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1405 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1406 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1410 if (id<0) return -1;
1411 if (ltrp->GetSide()!=side) return -1;
1417 if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1420 Float_t radius=TMath::Abs(ltrp->GetX());
1421 param.Rotate(ltrp->GetAlpha());
1422 AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
1424 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1425 Bool_t accept=kTRUE;
1427 // choose closer track
1429 AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1431 Float_t dist0=GetDistance(param0,ltrp);
1432 Float_t dist1=GetDistance(¶m,ltrp);
1433 if (dist0<dist1) accept=kFALSE;
1437 fClusterCounter[id]=countercl;
1438 fTracksEsdParam.AddAt(param.Clone(),id);
1439 fTracksEsd.AddAt(track,id);
1440 fTracksTPC.AddAt(seed,id);
1447 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1449 // Dump Laser info to the tree
1451 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1452 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1453 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1462 param->GetPxPyPz(pxyz);
1464 ltrp->GetPxPyPz(lpxyz);
1465 Float_t dist3D = GetDistance(param,ltrp);
1466 Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1467 Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1470 if (fStreamLevel>0){
1471 TTreeSRedirector *cstream = GetDebugStreamer();
1472 Int_t time = fESD->GetTimeStamp();
1473 Bool_t accept = AcceptLaser(id);
1475 (*cstream)<<"Track"<<
1477 "run="<<fRun<< // run number
1478 "event="<<fEvent<< // event number
1479 "time="<<fTime<< // time stamp of event
1480 "trigger="<<fTrigger<< // trigger
1481 "mag="<<fMagF<< // magnetic field
1485 "driftA.="<<fFitAside<<
1486 "driftC.="<<fFitCside<<
1490 "distphi="<<distphi<<
1493 "counter="<<fCounter[id]<<
1494 "clcounter="<<fClusterCounter[id]<<
1495 "satur="<<fClusterSatur[id]<<
1519 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1521 // Refit the track with different tracklet models:
1522 // 1. Per ROC using the kalman filter, different edge cuts
1523 // 2. Per ROC linear in y and z
1524 // 3. Per ROC quadratic in y and z
1525 // 4. Per track offset for each sector, linear for each sector, common quadratic
1526 // store x, y, z information for all models and the cluster to calculate the residuals
1529 // The clusters which do not fulfill given criteria are skipped
1531 // Cluters removed from fit
1532 // 1. Extended shape > kShapeCut
1533 // 2. In saturation Max > kMax
1534 // 3. Distance to edge < cutEdge
1536 // Clusters not used for the calibration:
1538 // 1. Extended shape > kShapeCut
1539 // 2. In saturation Max > kMax
1542 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1543 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1544 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1546 AliTPCclusterMI dummyCl;
1549 Int_t kMaxTracklets=2;
1550 Float_t kShapeCut = 1.3;
1551 Float_t kRatioCut = 2.;
1553 Float_t kMax = 900.;
1556 //=============================================//
1557 // Linear Fitters for the Different Approaches //
1558 //=============================================//
1559 //linear fit model in y and z; inner - outer sector, combined with offset
1560 static TLinearFitter fy1I(2,"hyp1");
1561 static TLinearFitter fy1O(2,"hyp1");
1562 static TLinearFitter fz1I(2,"hyp1");
1563 static TLinearFitter fz1O(2,"hyp1");
1564 static TLinearFitter fy1IO(3,"hyp2");
1565 static TLinearFitter fz1IO(3,"hyp2");
1566 //quadratic fit model in y and z; inner - sector
1567 static TLinearFitter fy2I(3,"hyp2");
1568 static TLinearFitter fy2O(3,"hyp2");
1569 static TLinearFitter fz2I(3,"hyp2");
1570 static TLinearFitter fz2O(3,"hyp2");
1571 //common quadratic fit for IROC and OROC in y and z
1572 static TLinearFitter fy4(5,"hyp4");
1573 static TLinearFitter fz4(5,"hyp4");
1583 //=============================//
1584 // Loop over all Tracklet Cuts //
1585 //=============================//
1586 for (Int_t icut=0; icut<fNcuts; icut++){
1587 Float_t xinMin = 2500, xinMax=-90;
1588 Float_t xoutMin=2500, xoutMax=-90;
1589 Float_t msigmaYIn=0;
1590 Float_t msigmaYOut=0;
1591 Float_t mqratioIn=0;
1592 Float_t mqratioOut=0;
1595 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1597 Double_t edgeCutX = fEdgeXcuts[icut];
1598 Double_t edgeCutY = fEdgeYcuts[icut];
1599 Int_t nclCut = (Int_t)fNClCuts[icut];
1600 //=========================//
1601 // Parameters to calculate //
1602 //=========================//
1603 //fit parameter inner, outer tracklet and combined fit
1604 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1605 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1607 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1608 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1609 TVectorD vecy4res(5),vecz4res(5);
1610 TVectorD vecy1resIO(3),vecz1resIO(3);
1611 // cluster and track positions for each row - used for residuals
1612 TVectorD vecgX(159); // global X
1613 TVectorD vecgY(159); // global Y
1614 TVectorD vecgZ(159); // global Z
1616 TVectorD vecX(159); // x is the same for all (row center)
1617 TVectorD vecYkalman(159); // y from kalman fit
1618 TVectorD vecZkalman(159); // z from kalman fit
1619 TVectorD vecY1(159); // y from pol1 fit per ROC
1620 TVectorD vecZ1(159); // z from pol1 fit per ROC
1621 TVectorD vecY1IO(159); // y from pol1 fit per ROC
1622 TVectorD vecZ1IO(159); // z from pol1 fit per ROC
1623 TVectorD vecY2(159); // y from pol2 fit per ROC
1624 TVectorD vecZ2(159); // z from pol2 fit per ROC
1625 TVectorD vecY4(159); // y from sector fit
1626 TVectorD vecZ4(159); // z from sector fit
1627 TVectorD vecClY(159); // y cluster position
1628 TVectorD vecClZ(159); // z cluster position
1629 TVectorD vecSec(159); // sector for each row
1630 TVectorD isReject(159); // flag - cluster to be rejected
1632 Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
1633 Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
1634 Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
1635 Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
1636 Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
1637 Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
1638 Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
1639 Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
1640 Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
1641 Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
1642 Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
1643 Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
1645 Int_t innerSector = -1; // number of inner sector
1646 Int_t outerSector = -1; // number of outer sector
1647 Int_t nclI=0; // number of clusters (inner)
1648 Int_t nclO=0; // number of clusters (outer)
1649 //=================================================//
1650 // Perform the Kalman Fit using the Tracklet Class //
1651 //=================================================//
1652 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1653 TObjArray tracklets=
1654 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1655 kFALSE,nclCut,kMaxTracklets);
1656 // tracklet pointers
1657 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1658 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1659 AliTPCTracklet *tr=0x0;
1660 AliTPCTracklet dummy;
1661 //continue if we didn't find a tracklet
1662 if ( !trInner && !trOuter ) continue;
1663 //================================================================================//
1664 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1665 //================================================================================//
1666 if ( trInner && trOuter ){
1667 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1668 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1672 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1676 if ( !trInner->GetInner() ) continue;
1678 if ( trInner->GetSector()>35 ){
1683 if ( !trOuter->GetInner() ) continue;
1685 if ( trOuter->GetSector()<36 ){
1691 innerSector = trInner->GetSector();
1692 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1693 outerSector = trOuter->GetSector();
1694 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1696 // array of clusters
1697 TClonesArray arrCl("AliTPCclusterMI",159);
1698 arrCl.ExpandCreateFast(159);
1699 //=======================================//
1700 // fill fitters with cluster information //
1701 //=======================================//
1702 AliDebug(3,"Filing Cluster Information");
1707 for (Int_t irow=158;irow>-1;--irow) {
1708 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1709 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1714 Float_t meanY=0, sumY=0;
1715 for (Int_t drow=-1;drow<=1;drow++) {
1716 if (irow+drow<0) continue;
1717 if (irow+drow>158) continue;
1718 AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1719 if (!ccurrent) continue;
1720 Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1721 if ( roc!=innerSector && roc!=outerSector ) continue;
1722 if (ccurrent->GetX()<10) continue;
1723 if (ccurrent->GetY()==0) continue;
1724 meanY+=ccurrent->GetY();
1727 if (sumY>0) meanY/=sumY;
1732 Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1733 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1736 Int_t roc = static_cast<Int_t>(c->GetDetector());
1737 if ( roc!=innerSector && roc!=outerSector ) continue;
1739 //store clusters in clones array
1742 if (c->GetMax()<4) continue; // noise cluster?
1743 if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
1745 vecX[irow] = c->GetX();
1746 vecClY[irow] = c->GetY();
1747 vecClZ[irow] = c->GetZ();
1750 // c->GetGlobalXYZ(gxyz);
1751 // vecgX[irow] = gxyz[0];
1752 // vecgY[irow] = gxyz[1];
1753 // vecgZ[irow] = gxyz[2];
1755 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1756 Double_t y = vecClY[irow];
1757 Double_t z = vecClZ[irow];
1759 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
1760 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1761 Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
1762 if ( roc == innerSector ) {
1763 x4[0]=1; //offset inner - outer sector
1764 x4[1]=x; //slope parameter inner sector
1767 x4[2]=x; //slope parameter outer sector
1769 x4[3]=x*x; //common parabolic shape
1770 if (pedgeX < fEdgeXcuts[icut]) continue;
1771 if (pedgeY < fEdgeYcuts[icut]) continue;
1772 if (c->GetMax()>900) continue; // cluster in saturation
1774 if ( roc==innerSector ){
1775 fy1I.AddPoint(x2,y);
1776 fz1I.AddPoint(x2,z);
1777 fy2I.AddPoint(x2,y);
1778 fz2I.AddPoint(x2,z);
1780 if (c->GetX()<xinMin) xinMin=c->GetX();
1781 if (c->GetX()>xinMax) xinMax=c->GetX();
1783 msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1784 mqratioIn +=c->GetMax()/c->GetQ();
1787 if ( roc==outerSector ){
1788 fy1O.AddPoint(x2,y);
1789 fz1O.AddPoint(x2,z);
1790 fy2O.AddPoint(x2,y);
1791 fz2O.AddPoint(x2,z);
1793 if (c->GetX()<xoutMin) xoutMin=c->GetX();
1794 if (c->GetX()>xoutMax) xoutMax=c->GetX();
1795 msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1796 mqratioOut +=c->GetMax()/c->GetQ();
1801 fy1IO.AddPoint(xIO,y);
1802 fz1IO.AddPoint(xIO,z);
1812 //======================================//
1813 // Evaluate and retrieve fit parameters //
1814 //======================================//
1815 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1817 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1822 fy1I.GetParameters(vecy1resInner);
1823 fz1I.GetParameters(vecz1resInner);
1824 fy2I.GetParameters(vecy2resInner);
1825 fz2I.GetParameters(vecz2resInner);
1826 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1827 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1828 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1829 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1832 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1837 fy1O.GetParameters(vecy1resOuter);
1838 fz1O.GetParameters(vecz1resOuter);
1839 fy2O.GetParameters(vecy2resOuter);
1840 fz2O.GetParameters(vecz2resOuter);
1841 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1842 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1843 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1844 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1847 if ( innerSector>0 && outerSector>0 ){
1848 if (fy4.GetNpoints()>0) {
1850 fy4.GetParameters(vecy4res);
1851 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1853 if (fz4.GetNpoints()>0) {
1855 fz4.GetParameters(vecz4res);
1856 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1858 if (fy1IO.GetNpoints()>0) {
1860 fy1IO.GetParameters(vecy1resIO);
1861 chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1863 if (fz1IO.GetNpoints()>0) {
1865 fz1IO.GetParameters(vecz1resIO);
1866 chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1870 fy4.ClearPoints(); fz4.ClearPoints();
1871 fy1I.ClearPoints(); fy1O.ClearPoints();
1872 fz1I.ClearPoints(); fz1O.ClearPoints();
1873 fy2I.ClearPoints(); fy2O.ClearPoints();
1874 fz2I.ClearPoints(); fz2O.ClearPoints();
1875 fy1IO.ClearPoints(); fz1IO.ClearPoints();
1876 //==============================//
1877 // calculate tracklet positions //
1878 //==============================//
1879 AliDebug(4,"Calculate tracklet positions");
1880 for (Int_t irow=158;irow>-1;--irow) {
1882 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1883 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1887 if (!c) { //no cluster
1890 if (c->GetMax()>kMax){ //saturation
1893 if ( vecSec[irow] == outerSector ) { //extended shape
1894 if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
1895 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
1897 if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
1898 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
1904 if ( vecSec[irow]==-1 ) continue; //no cluster info
1905 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1907 Double_t x=vecX[irow];
1908 Double_t xref=x-133.4;
1910 Double_t yoffInner=0;
1911 Double_t zoffInner=0;
1912 Double_t yoffInner1=0;
1913 Double_t zoffInner1=0;
1914 Double_t yslopeInner=0;
1915 Double_t yslopeOuter=0;
1916 Double_t zslopeInner=0;
1917 Double_t zslopeOuter=0;
1918 //positions of hyperplane fits
1919 if ( vecSec[irow] == outerSector ) {
1921 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1922 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1923 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1924 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1925 yslopeOuter=vecy4res[3];
1926 zslopeOuter=vecz4res[3];
1929 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1930 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1931 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1932 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1933 yoffInner=vecy4res[1];
1934 zoffInner=vecz4res[1];
1935 yoffInner1=vecy1resIO[1];
1936 zoffInner1=vecz1resIO[1];
1937 yslopeInner=vecy4res[2];
1938 zslopeInner=vecz4res[2];
1940 vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1941 vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1942 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1943 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1944 //positions of kalman fits
1945 Double_t gxyz[3],xyz[3];
1946 AliExternalTrackParam *param = 0x0;
1948 param=tr->GetInner();
1950 param->GetXYZ(gxyz);
1951 Float_t bz = AliTracker::GetBz(gxyz);
1952 param->GetYAt(x, bz, xyz[1]);
1953 param->GetZAt(x, bz, xyz[2]);
1954 vecYkalman[irow]=xyz[1];
1955 vecZkalman[irow]=xyz[2];
1962 //=====================================================================//
1963 // write results from the different tracklet fits with debug streamers //
1964 //=====================================================================//
1965 if (fStreamLevel>4){
1966 TTreeSRedirector *cstream = GetDebugStreamer();
1968 Float_t dedx = track->GetdEdx();
1969 (*cstream)<<"FitModels"<<
1970 "run="<<fRun<< // run number
1971 "event="<<fEvent<< // event number
1972 "time="<<fTime<< // time stamp of event
1973 "trigger="<<fTrigger<< // trigger
1974 "mag="<<fMagF<< // magnetic field
1977 "edgeCutX=" << edgeCutX <<
1978 "edgeCutY=" << edgeCutY <<
1979 "nclCut=" << nclCut <<
1980 "innerSector="<< innerSector <<
1981 "outerSector="<< outerSector <<
1984 "Tr.=" << extparam <<
1985 "yPol1In.=" << &vecy1resInner <<
1986 "zPol1In.=" << &vecz1resInner <<
1987 "yPol1InOut.="<< &vecy1resIO <<
1988 "zPol1InOut.="<< &vecz1resIO <<
1989 "yPol2In.=" << &vecy2resInner <<
1990 "zPol2In.=" << &vecz2resInner <<
1991 "yPol1Out.=" << &vecy1resOuter <<
1992 "zPol1Out.=" << &vecz1resOuter <<
1993 "yPol2Out.=" << &vecy2resOuter <<
1994 "zPol2Out.=" << &vecz2resOuter <<
1995 "yInOut.=" << &vecy4res <<
1996 "zInOut.=" << &vecz4res <<
1997 "chi2y1In=" << chi2I1y <<
1998 "chi2z1In=" << chi2I1z <<
1999 "chi2y1InOut="<< chi2IO1y <<
2000 "chi2z1InOut="<< chi2IO1z <<
2001 "chi2y1Out=" << chi2O1y <<
2002 "chi2z1Out=" << chi2O1z <<
2003 "chi2y2In=" << chi2I2y <<
2004 "chi2z2In=" << chi2I2z <<
2005 "chi2y2Out=" << chi2O2y <<
2006 "chi2z2Out=" << chi2O2z <<
2007 "chi2yInOut=" << chi2IOy <<
2008 "chi2zInOut=" << chi2IOz <<
2009 "trletIn.=" << trInner <<
2010 "trletOut.=" << trOuter <<
2013 "xinMin=" << xinMin<<
2014 "xinMax=" << xinMax<<
2015 "xoutMin=" << xoutMin<<
2016 "xoutMax=" << xoutMax<<
2017 "msigmaYIn=" <<msigmaYIn<<
2018 "msigmaYOut=" <<msigmaYOut<<
2019 "mqratioIn=" <<mqratioIn<<
2020 "mqratioOut=" << mqratioOut <<
2025 // wirte residuals information
2026 if (fStreamLevel>5){
2027 TTreeSRedirector *cstream = GetDebugStreamer();
2029 Float_t dedx = track->GetdEdx();
2030 (*cstream)<<"Residuals"<<
2031 "run="<<fRun<< // run number
2032 "event="<<fEvent<< // event number
2033 "time="<<fTime<< // time stamp of event
2034 "trigger="<<fTrigger<< // trigger
2035 "mag="<<fMagF<< // magnetic field
2038 "edgeCutX=" << edgeCutX <<
2039 "edgeCutY=" << edgeCutY <<
2040 "nclCut=" << nclCut <<
2042 "Tr.=" << extparam<<
2045 "vX.=" << &vecgX<< // global x
2046 "vY.=" << &vecgY<< // global y
2047 "vZ.=" << &vecgZ<< // global z
2049 "TrYpol1.=" << &vecY1 <<
2050 "TrZpol1.=" << &vecZ1 <<
2051 "TrYpol2.=" << &vecY2 <<
2052 "TrZpol2.=" << &vecZ2 <<
2053 "TrYpol1InOut.="<< &vecY1IO <<
2054 "TrZpol1InOut.="<< &vecZ1IO <<
2055 "TrYInOut.=" << &vecY4 <<
2056 "TrZInOut.=" << &vecZ4 <<
2057 "ClY.=" << &vecClY <<
2058 "ClZ.=" << &vecClZ <<
2059 "isReject.=" << &isReject<<
2060 "sec.=" << &vecSec <<
2063 "xinMin=" << xinMin<<
2064 "xinMax=" << xinMax<<
2065 "xoutMin=" << xoutMin<<
2066 "xoutMax=" << xoutMax<<
2067 "msigmaYIn=" <<msigmaYIn<<
2068 "msigmaYOut=" <<msigmaYOut<<
2069 "mqratioIn=" <<mqratioIn<<
2070 "mqratioOut=" << mqratioOut <<
2071 "yInOut.=" << &vecy4res <<
2072 "zInOut.=" << &vecz4res <<
2074 "chi2y1In=" << chi2I1y << //
2075 "chi2z1In=" << chi2I1z <<
2076 "chi2y1Out=" << chi2O1y <<
2077 "chi2z1Out=" << chi2O1z <<
2078 "chi2y1InOut="<< chi2IO1y <<
2079 "chi2z1InOut="<< chi2IO1z <<
2080 "chi2y2In=" << chi2I2y <<
2081 "chi2z2In=" << chi2I2z <<
2082 "chi2y2Out=" << chi2O2y <<
2083 "chi2z2Out=" << chi2O2z <<
2084 "chi2yInOut=" << chi2IOy <<
2085 "chi2zInOut=" << chi2IOz <<
2087 "yPol1In.=" << &vecy1resInner <<
2088 "zPol1In.=" << &vecz1resInner <<
2089 "yPol2In.=" << &vecy2resInner <<
2090 "zPol2In.=" << &vecz2resInner <<
2091 "yPol1Out.=" << &vecy1resOuter <<
2092 "zPol1Out.=" << &vecz1resOuter <<
2093 "yPol1InOut.="<< &vecy1resIO <<
2094 "zPol1InOut.="<< &vecz1resIO <<
2095 "yPol2Out.=" << &vecy2resOuter <<
2096 "zPol2Out.=" << &vecz2resOuter <<
2102 //==========================//
2103 // Fill Residual Histograms //
2104 //==========================//
2105 if (!fHisNclIn) MakeFitHistos();
2107 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2108 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2109 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2110 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2111 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2112 TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2113 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2114 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2116 for (Int_t irow=158;irow>-1;--irow) {
2117 if (vecSec[irow]==-1)continue; // no cluster info
2118 if (isReject[irow]>0.5) continue; //
2119 Double_t ycl = vecClY[irow];
2120 Double_t yfit = vecY1[irow];
2121 Double_t yfit2 = vecY2[irow];
2122 Double_t x = vecX[irow];
2123 Double_t yabsbeam = -1000;
2124 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2125 yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2126 else if(innerSector==fBeamSectorInner[id])
2127 yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2129 // Double_t yfit3 = vecY2[irow];
2130 Double_t zcl = vecClZ[irow];
2131 Double_t zfit = vecZ1[irow];
2132 Double_t zfit2 = vecZ2[irow];
2133 //Double_t zfit3 = vecZ2[irow];
2136 // The expressions for zcorrected has been obtained by
2137 // inverting the fits in the FitDriftV() method (ignoring the
2138 // global y dependence for now):
2140 // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2142 // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2145 // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
2147 // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2149 Double_t dzabs = -1000;
2150 Double_t zcorrected = -1000;
2151 if (ltrp->GetSide()==0){
2152 if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2153 // ignore global y dependence for now
2155 if(!fUseFixedDriftV)
2156 zcorrected = (zcl + (*fFitAside)[0] -
2157 (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2159 zcorrected = (zcl + fFixedFitAside0 -
2160 (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2161 // zcorrected = zcl;
2162 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2163 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2164 else if(innerSector==fBeamSectorInner[id])
2165 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2168 if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2170 if(!fUseFixedDriftV)
2171 zcorrected = (zcl - (*fFitCside)[0] +
2172 (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2174 zcorrected = (zcl - fFixedFitCside0 +
2175 (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2177 // zcorrected = zcl;
2178 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2179 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2180 else if(innerSector==fBeamSectorInner[id])
2181 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2185 if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2187 profy->Fill(irow,ycl-yfit);
2188 profy2->Fill(irow,ycl-yfit2);
2190 fHisYAbsErrors->Fill(id);
2191 // profyabs->Fill(irow,-0.99);
2193 profyabs->Fill(irow,ycl-yabsbeam);
2195 // profy3->Fill(irow,ycl-yfit3);
2198 profz->Fill(irow,zcl-zfit);
2199 profz2->Fill(irow,zcl-zfit2);
2200 //profz3->Fill(irow,zcl-zfit3);
2203 fHisZAbsErrors->Fill(id);
2205 profzabs->Fill(irow,dzabs);
2211 // Fill laser fit histograms
2213 Float_t dedx = track->GetdEdx();
2214 if (nclI>20&&nclO>20){
2215 fHisNclIn->Fill(id,nclI); //->Number of clusters inner
2216 fHisNclOut->Fill(id,nclO); //->Number of clusters outer
2217 fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
2219 fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
2220 fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
2221 fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
2223 fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2224 fHisdZfit->Fill(id,fFitZ[id]);
2227 fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
2228 fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
2229 fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
2230 fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
2231 fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
2234 fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
2235 fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
2236 fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
2237 fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
2238 fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
2241 fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
2242 fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
2243 fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
2245 fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
2246 fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
2247 fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
2249 fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
2250 fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
2251 fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
2252 fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
2254 fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
2255 fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
2256 fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
2257 fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
2260 fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
2261 fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
2265 fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
2266 fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
2271 // Fill raw THnSparses
2273 for (Int_t irow=0;irow<159;irow++) {
2274 AliTPCclusterMI *c=track->GetClusterPointer(irow);
2276 if (c->GetMax()>800) continue; // saturation cut
2277 //if (TMath::Sqrt(TMath::Abs(c->GetSigmaY2()))>1) continue; // saturation cut
2279 Double_t deltaY=c->GetY()-(*ltrp->GetVecLY())[irow];
2280 Double_t deltaZ=c->GetZ()-(*ltrp->GetVecLZ())[irow];
2281 //TString axisName[6]={"Delta","bin", "rms shape", "Q", "row","trackID"}
2282 Double_t xyz[6]={0, 0, 0,TMath::Sqrt(c->GetMax()),irow,id};
2285 xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaY2()));
2286 fHisLaserPad->Fill(xyz);
2288 xyz[1]=c->GetTimeBin();
2289 xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaZ2()));
2290 fHisLaserTime->Fill(xyz);
2296 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2298 // Dump information about laser beams
2299 // isOK variable indicates usability of the beam
2300 // Beam is not usable if:
2301 // a. No entries in range (krmsCut0)
2302 // b. Big sperad (krmscut1)
2303 // c. RMSto fit sigma bigger then (kmultiCut)
2304 // d. Too big angular spread
2307 const Float_t krmsCut0=0.001;
2308 const Float_t krmsCut1=0.16;
2309 const Float_t kmultiCut=2;
2310 const Float_t kcutP0=0.002;
2311 AliMagF* magF= dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
2312 Double_t xyz[3]={90,0,10}; // tmp. global position
2313 Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical
2314 Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local
2316 AliTPCcalibLaser *laser = this;
2317 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2318 TF1 fg("fg","gaus");
2319 AliTPCParam * tpcparam = 0;
2320 // start set up for absolute residuals analysis
2322 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2323 tpcparam = calib->GetParameters();
2324 if (!tpcparam) tpcparam = new AliTPCParamSR;
2326 AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2328 Float_t bfield = 0, bz=0;
2331 Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
2332 current = grp->GetL3Current((AliGRPObject::Stats)0);
2333 bfield = polarity*5*current/30000.;
2334 bz = polarity*5*current/30000.;
2335 printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2338 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2339 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2340 TLinearFitter lfabsyInner(2);
2341 lfabsyInner.SetFormula("1 ++ x");
2342 TLinearFitter lfabszInner(2);
2343 lfabszInner.SetFormula("1 ++ x");
2345 TLinearFitter lfabsyOuter(2);
2346 lfabsyOuter.SetFormula("1 ++ x");
2347 TLinearFitter lfabszOuter(2);
2348 lfabszOuter.SetFormula("1 ++ x");
2349 // end set up for absolute residuals analysis
2353 for (Int_t id=0; id<336; id++){
2355 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
2356 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2357 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
2358 TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
2359 TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
2360 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2361 //if (!hisphi) continue;
2362 Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2363 //if (entries<minEntries) continue;
2365 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2367 AliTPCLaserTrack::LoadTracks();
2368 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2370 ltrp->UpdatePoints();
2371 pcstream->GetFile()->cd();
2372 if (hisphi) hisphi->Write();
2373 if (hisphiP) hisphiP->Write();
2374 if (hisZ) hisZ->Write();
2375 if (hisP3) hisP3->Write();
2376 if (hisP4) hisP4->Write();
2378 Float_t meanphi = hisphi->GetMean();
2379 Float_t rmsphi = hisphi->GetRMS();
2381 Float_t meanphiP = hisphiP->GetMean();
2382 Float_t rmsphiP = hisphiP->GetRMS();
2383 Float_t meanZ = hisZ->GetMean();
2384 Float_t rmsZ = hisZ->GetRMS();
2385 if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2386 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2387 Double_t gphi1 = fg.GetParameter(1);
2388 Double_t gphi2 = fg.GetParameter(2);
2389 if (hisphiP->GetRMS()>0)
2390 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2391 Double_t gphiP1 = fg.GetParameter(1);
2392 Double_t gphiP2 = fg.GetParameter(2);
2394 if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2395 hisZ->Fit(&fg,"","");
2396 Double_t gz1 = fg.GetParameter(1);
2397 Double_t gz2 = fg.GetParameter(2);
2399 if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2400 hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2401 Double_t gp31 = fg.GetParameter(1);
2402 Double_t gp32 = fg.GetParameter(2);
2403 Double_t meanp3 = hisP3->GetMean();
2404 Double_t rmsp3 = hisP3->GetRMS();
2406 if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2407 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2408 Double_t gp41 = fg.GetParameter(1);
2409 Double_t gp42 = fg.GetParameter(2);
2410 Double_t meanp4 = hisP4->GetMean();
2411 Double_t rmsp4 = hisP4->GetRMS();
2413 Float_t meanS=hisS->GetMean();
2418 ltrp->GetPxPyPz(lpxyz);
2420 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2421 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2422 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2423 if (gphiP2>kcutP0) isOK=kFALSE;
2429 his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2430 Float_t mnclIn = his->GetMean();
2432 his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2433 Float_t mnclOut = his->GetMean();
2435 his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2436 Float_t mnclIO = his->GetMean();
2438 his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2439 Float_t mLclIn = his->GetMean();
2441 his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2442 Float_t mLclOut = his->GetMean();
2444 his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2445 Float_t mLclIO = his->GetMean();
2448 his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2449 Float_t mdEdx = his->GetMean();
2455 his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2456 Float_t mChi2YIn1= his->GetMean();
2458 his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2459 Float_t mChi2YOut1 = his->GetMean();
2461 his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2462 Float_t mChi2YIn2 = his->GetMean();
2464 his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2465 Float_t mChi2YOut2 = his->GetMean();
2467 his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2468 Float_t mChi2YIO1 = his->GetMean();
2470 his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2471 Float_t mChi2ZIn1 = his->GetMean();
2473 his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2474 Float_t mChi2ZOut1 = his->GetMean();
2476 his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2477 Float_t mChi2ZIn2 = his->GetMean();
2479 his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2480 Float_t mChi2ZOut2 = his->GetMean();
2482 his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2483 Float_t mChi2ZIO1 = his->GetMean();
2488 his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2489 Float_t edZfit = his->GetEntries();
2490 Float_t mdZfit = his->GetMean();
2491 Float_t rdZfit = his->GetRMS();
2494 his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2495 Float_t ePy1vP0 = his->GetEntries();
2496 Float_t mPy1vP0 = his->GetMean();
2497 Float_t rPy1vP0 = his->GetRMS();
2500 his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2501 Float_t ePy2vP0 = his->GetEntries();
2502 Float_t mPy2vP0 = his->GetMean();
2503 Float_t rPy2vP0 = his->GetRMS();
2506 his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2507 Float_t ePy3vP0 = his->GetEntries();
2508 Float_t mPy3vP0 = his->GetMean();
2509 Float_t rPy3vP0 = his->GetRMS();
2512 his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2513 Float_t ePy1vP1 = his->GetEntries();
2514 Float_t mPy1vP1 = his->GetMean();
2515 Float_t rPy1vP1 = his->GetRMS();
2518 his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2519 Float_t ePy2vP1 = his->GetEntries();
2520 Float_t mPy2vP1 = his->GetMean();
2521 Float_t rPy2vP1 = his->GetRMS();
2524 his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2525 Float_t ePy3vP1 = his->GetEntries();
2526 Float_t mPy3vP1 = his->GetMean();
2527 Float_t rPy3vP1 = his->GetRMS();
2530 his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2531 Float_t ePy2vP2In = his->GetEntries();
2532 Float_t mPy2vP2In = his->GetMean();
2533 Float_t rPy2vP2In = his->GetRMS();
2536 his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2537 Float_t ePy2vP2Out = his->GetEntries();
2538 Float_t mPy2vP2Out = his->GetMean();
2539 Float_t rPy2vP2Out = his->GetRMS();
2542 his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2543 Float_t ePy3vP2IO = his->GetEntries();
2544 Float_t mPy3vP2IO = his->GetMean();
2545 Float_t rPy3vP2IO = his->GetRMS();
2550 his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2551 Float_t ePz1vP0 = his->GetEntries();
2552 Float_t mPz1vP0 = his->GetMean();
2553 Float_t rPz1vP0 = his->GetRMS();
2556 his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2557 Float_t ePz2vP0 = his->GetEntries();
2558 Float_t mPz2vP0 = his->GetMean();
2559 Float_t rPz2vP0 = his->GetRMS();
2562 his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2563 Float_t ePz3vP0 = his->GetEntries();
2564 Float_t mPz3vP0 = his->GetMean();
2565 Float_t rPz3vP0 = his->GetRMS();
2568 his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2569 Float_t ePz1vP1 = his->GetEntries();
2570 Float_t mPz1vP1 = his->GetMean();
2571 Float_t rPz1vP1 = his->GetRMS();
2574 his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2575 Float_t ePz2vP1 = his->GetEntries();
2576 Float_t mPz2vP1 = his->GetMean();
2577 Float_t rPz2vP1 = his->GetRMS();
2580 his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2581 Float_t ePz3vP1 = his->GetEntries();
2582 Float_t mPz3vP1 = his->GetMean();
2583 Float_t rPz3vP1 = his->GetRMS();
2586 his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2587 Float_t ePz2vP2In = his->GetEntries();
2588 Float_t mPz2vP2In = his->GetMean();
2589 Float_t rPz2vP2In = his->GetRMS();
2592 his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2593 Float_t ePz2vP2Out = his->GetEntries();
2594 Float_t mPz2vP2Out = his->GetMean();
2595 Float_t rPz2vP2Out = his->GetRMS();
2598 his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2599 Float_t ePz3vP2IO = his->GetEntries();
2600 Float_t mPz3vP2IO = his->GetMean();
2601 Float_t rPz3vP2IO = his->GetRMS();
2604 // Fit absolute laser residuals
2605 TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2606 TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2608 Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2609 Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2611 TVectorD vecX(159); // X
2612 TVectorD vecY(159); // Y
2613 TVectorD vecR(159); // R
2614 TVectorD vecDY(159); // absolute residuals in Y
2615 TVectorD vecDZ(159); // absolute residuals in Z
2616 TVectorD vecN(159); // number of clusters
2617 TVectorD vecEy(159); //error y
2618 TVectorD vecEz(159); //error z
2619 TVectorD vecPhi(159); // local tangent
2620 TVectorD vecPhiR(159); // local tangent
2621 // magnetic field integrals
2622 TVectorD vecIBR(159); // radial
2623 TVectorD vecIBRPhi(159); // r-phi
2624 TVectorD vecIBLX(159); // local x
2625 TVectorD vecIBLY(159); // local y
2626 TVectorD vecIBGX(159); // local x
2627 TVectorD vecIBGY(159); // local y
2628 TVectorD vecIBZ(159); // z
2630 for (Int_t irow=0;irow<159;irow++){
2638 Double_t gx =(*(ltrp->fVecGX))[irow];
2639 Double_t gy =(*(ltrp->fVecGY))[irow];
2640 Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
2641 Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2642 Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2643 xyz[2]=(*(ltrp->fVecGZ))[irow];
2644 xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2645 xyz[1]=TMath::ATan2(gy,gx);
2646 Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2648 magF->GetTPCIntCyl(xyz,bxyz);
2649 magF->GetTPCInt(gxyz,bgxyz);
2650 vecIBR[irow]=bxyz[0];
2651 vecIBRPhi[irow]=bxyz[1];
2653 vecIBGX[irow]=bgxyz[0];
2654 vecIBGY[irow]=bgxyz[1];
2656 vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
2657 vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2660 vecIBZ[irow]=bxyz[2];
2665 lfabsyInner.ClearPoints();
2666 lfabszInner.ClearPoints();
2667 lfabsyOuter.ClearPoints();
2668 lfabszOuter.ClearPoints();
2671 Float_t yAbsInOffset = -100;
2672 Float_t yAbsInSlope = -100;
2673 Float_t yAbsInDeltay = -100;
2675 Float_t zAbsInOffset = -100;
2676 Float_t zAbsInSlope = -100;
2677 Float_t zAbsInDeltaz = -100;
2679 Float_t yAbsOutOffset = -100;
2680 Float_t yAbsOutSlope = -100;
2681 Float_t yAbsOutDeltay = -100;
2683 Float_t zAbsOutOffset = -100;
2684 Float_t zAbsOutSlope = -100;
2685 Float_t zAbsOutDeltaz = -100;
2687 Float_t lasTanPhiLocIn = -100;
2688 Float_t lasTanPhiLocOut = -100;
2690 if(histAbsY && histAbsY->GetEntries()>0) {
2692 Double_t rotAngOut = 10;
2693 Double_t rotAngIn = 10;
2694 if((secInner%36)!=(secOuter%36))
2695 rotAngIn += 20; // 30 degrees
2697 // Calculate laser mirror X position in local frame
2698 Double_t laserposOut =
2699 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2700 Double_t laserposIn =
2701 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2703 // Calculate laser tan phi in local frame
2704 lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2705 if(lasTanPhiLocOut<0) {
2706 lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2707 lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2710 lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2711 lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2714 lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2715 lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2717 TProfile* yprof = histAbsY->ProfileX("yprof");
2718 TProfile* zprof = histAbsZ->ProfileX("zprof");
2720 for(Int_t bin = 1; bin<=159; bin++) {
2722 if(yprof->GetBinEntries(bin)<5&&
2723 zprof->GetBinEntries(bin)<5) {
2727 // There is a problem in defining inner and outer sectors for
2728 // the outer beams (0 and 6) where both sectors are OROCs. To
2729 // make sure there is no overlap row 94 to 99 are cutted.
2730 if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2733 Int_t row = (bin-1);
2737 Bool_t isOuter = kTRUE;
2738 Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2741 (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2744 sector = TMath::Nint(fBeamSectorInner[id]);
2748 Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2749 vecN[bin-1] =yprof->GetBinEntries(bin);
2750 vecEy[bin-1]=yprof->GetBinError(bin);
2751 vecEz[bin-1]=zprof->GetBinError(bin);
2753 vecDY[bin-1] = yprof->GetBinContent(bin);
2754 vecDZ[bin-1] = zprof->GetBinContent(bin);
2755 if (bin>0&&bin<159){
2757 //truncated mean - skip first and the last pad row
2759 Int_t firstBin=TMath::Max(bin-5,0);
2760 Int_t lastBin =TMath::Min(bin+5,158);
2761 histAbsY->GetXaxis()->SetRangeUser(firstBin,lastBin);
2762 histAbsY->GetYaxis()->SetRangeUser(-2,2);
2763 vecEy[bin-1]=histAbsY->GetRMS(2);
2764 vecDY[bin-1]=histAbsY->GetMean(2);
2765 histAbsY->GetXaxis()->SetRangeUser(firstBin+2,lastBin-2);//use+-2 bins
2766 histAbsY->GetYaxis()->SetRangeUser(vecDY[bin-1]-4*vecEy[bin-1],vecDY[bin-1]+4*vecEy[bin-1]);
2767 if (yprof->GetBinEntries(bin-1)>0) vecEy[bin-1]=histAbsY->GetRMS(2)/TMath::Sqrt(yprof->GetBinEntries(bin-1));
2768 vecDY[bin-1]=histAbsY->GetMean(2);
2771 if(!isOuter) { // inner
2772 vecPhi[bin-1]=lasTanPhiLocIn;
2773 // Calculate local y from residual and database
2774 Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2777 Double_t r = TMath::Sqrt(x*x + y*y);
2779 // Find angle between laser vector and R vector
2780 // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2781 Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2782 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2783 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2784 if(lasTanPhiLocIn<0)
2785 vecPhiR[bin-1]*=-1; // to have the same sign
2787 if(yprof->GetBinEntries(bin)>=10) {
2788 lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2789 TMath::Max(yprof->GetBinError(bin), 0.001));
2791 if(zprof->GetBinEntries(bin)>=10) {
2792 lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2793 TMath::Max(zprof->GetBinError(bin), 0.001));
2796 vecPhi[bin-1]=lasTanPhiLocOut;
2797 // Calculate local y from residual and database
2798 Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2801 Double_t r = TMath::Sqrt(x*x + y*y);
2804 Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2805 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2806 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2807 if(lasTanPhiLocOut<0)
2808 vecPhiR[bin-1]*=-1; // to have the same sign
2810 if(yprof->GetBinEntries(bin)>=10) {
2811 lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2812 TMath::Max(yprof->GetBinError(bin), 0.001));
2814 if(zprof->GetBinEntries(bin)>=10) {
2815 lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2816 TMath::Max(zprof->GetBinError(bin), 0.001));
2823 delete yprof; delete zprof;
2826 // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2827 nClInY = lfabsyInner.GetNpoints();
2828 if(lfabsyInner.GetNpoints()>10) {
2829 lfabsyInner.EvalRobust(0.95);
2832 lfabsyInner.GetParameters(result);
2833 yAbsInOffset = result[0];
2834 yAbsInSlope = result[1];
2835 yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2837 nClInZ = lfabszInner.GetNpoints();
2838 if(lfabszInner.GetNpoints()>10) {
2839 lfabszInner.EvalRobust(0.95);
2842 lfabszInner.GetParameters(result);
2843 zAbsInOffset = result[0];
2844 zAbsInSlope = result[1];
2845 zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2847 nClOutY = lfabsyOuter.GetNpoints();
2848 if(lfabsyOuter.GetNpoints()>10) {
2849 lfabsyOuter.EvalRobust(0.95);
2852 lfabsyOuter.GetParameters(result);
2853 yAbsOutOffset = result[0];
2854 yAbsOutSlope = result[1];
2855 yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2857 nClOutZ = lfabszOuter.GetNpoints();
2858 if(lfabszOuter.GetNpoints()>10) {
2859 lfabszOuter.EvalRobust(0.95);
2862 lfabszOuter.GetParameters(result);
2863 zAbsOutOffset = result[0];
2864 zAbsOutSlope = result[1];
2865 zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2871 Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2872 Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2873 Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2874 Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2876 Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2877 Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2879 Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2880 Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2883 (*pcstream)<<"Mean"<<
2896 "id="<<id<< // track id
2897 "entries="<<entries<< // number of entries
2898 "bz="<<bfield<< // bfield
2899 "LTr.="<<ltrp<< // refernece track
2901 "mnclIn="<<mnclIn<< // mean number of clusters in inner
2902 "mnclOut="<<mnclOut<< // mean number of clusters in outer
2903 "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2904 "mLclIn="<<mLclIn<< // mean number of clusters in inner
2905 "mLclOut="<<mLclOut<< // mean number of clusters in outer
2906 "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2907 "mdEdx="<<mdEdx<< // mean dEdx
2908 "edZfit="<<edZfit<< // entries z fit
2909 "mdZfit="<<mdZfit<< // mean z fit
2910 "rdZfit="<<rdZfit<< // RMS z fit
2913 "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2914 "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2915 "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2916 "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2917 "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2918 "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2919 "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2920 "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2921 "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2922 "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2925 "ePy1vP0="<<ePy1vP0<<
2926 "mPy1vP0="<<mPy1vP0<<
2927 "rPy1vP0="<<rPy1vP0<<
2928 "ePy2vP0="<<ePy2vP0<<
2929 "mPy2vP0="<<mPy2vP0<<
2930 "rPy2vP0="<<rPy2vP0<<
2931 "ePy3vP0="<<ePy3vP0<<
2932 "mPy3vP0="<<mPy3vP0<<
2933 "rPy3vP0="<<rPy3vP0<<
2934 "ePy1vP1="<<ePy1vP1<<
2935 "mPy1vP1="<<mPy1vP1<<
2936 "rPy1vP1="<<rPy1vP1<<
2937 "ePy2vP1="<<ePy2vP1<<
2938 "mPy2vP1="<<mPy2vP1<<
2939 "rPy2vP1="<<rPy2vP1<<
2940 "ePy3vP1="<<ePy3vP1<<
2941 "mPy3vP1="<<mPy3vP1<<
2942 "rPy3vP1="<<rPy3vP1<<
2943 "ePy2vP2In="<<ePy2vP2In<<
2944 "mPy2vP2In="<<mPy2vP2In<<
2945 "rPy2vP2In="<<rPy2vP2In<<
2946 "ePy2vP2Out="<<ePy2vP2Out<<
2947 "mPy2vP2Out="<<mPy2vP2Out<<
2948 "rPy2vP2Out="<<rPy2vP2Out<<
2949 "ePy3vP2IO="<<ePy3vP2IO<<
2950 "mPy3vP2IO="<<mPy3vP2IO<<
2951 "rPy3vP2IO="<<rPy3vP2IO<<
2954 "ePz1vP0="<<ePz1vP0<<
2955 "mPz1vP0="<<mPz1vP0<<
2956 "rPz1vP0="<<rPz1vP0<<
2957 "ePz2vP0="<<ePz2vP0<<
2958 "mPz2vP0="<<mPz2vP0<<
2959 "rPz2vP0="<<rPz2vP0<<
2960 "ePz3vP0="<<ePz3vP0<<
2961 "mPz3vP0="<<mPz3vP0<<
2962 "rPz3vP0="<<rPz3vP0<<
2963 "ePz1vP1="<<ePz1vP1<<
2964 "mPz1vP1="<<mPz1vP1<<
2965 "rPz1vP1="<<rPz1vP1<<
2966 "ePz2vP1="<<ePz2vP1<<
2967 "mPz2vP1="<<mPz2vP1<<
2968 "rPz2vP1="<<rPz2vP1<<
2969 "ePz3vP1="<<ePz3vP1<<
2970 "mPz3vP1="<<mPz3vP1<<
2971 "rPz3vP1="<<rPz3vP1<<
2972 "ePz2vP2In="<<ePz2vP2In<<
2973 "mPz2vP2In="<<mPz2vP2In<<
2974 "rPz2vP2In="<<rPz2vP2In<<
2975 "ePz2vP2Out="<<ePz2vP2Out<<
2976 "mPz2vP2Out="<<mPz2vP2Out<<
2977 "rPz2vP2Out="<<rPz2vP2Out<<
2978 "ePz3vP2IO="<<ePz3vP2IO<<
2979 "mPz3vP2IO="<<mPz3vP2IO<<
2980 "rPz3vP2IO="<<rPz3vP2IO<<
2984 "lx0="<<lxyz[0]<< // reference x
2985 "lx1="<<lxyz[1]<< // reference y
2986 "lx2="<<lxyz[2]<< // refernece z
2987 "lpx0="<<lpxyz[0]<< // reference x
2988 "lpx1="<<lpxyz[1]<< // reference y
2989 "lpx2="<<lpxyz[2]<< // refernece z
2993 "mphi="<<meanphi<< //
2994 "rmsphi="<<rmsphi<< //
2998 "mphiP="<<meanphiP<< //
2999 "rmsphiP="<<rmsphiP<< //
3008 "gp31="<<gp31<< //gaus mean - tgl
3009 "gp32="<<gp32<< //gaus rms -tgl
3010 "meanp3="<<meanp3<< //mean - tgl
3011 "rmsp3="<<rmsp3<< //rms -tgl
3012 "gp41="<<gp41<< //gaus mean - P4
3013 "gp42="<<gp42<< //gaus rms - P4
3014 "meanp4="<<meanp4<< //mean - P4
3015 "rmsp4="<<rmsp4<< //rms - P4
3016 // Parameters from abs res analysis
3017 "SecIn="<<secInner<< // inner sector
3018 "SecOut="<<secOuter<< // outer sector
3019 "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
3020 "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
3021 "ibr.="<<&vecIBR<< // radial filed integral
3022 "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
3023 "ibr.="<<&vecIBR<< // radial filed integral
3024 "ibz.="<<&vecIBZ<< // z filed integral
3026 "iblx.="<<&vecIBLX<< // local bx integral
3027 "ibly.="<<&vecIBLY<< // local by integral
3028 "ibgx.="<<&vecIBGX<< // global bx integral
3029 "ibgy.="<<&vecIBGY<< // global by integral
3031 "X.="<<&vecX<< // local x
3032 "Y.="<<&vecY<< // local y
3033 "R.="<<&vecR<< // radius
3034 "dY.="<<&vecDY<< // abs y residuals
3035 "dZ.="<<&vecDZ<< // abs z residuals
3036 "eY.="<<&vecEy<< // error of y residuals
3037 "eZ.="<<&vecEz<< // error of z residuals
3038 "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
3039 "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
3040 "nCl.="<<&vecN<< // number of clusters
3042 "nClInY="<<nClInY<< // Number of clusters for inner y
3043 "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
3044 "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
3045 "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
3046 "nClInZ="<<nClInZ<< // Number of clusters for inner z
3047 "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
3048 "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
3049 "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
3051 "nClOutY="<<nClOutY<< // Number of clusters for outer y
3052 "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
3053 "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
3054 "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
3055 "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
3056 "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
3057 "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
3058 "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
3065 TFile fmean("laserMean.root");
3075 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3079 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3080 TFile * f = pcstream->GetFile();
3082 f->mkdir("dirphiP");
3084 TF1 fp("p1","pol1");
3089 char grnamefull[1000];
3092 Double_t mphiP[100];
3093 Double_t smphi[100];
3094 Double_t smphiP[100];
3105 for (Int_t id=0; id<336; id++){
3107 snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
3108 Int_t entries = chain->Draw("bz",cut,"goff");
3109 if (entries<3) continue;
3110 AliTPCLaserTrack *ltrp = 0;
3111 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3112 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3116 ltrp->GetPxPyPz(lpxyz);
3118 chain->Draw("bz",cut,"goff");
3119 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3120 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3121 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3123 chain->Draw("gphi1",cut,"goff");
3124 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3125 chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3126 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3128 chain->Draw("gphiP1",cut,"goff");
3129 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3130 chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3131 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3133 chain->Draw("gz1",cut,"goff");
3134 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3135 chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3136 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3139 snprintf(grnamefull,1000,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3140 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3147 for (Int_t ientry=0; ientry<entries; ientry++){
3148 if (TMath::Abs(bz[ientry])<0.05){
3149 phiB0 = mphi[ientry];
3150 phiPB0 = mphiP[ientry];
3154 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3157 pphi[0] = fp.GetParameter(0); // offset
3158 pphi[1] = fp.GetParameter(1); // slope
3159 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3160 snprintf(grname,1000,"phi_id%d",id);
3161 grphi->SetName(grname); grphi->SetTitle(grnamefull);
3162 grphi->GetXaxis()->SetTitle("b_{z} (T)");
3163 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3164 grphi->SetMaximum(1.2);
3165 grphi->SetMinimum(-1.2);
3169 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3172 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3175 pphiP[0] = fp.GetParameter(0); // offset
3176 pphiP[1] = fp.GetParameter(1); // slope
3177 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3178 snprintf(grname,1000,"phiP_id%d",id);
3179 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3180 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3181 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3182 grphiP->SetMaximum(pphiP[0]+0.005);
3183 grphiP->SetMinimum(pphiP[0]-0.005);
3185 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3190 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3193 pmZ[0] = fp.GetParameter(0); // offset
3194 pmZ[1] = fp.GetParameter(1); // slope
3195 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3196 snprintf(grname,1000,"mZ_id%d",id);
3197 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3198 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3199 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3201 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3207 for (Int_t ientry=0; ientry<entries; ientry++){
3208 (*pcstream)<<"Mean"<<
3211 "entries="<<entries<<
3213 "lx0="<<lxyz[0]<< // reference x
3214 "lx1="<<lxyz[1]<< // reference y
3215 "lx2="<<lxyz[2]<< // refernece z
3216 "lpx0="<<lpxyz[0]<< // reference x
3217 "lpx1="<<lpxyz[1]<< // reference y
3218 "lpx2="<<lpxyz[2]<< // refernece z
3220 "phiB0="<<phiB0<< // position shift at 0 field
3221 "phiPB0="<<phiPB0<< // angular shift at 0 field
3222 "zB0="<<zB0<< // z shift for 0 field
3224 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3225 "pphi0="<<pphi[0]<< // offset
3226 "pphi1="<<pphi[1]<< // slope
3227 "pphi2="<<pphi[2]<< // norm chi2
3229 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3230 "pphiP0="<<pphiP[0]<< // offset
3231 "pphiP1="<<pphiP[1]<< // slope
3232 "pphiP2="<<pphiP[2]<< // norm chi2
3234 "gz1="<<mZ[ientry]<<
3235 "pmZ0="<<pmZ[0]<< // offset
3236 "pmZ1="<<pmZ[1]<< // slope
3237 "pmZ2="<<pmZ[2]<< // norm chi2
3247 void AliTPCcalibLaser::Analyze(){
3254 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3256 TIterator* iter = li->MakeIterator();
3257 AliTPCcalibLaser* cal = 0;
3258 static Int_t counter0=0;
3259 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3260 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3261 // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3264 printf("Marging number %d\n", counter0);
3267 MergeFitHistos(cal);
3272 // TProfile *hp=0x0;
3273 //TProfile *hpm=0x0;
3275 for (Int_t id=0; id<336; id++){
3276 // merge fDeltaZ histograms
3277 hm = (TH1F*)cal->fDeltaZ.At(id);
3278 h = (TH1F*)fDeltaZ.At(id);
3279 if (!h &&hm &&hm->GetEntries()>0) {
3280 h=(TH1F*)hm->Clone();
3282 fDeltaZ.AddAt(h,id);
3284 if (h && hm) h->Add(hm);
3286 // merge fP3 histograms
3287 hm = (TH1F*)cal->fDeltaP3.At(id);
3288 h = (TH1F*)fDeltaP3.At(id);
3289 if (!h&&hm &&hm->GetEntries()>0) {
3290 h=(TH1F*)hm->Clone();
3292 fDeltaP3.AddAt(h,id);
3294 if (h && hm) h->Add(hm);
3295 // merge fP4 histograms
3296 hm = (TH1F*)cal->fDeltaP4.At(id);
3297 h = (TH1F*)fDeltaP4.At(id);
3298 if (!h &&hm &&hm->GetEntries()>0) {
3299 h=(TH1F*)hm->Clone();
3301 fDeltaP4.AddAt(h,id);
3303 if (h&&hm) h->Add(hm);
3306 // merge fDeltaPhi histograms
3307 hm = (TH1F*)cal->fDeltaPhi.At(id);
3308 h = (TH1F*)fDeltaPhi.At(id);
3309 if (!h &&hm &&hm->GetEntries()>0) {
3310 h= (TH1F*)hm->Clone();
3312 fDeltaPhi.AddAt(h,id);
3314 if (h&&hm) h->Add(hm);
3315 // merge fDeltaPhiP histograms
3316 hm = (TH1F*)cal->fDeltaPhiP.At(id);
3317 h = (TH1F*)fDeltaPhiP.At(id);
3318 if (!h&&hm &&hm->GetEntries()>0) {
3319 h=(TH1F*)hm->Clone();
3321 fDeltaPhiP.AddAt(h,id);
3323 if (h&&hm) h->Add(hm);
3324 // merge fSignals histograms
3325 hm = (TH1F*)cal->fSignals.At(id);
3326 h = (TH1F*)fSignals.At(id);
3327 if (!h&&hm &&hm->GetEntries()>0) {
3328 h=(TH1F*)hm->Clone();
3330 fSignals.AddAt(h,id);
3332 if (h&&hm) h->Add(hm);
3335 // merge ProfileY histograms -0
3336 h2m = (TH2F*)cal->fDeltaYres.At(id);
3337 h2 = (TH2F*)fDeltaYres.At(id);
3338 if (h2m&&h2) h2->Add(h2m);
3340 h2m = (TH2F*)cal->fDeltaZres.At(id);
3341 h2 = (TH2F*)fDeltaZres.At(id);
3342 if (h2m&&h2) h2->Add(h2m);
3343 // merge ProfileY histograms - 2
3344 h2m = (TH2F*)cal->fDeltaYres2.At(id);
3345 h2 = (TH2F*)fDeltaYres2.At(id);
3346 if (h2m&&h2) h2->Add(h2m);
3348 h2m = (TH2F*)cal->fDeltaZres2.At(id);
3349 h2 = (TH2F*)fDeltaZres2.At(id);
3350 if (h2m&&h2) h2->Add(h2m);
3352 // merge ProfileY histograms - abs
3353 h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3354 h2 = (TH2F*)fDeltaYresAbs.At(id);
3355 if (h2m&&h2) h2->Add(h2m);
3356 if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaYresAbs.AddAt(h2,id);}
3357 h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3358 h2 = (TH2F*)fDeltaZresAbs.At(id);
3359 if (h2m&&h2) h2->Add(h2m);
3360 if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaZresAbs.AddAt(h2,id);}
3361 // merge ProfileY histograms - 3
3362 //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3363 //h2 = (TH2F*)fDeltaYres3.At(id);
3364 //if (h2m) h2->Add(h2m);
3366 //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3367 //h2 = (TH2F*)fDeltaZres3.At(id);
3368 //if (h2m) h->Add(h2m);
3376 void AliTPCcalibLaser::MakeFitHistos(){
3378 // Make a fit histograms
3380 // Number of clusters
3382 //TH2F *fHisNclIn; //->Number of clusters inner
3383 //TH2F *fHisNclOut; //->Number of clusters outer
3384 //TH2F *fHisNclIO; //->Number of cluster inner outer
3385 // TH2F *fHisdEdx; //->dEdx histo
3386 fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3387 fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3388 fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3390 fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3391 fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3392 fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3394 fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3395 fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3400 // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3401 // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3402 // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3403 // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3404 // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3405 fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3406 fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3407 fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3408 fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3409 fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3410 // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3411 // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3412 // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3413 // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3414 // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3415 fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3416 fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3417 fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3418 fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3419 fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3424 // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3425 // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3426 // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3427 // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3428 // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3429 // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3430 // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3431 // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3432 // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3433 fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3434 fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3435 fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3436 fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3437 fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3438 fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3439 fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3440 fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3441 fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3444 // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3445 // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3446 // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3447 // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3448 // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3449 // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3450 // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3451 // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3452 // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3453 fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3454 fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3455 fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3456 fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3457 fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3458 fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3459 fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3460 fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3461 fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3463 fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3464 fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3466 fHisNclIn->SetDirectory(0); //->Number of clusters inner
3467 fHisNclOut->SetDirectory(0); //->Number of clusters outer
3468 fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3469 fHisLclIn->SetDirectory(0); //->Level arm inner
3470 fHisLclOut->SetDirectory(0); //->Level arm outer
3471 fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3472 fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3473 fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3476 fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3477 fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3478 fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3479 fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3480 fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3481 fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3482 fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3483 fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3484 fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3485 fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3488 fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3489 fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3490 fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3491 fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3492 fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3493 fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3494 fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3495 fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3496 fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3499 fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3500 fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3501 fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3502 fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3503 fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3504 fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3505 fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3506 fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3507 fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3509 fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3510 fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3517 for (Int_t id=0; id<336;id++){
3518 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3519 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3520 //TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3522 TH2F *profz2 = 0;//(TH2F*)fDeltaZres2.UncheckedAt(id);
3523 TH2F *profyabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
3524 TH2F *profzabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
3525 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3526 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3528 profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3529 profy->SetDirectory(0);
3530 fDeltaYres.AddAt(profy,id);
3531 profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3532 profy2->SetDirectory(0);
3533 fDeltaYres2.AddAt(profy2,id);
3534 if(!fUseFixedDriftV)
3535 profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3537 profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
3538 profyabs->SetDirectory(0);
3539 fDeltaYresAbs.AddAt(profyabs,id);
3540 //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3541 //profy3->SetDirectory(0);
3542 //fDeltaYres3.AddAt(profy3,id);
3545 profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3546 profz->SetDirectory(0);
3547 fDeltaZres.AddAt(profz,id);
3548 profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3549 profz2->SetDirectory(0);
3550 fDeltaZres2.AddAt(profz2,id);
3551 if(!fUseFixedDriftV)
3552 profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3554 profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
3555 profzabs->SetDirectory(0);
3556 fDeltaZresAbs.AddAt(profzabs,id);
3557 //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3558 //profz3->SetDirectory(0);
3559 //fDeltaZres3.AddAt(profz3,id);
3564 for (Int_t id=0; id<336;id++){
3565 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3566 //TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3570 TH1F * hisdphi = 0;//(TH1F*)fDeltaPhi.At(id);
3571 TH1F * hisdphiP = 0;//(TH1F*)fDeltaPhiP.At(id);
3572 TH1F * hisSignal = 0; //(TH1F*)fSignals.At(id);
3575 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3576 hisdz->SetDirectory(0);
3577 fDeltaZ.AddAt(hisdz,id);
3579 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3580 hisP3->SetDirectory(0);
3581 fDeltaP3.AddAt(hisP3,id);
3583 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3584 hisP4->SetDirectory(0);
3585 fDeltaP4.AddAt(hisP4,id);
3588 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3589 hisdphi->SetDirectory(0);
3590 fDeltaPhi.AddAt(hisdphi,id);
3592 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3593 hisdphiP->SetDirectory(0);
3594 fDeltaPhiP.AddAt(hisdphiP,id);
3595 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3596 hisSignal->SetDirectory(0);
3597 fSignals.AddAt(hisSignal,id);
3601 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3602 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3603 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3604 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3608 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3609 Int_t binsLaser[12]= {336, //id
3621 Double_t xminLaser[12]= {0, //id
3633 Double_t xmaxLaser[12]= {336, //id
3646 TString nameLaser[12]= {"id",
3658 TString titleLaser[12]= {"id",
3670 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3671 for (Int_t iaxis=1; iaxis<12; iaxis++){
3672 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3673 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3677 // Pad SigmaShape Q charge pad row trackID
3678 Int_t binsRow[6]={200, 10000, 20, 30, 159, 336};
3679 Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
3680 Double_t axisMax[6]={ 1, 1000, 1, 30, 159, 336};
3681 TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
3686 fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3695 fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3697 for (Int_t iaxis=0; iaxis<6; iaxis++){
3698 fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
3699 fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
3703 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3705 // Merge content of histograms
3707 // Only first histogram is checked - all other should be the same
3708 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3709 if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
3710 if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
3711 if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
3712 if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
3714 if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
3715 if (!fHisNclIn) MakeFitHistos();
3716 if (fHisNclIn->GetEntries()<1) MakeFitHistos();
3718 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3719 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3720 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3721 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3722 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3723 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3724 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3725 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3728 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3729 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3730 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3731 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3732 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3733 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3734 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3735 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3736 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3737 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3740 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3741 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3742 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3743 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3744 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3745 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3746 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3747 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3748 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3751 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3752 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3753 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3754 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3755 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3756 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3757 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3758 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3759 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3760 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3761 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3774 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3776 // Dump fit information - collect information from the streamers
3779 TChain * chainFit=0;
3780 TChain * chainTrack=0;
3783 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3784 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3785 AliXRDPROOFtoolkit tool;
3786 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3787 chainTrack->Lookup();
3788 chainTrack->SetProof(kTRUE);
3789 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3790 chainDrift->Lookup();
3791 chainDrift->SetProof(kTRUE);
3793 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3795 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3797 chainFit->SetProof(kTRUE);
3798 chain->SetProof(kTRUE);
3799 AliTPCLaserTrack::LoadTracks();
3800 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3806 TCut cutP3("abs(Tr.fP[3])<0.1");
3807 TCut cutP4("abs(Tr.fP[4])<0.5");
3808 TCut cutPx = cutP3+cutP4;
3809 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3810 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3811 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3812 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3814 TCut cutdEdx("sqrt(dEdx)>3");
3815 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3816 TCut cutN("nclO>20&&nclI>20");
3817 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3821 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3822 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3823 TCut cutClX("abs(Cl[].fX)>10");
3824 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3825 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3826 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3827 TCut cutQ("sqrt(Cl[].fMax)>4");
3828 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3832 TH1F * phisAccept = 0;
3834 TProfile * pdEdx = 0;
3842 TProfile * pNclI = 0;
3843 TProfile * pNclO = 0;
3845 TProfile * pchi2YIn =0;
3846 TProfile * pchi2ZIn =0;
3847 TProfile * pchi2YOut =0;
3848 TProfile * pchi2ZOut =0;
3849 TProfile * pchi2YInOut =0;
3850 TProfile * pchi2ZInOut =0;;
3852 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3853 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3854 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3855 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3856 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3857 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3859 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3860 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3863 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3864 pP0 = (TProfile*)gROOT->FindObject("hP0");
3865 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3866 pP1 = (TProfile*)gROOT->FindObject("hP1");
3867 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3868 pP2 = (TProfile*)gROOT->FindObject("hP2");
3869 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3870 pP3 = (TProfile*)gROOT->FindObject("hP3");
3871 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3872 pP4 = (TProfile*)gROOT->FindObject("hP4");
3874 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3875 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3876 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3877 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3880 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3881 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3882 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3883 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3884 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3885 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3886 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3887 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3888 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3889 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3890 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3891 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3893 // second derivatives
3895 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3896 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3897 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3898 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3899 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3900 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3902 phisPy2In->FitSlicesY();
3903 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3904 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3905 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3907 phisPy2Out->FitSlicesY();
3908 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3909 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3910 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3912 phisPy2InOut->FitSlicesY();
3913 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3914 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3915 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3917 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3918 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3919 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3920 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3921 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3922 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3924 phisPz2In->FitSlicesY();
3925 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3926 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3927 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3929 phisPz2Out->FitSlicesY();
3930 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3931 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3932 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3934 phisPz2InOut->FitSlicesY();
3935 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3936 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3937 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3944 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3945 for (Int_t ilaser=0; ilaser<336; ilaser++){
3946 Float_t all=phisAl->GetBinContent(ilaser+1);
3947 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3948 Float_t out=phisOut->GetBinContent(ilaser+1);
3949 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3950 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3951 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3952 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3953 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3954 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3957 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3958 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3960 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3961 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3962 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3963 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3964 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3965 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3967 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3968 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3969 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3971 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3972 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3973 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3975 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3976 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3977 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3979 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3980 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3981 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3983 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3984 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3985 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3987 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3988 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3989 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3991 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3992 (*pcstream)<<"Scan"<<
4006 "chi2YIn="<<chi2YIn<<
4007 "chi2YOut="<<chi2YOut<<
4008 "chi2YInOut="<<chi2YInOut<<
4009 "chi2ZIn="<<chi2ZIn<<
4010 "chi2ZOut="<<chi2ZOut<<
4011 "chi2ZInOut="<<chi2ZInOut<<
4013 "nPy2In="<<entriesPy2In<<
4014 "mPy2In="<<meanPy2In<<
4015 "sPy2In="<<sigmaPy2In<<
4017 "nPy2Out="<<entriesPy2Out<<
4018 "mPy2Out="<<meanPy2Out<<
4019 "sPy2Out="<<sigmaPy2Out<<
4021 "nPy2InOut="<<entriesPy2InOut<<
4022 "mPy2InOut="<<meanPy2InOut<<
4023 "sPy2InOut="<<sigmaPy2InOut<<
4025 "nPz2In="<<entriesPz2In<<
4026 "mPz2In="<<meanPz2In<<
4027 "sPz2In="<<sigmaPz2In<<
4029 "nPz2Out="<<entriesPz2Out<<
4030 "mPz2Out="<<meanPz2Out<<
4031 "sPz2Out="<<sigmaPz2Out<<
4033 "nPz2InOut="<<entriesPz2InOut<<
4034 "mPz2InOut="<<meanPz2InOut<<
4035 "sPz2InOut="<<sigmaPz2InOut<<
4043 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
4044 TVectorD& meanSlope,
4045 TVectorD& sectorArray,
4048 // This method should ideally go in AliTPCLaser
4049 // option == 0 (pads outer - closest to beam)
4050 // option == 1 (pads inner)
4051 // option == 2 (time outer)
4052 // option == 3 (time inner)
4053 Int_t nFailures = 0;
4055 for(Int_t id = 0; id < 336; id++) {
4057 if (!AliTPCLaserTrack::GetTracks())
4058 AliTPCLaserTrack::LoadTracks();
4059 AliTPCLaserTrack *ltrp =
4060 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4062 AliExternalTrackParam trackParam(*ltrp);
4064 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
4065 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
4066 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
4068 Double_t angle = trackParam.GetAlpha();
4070 angle += 2*TMath::Pi();
4071 if(trackParam.GetSnp()>0) // track points to sector "before"
4072 angle -= deltaangle*TMath::DegToRad();
4073 else // track points to sector "after"
4074 angle += deltaangle*TMath::DegToRad();
4076 Bool_t success = trackParam.Rotate(angle);
4079 // cout << "WARNING: Rotate failed for ID: " << id << endl;
4083 angle *= TMath::RadToDeg();
4085 Int_t sector = TMath::Nint((angle-10.0)/20.0);
4090 if(ltrp->GetSide()==1) // C side
4092 if(option==0 || option==2)
4094 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
4097 sectorArray[id] = sector;
4099 const Double_t x0 = 0;
4101 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
4102 Double_t slopez = trackParam.GetTgl();
4103 // One needs a factor sqrt(1+slopey**2) to take into account the
4104 // longer path length
4105 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
4106 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
4108 // Double_t offsetz = trackParam.GetZ();
4109 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4110 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4111 if(option==2 || option==3) {
4112 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4114 meanOffset[id] = offsety; meanSlope[id] = slopey;
4119 AliWarning(Form("Rotate method failed %d times", nFailures));
4124 void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
4127 //input="TPCLaserObjects.root"
4129 // 0. OBJ: TAxis Delta
4130 // 1. OBJ: TAxis bin
4131 // 2. OBJ: TAxis rms shape
4132 // 3. OBJ: TAxis sqrt(Q)
4133 // 4. OBJ: TAxis row
4134 // 5. OBJ: TAxis trackID
4136 const Double_t kSigma=4.;
4138 AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
4139 THnSparse * hisPadInput = laserTPC->fHisLaserPad;
4140 THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
4141 TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
4142 TVectorD meanY(159), sigmaY(159);
4143 TVectorD meanZ(159), sigmaZ(159);
4144 TVectorD meanPad(159), sigmaPad(159);
4145 TVectorD meanTime(159), sigmaTime(159);
4146 TVectorD meanDPad(159), sigmaDPad(159);
4147 TVectorD meanDTime(159), sigmaDTime(159);
4148 TVectorD meandEdx(159), sigmadEdx(159);
4149 TVectorD meanSTime(159), sigmaSTime(159);
4150 TVectorD meanSPad(159), sigmaSPad(159);
4151 TVectorD entries(159);
4153 Int_t indexes[10]={0,1,2,3,4,5,6};
4155 AliTPCLaserTrack::LoadTracks();
4157 for (Int_t id=0; id<336; id++){ // llop over laser beams
4158 printf("id=\t%d\n",id);
4160 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4162 hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
4163 hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
4165 his=hisTimeInput->Projection(3);
4166 Int_t firstBindEdx=his->FindFirstBinAbove(0);
4167 Int_t lastBindEdx=his->FindLastBinAbove(0);
4168 hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4169 hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4172 his=hisTimeInput->Projection(1);
4173 // Int_t firstBinTime=his->FindFirstBinAbove(0);
4174 //Int_t lastBinTime=his->FindLastBinAbove(0);
4175 //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
4179 his=hisTimeInput->Projection(2);
4180 //Int_t firstBinZ=his->FindFirstBinAbove(0);
4181 //Int_t lastBinZ=his->FindLastBinAbove(0);
4182 //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
4185 his=hisPadInput->Projection(2);
4186 // Int_t firstBinY=his->FindFirstBinAbove(0);
4187 //Int_t lastBinY=his->FindLastBinAbove(0);
4188 //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
4193 THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
4194 THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
4197 for (Int_t irow=0; irow<159; irow++){
4199 if ((*(ltrp->GetVecSec()))[irow] <0) continue;
4200 if ((*(ltrp->GetVecLX()))[irow] <80) continue;
4202 hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
4203 hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
4204 //THnSparse *hisPad = hisPad0->Projection(4,indexes);
4205 //THnSparse *hisTime = hisTime0->Projection(4,indexes);
4206 THnSparse *hisPad = hisPad0;
4207 THnSparse *hisTime = hisTime0;
4209 // Get mean value of QA variables
4212 his=hisTime->Projection(3);
4213 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4214 meandEdx[irow] =his->GetMean();
4215 sigmadEdx[irow]=his->GetRMS();
4216 // Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
4217 //Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
4218 // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
4219 //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
4224 his=hisTime->Projection(2);
4225 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
4226 meanSTime[irow] =his->GetMean();
4227 sigmaSTime[irow]=his->GetRMS();
4228 //Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4229 //Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4230 // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
4234 his=hisPad->Projection(2);
4235 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4236 meanSPad[irow] =his->GetMean();
4237 sigmaSPad[irow]=his->GetRMS();
4238 // Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4239 //Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4240 // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
4243 // apply selection on QA variables
4248 his=hisPad->Projection(0);
4249 entries[irow]=his->GetEntries();
4250 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4251 meanY[irow] =his->GetMean();
4252 sigmaY[irow]=his->GetRMS();
4255 his=hisTime->Projection(0);
4256 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4257 meanZ[irow] =his->GetMean();
4258 sigmaZ[irow]=his->GetRMS();
4261 his=hisPad->Projection(1);
4262 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4263 meanPad[irow] =his->GetMean();
4264 meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
4265 sigmaPad[irow]=his->GetRMS();
4268 his=hisTime->Projection(1);
4269 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4270 meanTime[irow] = his->GetMean();
4271 meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
4272 sigmaTime[irow]=his->GetRMS();
4281 (*pcstream)<<"laserClusters"<<
4286 "entries.="<<&entries<<
4287 "my.="<<&meanY<< //mean delta y
4288 "rmsy.="<<&sigmaY<< //rms deltay
4289 "mz.="<<&meanZ<< //mean deltaz
4290 "rmsz.="<<&sigmaZ<< //rms z
4292 "mPad.="<<&meanPad<< // mean pad
4293 "mDPad.="<<&meanDPad<< // mead dpad
4294 "rmsPad.="<<&sigmaPad<< // rms pad
4295 "mTime.="<<&meanTime<<
4296 "mDTime.="<<&meanTime<<
4297 "rmsTime.="<<&sigmaTime<<
4299 "mdEdx.="<<&meandEdx<< //mean dedx
4300 "rmsdEdx.="<<&sigmadEdx<< //rms dedx
4301 "mSPad.="<<&meanSPad<< //mean sigma pad
4302 "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
4303 "mSTime.="<<&meanSTime<< //mean sigma time
4304 "rmsSTime.="<<&sigmaSTime<<
4317 void AliTPCcalibLaser::FitLaserClusters(Int_t run){
4320 //input="TPCLaserObjects.root"
4322 // 1. Select cluster candidates, remove outlyers
4324 // clusters with atypical spread (e.g due track overlaps)
4325 // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
4326 // 2. Fit the tracklets -per sector - in pad and time coordinate frame
4328 // Store info distance of track to pad, time center
4329 // Fit the correction for distance to the center (sin,cos)
4331 const Double_t kEpsilon=0.000001;
4332 const Int_t kMinClusters=20;
4333 const Double_t kEdgeCut=3;
4334 const Double_t kDistCut=1.5; // cut distance to the ideal track
4335 const Double_t kDistCutFit=0.5;
4336 const Double_t kDistCutFitPad=0.25;
4337 const Double_t kDistCutFitTime=0.25;
4338 const Int_t kSmoothRow=5;
4339 TFile f("hisLasers.root"); // Input file
4340 TTree * treeInput=(TTree*)f.Get("laserClusters");
4341 TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
4346 TVectorD *vecTime=0;
4349 TVectorD *meandEdx=0;
4351 TVectorD fitPad(159);
4352 TVectorD fitTime(159);
4353 TVectorD fitPadLocal(159);
4354 TVectorD fitTimeLocal(159);
4355 TVectorD fitDPad(159);
4356 TVectorD fitDTime(159);
4357 TVectorD fitIPad(159);
4358 TVectorD fitITime(159);
4359 Double_t chi2PadIROC=0;
4360 Double_t chi2PadOROC=0;
4362 treeInput->SetBranchAddress("my.",&vecMY);
4363 treeInput->SetBranchAddress("mz.",&vecMZ);
4364 treeInput->SetBranchAddress("mPad.",&vecPad);
4365 treeInput->SetBranchAddress("mTime.",&vecTime);
4366 treeInput->SetBranchAddress("rmsy.",&vecSY);
4367 treeInput->SetBranchAddress("rmsz.",&vecSZ);
4368 treeInput->SetBranchAddress("entries.",&vecN);
4369 treeInput->SetBranchAddress("mdEdx.",&meandEdx);
4371 AliTPCLaserTrack::LoadTracks();
4374 TVectorD fitPadIROC(3), fitPadOROC(3);
4375 TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
4376 TVectorD fitTimeIROC(3), fitTimeOROC(3);
4377 TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
4379 AliTPCROC * roc = AliTPCROC::Instance();
4380 Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
4383 for (Int_t id=0; id<336; id++){
4385 treeInput->GetEntry(id);
4386 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4387 Int_t medianEntries = TMath::Nint(TMath::Median(159,vecN->GetMatrixArray()));
4388 Double_t medianRMSY = TMath::Median(159,vecSY->GetMatrixArray());
4389 Double_t rmsRMSY = TMath::RMS(159,vecSY->GetMatrixArray());
4390 Double_t medianRMSZ = TMath::Median(159,vecSZ->GetMatrixArray());
4391 Double_t rmsRMSZ = TMath::RMS(159,vecSZ->GetMatrixArray());
4392 Double_t mdEdx = TMath::Median(159,meandEdx->GetMatrixArray());
4393 Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
4394 Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
4395 TLinearFitter fitterY(2,"pol1");
4396 TLinearFitter fitterZ(2,"pol1");
4397 TLinearFitter fitterPad(2,"pol1");
4398 TLinearFitter fitterTime(2,"pol1");
4399 TLinearFitter fitterPadSin(2,"hyp1");
4400 TLinearFitter fitterTimeSin(3,"hyp2");
4403 for (UInt_t irow=0; irow<159; irow++){
4404 fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
4405 fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
4406 Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
4410 Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
4411 Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
4412 (*vecPad)[irow]-=npads*0.5;
4414 if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4415 if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4417 if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
4418 if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
4419 if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
4420 if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
4421 if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
4422 if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
4423 if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
4424 Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
4425 if (TMath::Abs(dEdge)<kEdgeCut) continue;
4426 if (irow<roc->GetNRows(0)){
4427 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
4429 if (irow>roc->GetNRows(0)){
4430 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
4436 //fit OROC - get delta pad and delta time
4438 fitterPad.ClearPoints();
4439 fitterTime.ClearPoints();
4440 fitterPadSin.ClearPoints();
4441 fitterTimeSin.ClearPoints();
4442 {for (Int_t irow=2; irow<157; irow++){
4443 if (isOK[irow]<0.5) continue;
4444 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4445 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4446 Double_t y=(*vecPad)[irow];
4447 Double_t z=(*vecTime)[irow];
4448 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4449 fitterPad.AddPoint(&x,y);
4450 fitterTime.AddPoint(&x,z);
4453 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4456 chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4457 for (Int_t irow=2; irow<157; irow++){
4458 if (isOK[irow]<0.5) continue;
4459 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4460 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4461 Double_t y=(*vecPad)[irow];
4462 Double_t z=(*vecTime)[irow];
4463 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4464 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4465 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4466 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4467 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4468 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4469 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4470 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4471 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4472 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4473 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4475 Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
4476 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4477 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4478 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4479 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4482 fitterPadSin.Eval();
4483 fitterTimeSin.Eval();
4484 fitterPadSin.FixParameter(0,0);
4485 fitterTimeSin.FixParameter(0,0);
4486 fitterPadSin.Eval();
4487 fitterTimeSin.Eval();
4489 fitterPad.GetParameters(fitPadOROC);
4490 fitterTime.GetParameters(fitTimeOROC);
4491 fitterPadSin.GetParameters(fitPadOROCSin);
4492 fitterTimeSin.GetParameters(fitTimeOROCSin);
4498 fitterPad.ClearPoints();
4499 fitterTime.ClearPoints();
4500 fitterPadSin.ClearPoints();
4501 fitterTimeSin.ClearPoints();
4502 for (Int_t irow=2; irow<157; irow++){
4503 if (isOK[irow]<0.5) continue;
4504 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4505 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4506 Double_t y=(*vecPad)[irow];
4507 Double_t z=(*vecTime)[irow];
4508 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4509 fitterPad.AddPoint(&x,y);
4510 fitterTime.AddPoint(&x,z);
4513 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4516 chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4517 for (Int_t irow=2; irow<157; irow++){
4518 if (isOK[irow]<0.5) continue;
4519 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4520 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4521 Double_t y=(*vecPad)[irow];
4522 Double_t z=(*vecTime)[irow];
4523 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4524 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4525 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4526 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4527 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4528 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4529 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4530 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4531 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4532 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4533 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4534 if (isOK[irow]>0.5){
4535 Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
4536 TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
4537 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4538 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4539 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4540 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4543 fitterPadSin.Eval();
4544 fitterTimeSin.Eval();
4545 fitterPadSin.FixParameter(0,0);
4546 fitterTimeSin.FixParameter(0,0);
4547 fitterPadSin.Eval();
4548 fitterTimeSin.Eval();
4549 fitterPad.GetParameters(fitPadIROC);
4550 fitterTime.GetParameters(fitTimeIROC);
4551 fitterPadSin.GetParameters(fitPadIROCSin);
4552 fitterTimeSin.GetParameters(fitTimeIROCSin);
4554 for (Int_t irow=0; irow<159; irow++){
4555 if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
4556 if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
4557 if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
4558 if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
4560 for (Int_t irow=kSmoothRow/2; irow<159-kSmoothRow/2; irow++){
4561 fitPadLocal[irow]=0;
4562 fitTimeLocal[irow]=0;
4563 if (isOK[irow]<0.5) continue;
4564 Int_t sector=(irow<Int_t(roc->GetNRows(0)))? sectorInner:sectorOuter;
4565 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
4567 TLinearFitter fitterPadLocal(2,"pol1");
4568 TLinearFitter fitterTimeLocal(2,"pol1");
4569 Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
4570 for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
4571 Int_t jrow=irow+delta;
4573 if (jrow>159) jrow=159;
4574 if (isOK[jrow]<0.5) continue;
4575 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
4576 Double_t y=(*vecPad)[jrow];
4577 Double_t z=(*vecTime)[jrow];
4578 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
4579 fitterPadLocal.AddPoint(&x,y);
4580 fitterTimeLocal.AddPoint(&x,z);
4582 if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
4583 fitterPadLocal.Eval();
4584 fitterTimeLocal.Eval();
4585 fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
4586 fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
4590 (*pcstream)<<"fit"<<
4593 "chi2PadIROC="<<chi2PadIROC<<
4594 "chi2PadOROC="<<chi2PadOROC<<
4598 // mean measured-ideal values
4601 // local coordinate fit
4603 "mTime.="<<vecTime<<
4604 "fitPad.="<<&fitPad<<
4605 "fitTime.="<<&fitTime<<
4606 "fitPadLocal.="<<&fitPadLocal<<
4607 "fitTimeLocal.="<<&fitTimeLocal<<
4608 "fitDPad.="<<&fitDPad<<
4609 "fitDTime.="<<&fitDTime<<
4610 "fitIPad.="<<&fitIPad<<
4611 "fitITime.="<<&fitITime<<
4613 "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
4614 "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
4615 "fitPadOROC.="<<&fitPadOROC<<
4616 "fitPadOROCSin.="<<&fitPadOROCSin<<
4618 "fitTimeIROC.="<<&fitTimeIROC<<
4619 "fitTimeIROCSin.="<<&fitTimeIROCSin<<
4620 "fitTimeOROC.="<<&fitTimeOROC<<
4621 "fitTimeOROCSin.="<<&fitTimeOROCSin<<