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) UpdateFitHistos();
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) UpdateFitHistos();
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);
3604 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3605 Int_t binsLaser[12]= {336, //id
3617 Double_t xminLaser[12]= {0, //id
3629 Double_t xmaxLaser[12]= {336, //id
3642 TString nameLaser[12]= {"id",
3654 TString titleLaser[12]= {"id",
3666 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3667 for (Int_t iaxis=1; iaxis<12; iaxis++){
3668 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3669 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3673 // Pad SigmaShape Q charge pad row trackID
3674 Int_t binsRow[6]={200, 10000, 20, 30, 159, 336};
3675 Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
3676 Double_t axisMax[6]={ 1, 1000, 1, 30, 159, 336};
3677 TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
3682 fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3691 fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3693 for (Int_t iaxis=0; iaxis<6; iaxis++){
3694 fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
3695 fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
3699 void AliTPCcalibLaser::UpdateFitHistos(){
3700 //create the fit histos and set the beam parameters(needs OCDB access)
3702 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3703 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3704 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3705 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3708 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3710 // Merge content of histograms
3712 // Only first histogram is checked - all other should be the same
3713 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3714 if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
3715 if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
3716 if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
3717 if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
3719 if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
3720 if (!fHisNclIn) MakeFitHistos();
3721 if (fHisNclIn->GetEntries()<1) MakeFitHistos();
3723 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3724 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3725 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3726 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3727 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3728 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3729 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3730 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3733 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3734 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3735 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3736 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3737 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3738 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3739 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3740 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3741 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3742 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3745 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3746 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3747 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3748 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3749 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3750 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3751 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3752 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3753 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3756 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3757 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3758 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3759 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3760 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3761 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3762 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3763 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3764 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3765 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3766 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3779 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3781 // Dump fit information - collect information from the streamers
3784 TChain * chainFit=0;
3785 TChain * chainTrack=0;
3788 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3789 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3790 AliXRDPROOFtoolkit tool;
3791 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3792 chainTrack->Lookup();
3793 chainTrack->SetProof(kTRUE);
3794 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3795 chainDrift->Lookup();
3796 chainDrift->SetProof(kTRUE);
3798 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3800 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3802 chainFit->SetProof(kTRUE);
3803 chain->SetProof(kTRUE);
3804 AliTPCLaserTrack::LoadTracks();
3805 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3811 TCut cutP3("abs(Tr.fP[3])<0.1");
3812 TCut cutP4("abs(Tr.fP[4])<0.5");
3813 TCut cutPx = cutP3+cutP4;
3814 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3815 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3816 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3817 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3819 TCut cutdEdx("sqrt(dEdx)>3");
3820 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3821 TCut cutN("nclO>20&&nclI>20");
3822 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3826 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3827 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3828 TCut cutClX("abs(Cl[].fX)>10");
3829 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3830 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3831 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3832 TCut cutQ("sqrt(Cl[].fMax)>4");
3833 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3837 TH1F * phisAccept = 0;
3839 TProfile * pdEdx = 0;
3847 TProfile * pNclI = 0;
3848 TProfile * pNclO = 0;
3850 TProfile * pchi2YIn =0;
3851 TProfile * pchi2ZIn =0;
3852 TProfile * pchi2YOut =0;
3853 TProfile * pchi2ZOut =0;
3854 TProfile * pchi2YInOut =0;
3855 TProfile * pchi2ZInOut =0;;
3857 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3858 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3859 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3860 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3861 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3862 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3864 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3865 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3868 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3869 pP0 = (TProfile*)gROOT->FindObject("hP0");
3870 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3871 pP1 = (TProfile*)gROOT->FindObject("hP1");
3872 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3873 pP2 = (TProfile*)gROOT->FindObject("hP2");
3874 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3875 pP3 = (TProfile*)gROOT->FindObject("hP3");
3876 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3877 pP4 = (TProfile*)gROOT->FindObject("hP4");
3879 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3880 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3881 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3882 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3885 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3886 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3887 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3888 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3889 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3890 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3891 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3892 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3893 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3894 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3895 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3896 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3898 // second derivatives
3900 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3901 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3902 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3903 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3904 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3905 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3907 phisPy2In->FitSlicesY();
3908 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3909 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3910 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3912 phisPy2Out->FitSlicesY();
3913 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3914 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3915 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3917 phisPy2InOut->FitSlicesY();
3918 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3919 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3920 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3922 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3923 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3924 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3925 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3926 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3927 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3929 phisPz2In->FitSlicesY();
3930 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3931 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3932 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3934 phisPz2Out->FitSlicesY();
3935 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3936 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3937 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3939 phisPz2InOut->FitSlicesY();
3940 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3941 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3942 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3949 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3950 for (Int_t ilaser=0; ilaser<336; ilaser++){
3951 Float_t all=phisAl->GetBinContent(ilaser+1);
3952 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3953 Float_t out=phisOut->GetBinContent(ilaser+1);
3954 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3955 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3956 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3957 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3958 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3959 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3962 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3963 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3965 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3966 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3967 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3968 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3969 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3970 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3972 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3973 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3974 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3976 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3977 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3978 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3980 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3981 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3982 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3984 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3985 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3986 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3988 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3989 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3990 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3992 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3993 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3994 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3996 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3997 (*pcstream)<<"Scan"<<
4011 "chi2YIn="<<chi2YIn<<
4012 "chi2YOut="<<chi2YOut<<
4013 "chi2YInOut="<<chi2YInOut<<
4014 "chi2ZIn="<<chi2ZIn<<
4015 "chi2ZOut="<<chi2ZOut<<
4016 "chi2ZInOut="<<chi2ZInOut<<
4018 "nPy2In="<<entriesPy2In<<
4019 "mPy2In="<<meanPy2In<<
4020 "sPy2In="<<sigmaPy2In<<
4022 "nPy2Out="<<entriesPy2Out<<
4023 "mPy2Out="<<meanPy2Out<<
4024 "sPy2Out="<<sigmaPy2Out<<
4026 "nPy2InOut="<<entriesPy2InOut<<
4027 "mPy2InOut="<<meanPy2InOut<<
4028 "sPy2InOut="<<sigmaPy2InOut<<
4030 "nPz2In="<<entriesPz2In<<
4031 "mPz2In="<<meanPz2In<<
4032 "sPz2In="<<sigmaPz2In<<
4034 "nPz2Out="<<entriesPz2Out<<
4035 "mPz2Out="<<meanPz2Out<<
4036 "sPz2Out="<<sigmaPz2Out<<
4038 "nPz2InOut="<<entriesPz2InOut<<
4039 "mPz2InOut="<<meanPz2InOut<<
4040 "sPz2InOut="<<sigmaPz2InOut<<
4048 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
4049 TVectorD& meanSlope,
4050 TVectorD& sectorArray,
4053 // This method should ideally go in AliTPCLaser
4054 // option == 0 (pads outer - closest to beam)
4055 // option == 1 (pads inner)
4056 // option == 2 (time outer)
4057 // option == 3 (time inner)
4058 Int_t nFailures = 0;
4060 for(Int_t id = 0; id < 336; id++) {
4062 if (!AliTPCLaserTrack::GetTracks())
4063 AliTPCLaserTrack::LoadTracks();
4064 AliTPCLaserTrack *ltrp =
4065 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4067 AliExternalTrackParam trackParam(*ltrp);
4069 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
4070 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
4071 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
4073 Double_t angle = trackParam.GetAlpha();
4075 angle += 2*TMath::Pi();
4076 if(trackParam.GetSnp()>0) // track points to sector "before"
4077 angle -= deltaangle*TMath::DegToRad();
4078 else // track points to sector "after"
4079 angle += deltaangle*TMath::DegToRad();
4081 Bool_t success = trackParam.Rotate(angle);
4084 // cout << "WARNING: Rotate failed for ID: " << id << endl;
4088 angle *= TMath::RadToDeg();
4090 Int_t sector = TMath::Nint((angle-10.0)/20.0);
4095 if(ltrp->GetSide()==1) // C side
4097 if(option==0 || option==2)
4099 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
4102 sectorArray[id] = sector;
4104 const Double_t x0 = 0;
4106 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
4107 Double_t slopez = trackParam.GetTgl();
4108 // One needs a factor sqrt(1+slopey**2) to take into account the
4109 // longer path length
4110 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
4111 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
4113 // Double_t offsetz = trackParam.GetZ();
4114 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4115 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4116 if(option==2 || option==3) {
4117 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4119 meanOffset[id] = offsety; meanSlope[id] = slopey;
4124 AliWarning(Form("Rotate method failed %d times", nFailures));
4129 void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
4132 //input="TPCLaserObjects.root"
4134 // 0. OBJ: TAxis Delta
4135 // 1. OBJ: TAxis bin
4136 // 2. OBJ: TAxis rms shape
4137 // 3. OBJ: TAxis sqrt(Q)
4138 // 4. OBJ: TAxis row
4139 // 5. OBJ: TAxis trackID
4141 const Double_t kSigma=4.;
4143 AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
4144 THnSparse * hisPadInput = laserTPC->fHisLaserPad;
4145 THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
4146 TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
4147 TVectorD meanY(159), sigmaY(159);
4148 TVectorD meanZ(159), sigmaZ(159);
4149 TVectorD meanPad(159), sigmaPad(159);
4150 TVectorD meanTime(159), sigmaTime(159);
4151 TVectorD meanDPad(159), sigmaDPad(159);
4152 TVectorD meanDTime(159), sigmaDTime(159);
4153 TVectorD meandEdx(159), sigmadEdx(159);
4154 TVectorD meanSTime(159), sigmaSTime(159);
4155 TVectorD meanSPad(159), sigmaSPad(159);
4156 TVectorD entries(159);
4158 Int_t indexes[10]={0,1,2,3,4,5,6};
4160 AliTPCLaserTrack::LoadTracks();
4162 for (Int_t id=0; id<336; id++){ // llop over laser beams
4163 printf("id=\t%d\n",id);
4165 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4167 hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
4168 hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
4170 his=hisTimeInput->Projection(3);
4171 Int_t firstBindEdx=his->FindFirstBinAbove(0);
4172 Int_t lastBindEdx=his->FindLastBinAbove(0);
4173 hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4174 hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4177 his=hisTimeInput->Projection(1);
4178 // Int_t firstBinTime=his->FindFirstBinAbove(0);
4179 //Int_t lastBinTime=his->FindLastBinAbove(0);
4180 //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
4184 his=hisTimeInput->Projection(2);
4185 //Int_t firstBinZ=his->FindFirstBinAbove(0);
4186 //Int_t lastBinZ=his->FindLastBinAbove(0);
4187 //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
4190 his=hisPadInput->Projection(2);
4191 // Int_t firstBinY=his->FindFirstBinAbove(0);
4192 //Int_t lastBinY=his->FindLastBinAbove(0);
4193 //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
4198 THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
4199 THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
4202 for (Int_t irow=0; irow<159; irow++){
4204 if ((*(ltrp->GetVecSec()))[irow] <0) continue;
4205 if ((*(ltrp->GetVecLX()))[irow] <80) continue;
4207 hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
4208 hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
4209 //THnSparse *hisPad = hisPad0->Projection(4,indexes);
4210 //THnSparse *hisTime = hisTime0->Projection(4,indexes);
4211 THnSparse *hisPad = hisPad0;
4212 THnSparse *hisTime = hisTime0;
4214 // Get mean value of QA variables
4217 his=hisTime->Projection(3);
4218 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4219 meandEdx[irow] =his->GetMean();
4220 sigmadEdx[irow]=his->GetRMS();
4221 // Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
4222 //Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
4223 // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
4224 //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
4229 his=hisTime->Projection(2);
4230 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
4231 meanSTime[irow] =his->GetMean();
4232 sigmaSTime[irow]=his->GetRMS();
4233 //Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4234 //Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4235 // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
4239 his=hisPad->Projection(2);
4240 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4241 meanSPad[irow] =his->GetMean();
4242 sigmaSPad[irow]=his->GetRMS();
4243 // Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4244 //Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4245 // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
4248 // apply selection on QA variables
4253 his=hisPad->Projection(0);
4254 entries[irow]=his->GetEntries();
4255 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4256 meanY[irow] =his->GetMean();
4257 sigmaY[irow]=his->GetRMS();
4260 his=hisTime->Projection(0);
4261 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4262 meanZ[irow] =his->GetMean();
4263 sigmaZ[irow]=his->GetRMS();
4266 his=hisPad->Projection(1);
4267 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4268 meanPad[irow] =his->GetMean();
4269 meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
4270 sigmaPad[irow]=his->GetRMS();
4273 his=hisTime->Projection(1);
4274 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4275 meanTime[irow] = his->GetMean();
4276 meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
4277 sigmaTime[irow]=his->GetRMS();
4286 (*pcstream)<<"laserClusters"<<
4291 "entries.="<<&entries<<
4292 "my.="<<&meanY<< //mean delta y
4293 "rmsy.="<<&sigmaY<< //rms deltay
4294 "mz.="<<&meanZ<< //mean deltaz
4295 "rmsz.="<<&sigmaZ<< //rms z
4297 "mPad.="<<&meanPad<< // mean pad
4298 "mDPad.="<<&meanDPad<< // mead dpad
4299 "rmsPad.="<<&sigmaPad<< // rms pad
4300 "mTime.="<<&meanTime<<
4301 "mDTime.="<<&meanTime<<
4302 "rmsTime.="<<&sigmaTime<<
4304 "mdEdx.="<<&meandEdx<< //mean dedx
4305 "rmsdEdx.="<<&sigmadEdx<< //rms dedx
4306 "mSPad.="<<&meanSPad<< //mean sigma pad
4307 "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
4308 "mSTime.="<<&meanSTime<< //mean sigma time
4309 "rmsSTime.="<<&sigmaSTime<<
4322 void AliTPCcalibLaser::FitLaserClusters(Int_t run){
4325 //input="TPCLaserObjects.root"
4327 // 1. Select cluster candidates, remove outlyers
4329 // clusters with atypical spread (e.g due track overlaps)
4330 // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
4331 // 2. Fit the tracklets -per sector - in pad and time coordinate frame
4333 // Store info distance of track to pad, time center
4334 // Fit the correction for distance to the center (sin,cos)
4336 const Double_t kEpsilon=0.000001;
4337 const Int_t kMinClusters=20;
4338 const Double_t kEdgeCut=3;
4339 const Double_t kDistCut=1.5; // cut distance to the ideal track
4340 const Double_t kDistCutFit=0.5;
4341 const Double_t kDistCutFitPad=0.25;
4342 const Double_t kDistCutFitTime=0.25;
4343 const Int_t kSmoothRow=5;
4344 TFile f("hisLasers.root"); // Input file
4345 TTree * treeInput=(TTree*)f.Get("laserClusters");
4346 TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
4351 TVectorD *vecTime=0;
4354 TVectorD *meandEdx=0;
4356 TVectorD fitPad(159);
4357 TVectorD fitTime(159);
4358 TVectorD fitPadLocal(159);
4359 TVectorD fitTimeLocal(159);
4360 TVectorD fitDPad(159);
4361 TVectorD fitDTime(159);
4362 TVectorD fitIPad(159);
4363 TVectorD fitITime(159);
4364 Double_t chi2PadIROC=0;
4365 Double_t chi2PadOROC=0;
4367 treeInput->SetBranchAddress("my.",&vecMY);
4368 treeInput->SetBranchAddress("mz.",&vecMZ);
4369 treeInput->SetBranchAddress("mPad.",&vecPad);
4370 treeInput->SetBranchAddress("mTime.",&vecTime);
4371 treeInput->SetBranchAddress("rmsy.",&vecSY);
4372 treeInput->SetBranchAddress("rmsz.",&vecSZ);
4373 treeInput->SetBranchAddress("entries.",&vecN);
4374 treeInput->SetBranchAddress("mdEdx.",&meandEdx);
4376 AliTPCLaserTrack::LoadTracks();
4379 TVectorD fitPadIROC(3), fitPadOROC(3);
4380 TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
4381 TVectorD fitTimeIROC(3), fitTimeOROC(3);
4382 TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
4384 AliTPCROC * roc = AliTPCROC::Instance();
4385 Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
4388 for (Int_t id=0; id<336; id++){
4390 treeInput->GetEntry(id);
4391 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4392 Int_t medianEntries = TMath::Nint(TMath::Median(159,vecN->GetMatrixArray()));
4393 Double_t medianRMSY = TMath::Median(159,vecSY->GetMatrixArray());
4394 Double_t rmsRMSY = TMath::RMS(159,vecSY->GetMatrixArray());
4395 Double_t medianRMSZ = TMath::Median(159,vecSZ->GetMatrixArray());
4396 Double_t rmsRMSZ = TMath::RMS(159,vecSZ->GetMatrixArray());
4397 Double_t mdEdx = TMath::Median(159,meandEdx->GetMatrixArray());
4398 Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
4399 Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
4400 TLinearFitter fitterY(2,"pol1");
4401 TLinearFitter fitterZ(2,"pol1");
4402 TLinearFitter fitterPad(2,"pol1");
4403 TLinearFitter fitterTime(2,"pol1");
4404 TLinearFitter fitterPadSin(2,"hyp1");
4405 TLinearFitter fitterTimeSin(3,"hyp2");
4408 for (UInt_t irow=0; irow<159; irow++){
4409 fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
4410 fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
4411 Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
4415 Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
4416 Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
4417 (*vecPad)[irow]-=npads*0.5;
4419 if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4420 if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4422 if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
4423 if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
4424 if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
4425 if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
4426 if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
4427 if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
4428 if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
4429 Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
4430 if (TMath::Abs(dEdge)<kEdgeCut) continue;
4431 if (irow<roc->GetNRows(0)){
4432 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
4434 if (irow>roc->GetNRows(0)){
4435 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
4441 //fit OROC - get delta pad and delta time
4443 fitterPad.ClearPoints();
4444 fitterTime.ClearPoints();
4445 fitterPadSin.ClearPoints();
4446 fitterTimeSin.ClearPoints();
4447 {for (Int_t irow=2; irow<157; irow++){
4448 if (isOK[irow]<0.5) continue;
4449 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4450 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4451 Double_t y=(*vecPad)[irow];
4452 Double_t z=(*vecTime)[irow];
4453 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4454 fitterPad.AddPoint(&x,y);
4455 fitterTime.AddPoint(&x,z);
4458 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4461 chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4462 for (Int_t irow=2; irow<157; irow++){
4463 if (isOK[irow]<0.5) continue;
4464 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4465 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4466 Double_t y=(*vecPad)[irow];
4467 Double_t z=(*vecTime)[irow];
4468 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4469 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4470 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4471 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4472 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4473 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4474 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4475 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4476 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4477 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4478 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4480 Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
4481 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4482 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4483 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4484 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4487 fitterPadSin.Eval();
4488 fitterTimeSin.Eval();
4489 fitterPadSin.FixParameter(0,0);
4490 fitterTimeSin.FixParameter(0,0);
4491 fitterPadSin.Eval();
4492 fitterTimeSin.Eval();
4494 fitterPad.GetParameters(fitPadOROC);
4495 fitterTime.GetParameters(fitTimeOROC);
4496 fitterPadSin.GetParameters(fitPadOROCSin);
4497 fitterTimeSin.GetParameters(fitTimeOROCSin);
4503 fitterPad.ClearPoints();
4504 fitterTime.ClearPoints();
4505 fitterPadSin.ClearPoints();
4506 fitterTimeSin.ClearPoints();
4507 for (Int_t irow=2; irow<157; irow++){
4508 if (isOK[irow]<0.5) continue;
4509 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4510 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4511 Double_t y=(*vecPad)[irow];
4512 Double_t z=(*vecTime)[irow];
4513 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4514 fitterPad.AddPoint(&x,y);
4515 fitterTime.AddPoint(&x,z);
4518 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4521 chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4522 for (Int_t irow=2; irow<157; irow++){
4523 if (isOK[irow]<0.5) continue;
4524 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4525 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4526 Double_t y=(*vecPad)[irow];
4527 Double_t z=(*vecTime)[irow];
4528 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4529 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4530 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4531 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4532 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4533 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4534 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4535 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4536 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4537 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4538 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4539 if (isOK[irow]>0.5){
4540 Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
4541 TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
4542 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4543 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4544 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4545 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4548 fitterPadSin.Eval();
4549 fitterTimeSin.Eval();
4550 fitterPadSin.FixParameter(0,0);
4551 fitterTimeSin.FixParameter(0,0);
4552 fitterPadSin.Eval();
4553 fitterTimeSin.Eval();
4554 fitterPad.GetParameters(fitPadIROC);
4555 fitterTime.GetParameters(fitTimeIROC);
4556 fitterPadSin.GetParameters(fitPadIROCSin);
4557 fitterTimeSin.GetParameters(fitTimeIROCSin);
4559 for (Int_t irow=0; irow<159; irow++){
4560 if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
4561 if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
4562 if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
4563 if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
4565 for (Int_t irow=kSmoothRow/2; irow<159-kSmoothRow/2; irow++){
4566 fitPadLocal[irow]=0;
4567 fitTimeLocal[irow]=0;
4568 if (isOK[irow]<0.5) continue;
4569 Int_t sector=(irow<Int_t(roc->GetNRows(0)))? sectorInner:sectorOuter;
4570 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
4572 TLinearFitter fitterPadLocal(2,"pol1");
4573 TLinearFitter fitterTimeLocal(2,"pol1");
4574 Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
4575 for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
4576 Int_t jrow=irow+delta;
4578 if (jrow>159) jrow=159;
4579 if (isOK[jrow]<0.5) continue;
4580 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
4581 Double_t y=(*vecPad)[jrow];
4582 Double_t z=(*vecTime)[jrow];
4583 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
4584 fitterPadLocal.AddPoint(&x,y);
4585 fitterTimeLocal.AddPoint(&x,z);
4587 if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
4588 fitterPadLocal.Eval();
4589 fitterTimeLocal.Eval();
4590 fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
4591 fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
4595 (*pcstream)<<"fit"<<
4598 "chi2PadIROC="<<chi2PadIROC<<
4599 "chi2PadOROC="<<chi2PadOROC<<
4603 // mean measured-ideal values
4606 // local coordinate fit
4608 "mTime.="<<vecTime<<
4609 "fitPad.="<<&fitPad<<
4610 "fitTime.="<<&fitTime<<
4611 "fitPadLocal.="<<&fitPadLocal<<
4612 "fitTimeLocal.="<<&fitTimeLocal<<
4613 "fitDPad.="<<&fitDPad<<
4614 "fitDTime.="<<&fitDTime<<
4615 "fitIPad.="<<&fitIPad<<
4616 "fitITime.="<<&fitITime<<
4618 "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
4619 "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
4620 "fitPadOROC.="<<&fitPadOROC<<
4621 "fitPadOROCSin.="<<&fitPadOROCSin<<
4623 "fitTimeIROC.="<<&fitTimeIROC<<
4624 "fitTimeIROCSin.="<<&fitTimeIROCSin<<
4625 "fitTimeOROC.="<<&fitTimeOROC<<
4626 "fitTimeOROCSin.="<<&fitTimeOROCSin<<