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"
134 ClassImp(AliTPCcalibLaser)
136 AliTPCcalibLaser::AliTPCcalibLaser():
143 fTracksEsdParam(336),
153 fHisLaser(0), // N dim histogram of laser
154 fHisNclIn(0), //->Number of clusters inner
155 fHisNclOut(0), //->Number of clusters outer
156 fHisNclIO(0), //->Number of cluster inner outer
157 fHisLclIn(0), //->Level arm inner
158 fHisLclOut(0), //->Level arm outer
159 fHisLclIO(0), //->Number of cluster inner outer
160 fHisdEdx(0), //->dEdx histo
161 fHisdZfit(0), //->distance to the mirror after linear fit
164 fHisChi2YIn1(0), //->chi2 y inner - line
165 fHisChi2YOut1(0), //->chi2 y inner - line
166 fHisChi2YIn2(0), //->chi2 y inner - parabola
167 fHisChi2YOut2(0), //->chi2 y inner - parabola
168 fHisChi2YIO1(0), //->chi2 y IO - common
169 fHisChi2ZIn1(0), //->chi2 z inner - line
170 fHisChi2ZOut1(0), //->chi2 z inner - line
171 fHisChi2ZIn2(0), //->chi2 z inner - parabola
172 fHisChi2ZOut2(0), //->chi2 z inner - parabola
173 fHisChi2ZIO1(0), //->chi2 z IO - common
176 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
177 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
178 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
179 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
180 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
181 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
182 fHisPy2vP2In(0), //-> Curv P2inner - parabola
183 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
184 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
187 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
188 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
189 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
190 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
191 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
192 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
193 fHisPz2vP2In(0), //-> Curv P2inner - parabola
194 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
195 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
197 fDeltaYres(336), //->2D histo of residuals
198 fDeltaZres(336), //->2D histo fo residuals
199 fDeltaYres2(336), //->2D histo of residuals
200 fDeltaZres2(336), //->2D histo fo residuals
201 fDeltaYresAbs(336), //->2D histo of residuals
202 fHisYAbsErrors(0), //-> total errors per beam in the abs res y analysis
203 fDeltaZresAbs(336), //->2D histo of residuals
204 fHisZAbsErrors(0), //-> total errors per beam in the abs res z analysis
205 //fDeltaYres3(336), //->2D histo of residuals
206 //fDeltaZres3(336), //->2D histo fo residuals
207 fFitAside(new TVectorD(5)),
208 fFitCside(new TVectorD(5)),
209 fFitACside(new TVectorD(6)),
214 fBeamSectorOuter(336),
215 fBeamSectorInner(336),
216 fBeamOffsetYOuter(336),
217 fBeamSlopeYOuter(336),
218 fBeamOffsetYInner(336),
219 fBeamSlopeYInner(336),
220 fBeamOffsetZOuter(336),
221 fBeamSlopeZOuter(336),
222 fBeamOffsetZInner(336),
223 fBeamSlopeZInner(336),
224 fInverseSlopeZ(kTRUE),
226 fFixedFitAside0(0.0),
227 fFixedFitAside1(1.0),
228 fFixedFitCside0(0.0),
234 fTracksEsdParam.SetOwner(kTRUE);
237 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
244 fTracksEsdParam(336),
248 fDeltaZ(336), // array of histograms of delta z for each track
249 fDeltaP3(336), // array of histograms of delta z for each track
250 fDeltaP4(336), // array of histograms of P3 for each track
251 fDeltaPhi(336), // array of histograms of P4 for each track
252 fDeltaPhiP(336), // array of histograms of delta z for each track
253 fSignals(336), // array of dedx signals
256 fHisLaser(0), // N dim histogram of laser
257 fHisNclIn(0), //->Number of clusters inner
258 fHisNclOut(0), //->Number of clusters outer
259 fHisNclIO(0), //->Number of cluster inner outer
260 fHisLclIn(0), //->Level arm inner
261 fHisLclOut(0), //->Level arm outer
262 fHisLclIO(0), //->Number of cluster inner outer
263 fHisdEdx(0), //->dEdx histo
264 fHisdZfit(0), //->distance to the mirror after linear fit
267 fHisChi2YIn1(0), //->chi2 y inner - line
268 fHisChi2YOut1(0), //->chi2 y inner - line
269 fHisChi2YIn2(0), //->chi2 y inner - parabola
270 fHisChi2YOut2(0), //->chi2 y inner - parabola
271 fHisChi2YIO1(0), //->chi2 y IO - common
272 fHisChi2ZIn1(0), //->chi2 z inner - line
273 fHisChi2ZOut1(0), //->chi2 z inner - line
274 fHisChi2ZIn2(0), //->chi2 z inner - parabola
275 fHisChi2ZOut2(0), //->chi2 z inner - parabola
276 fHisChi2ZIO1(0), //->chi2 z IO - common
279 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
280 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
281 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
282 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
283 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
284 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
285 fHisPy2vP2In(0), //-> Curv P2inner - parabola
286 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
287 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
290 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
291 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
292 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
293 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
294 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
295 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
296 fHisPz2vP2In(0), //-> Curv P2inner - parabola
297 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
298 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
312 fFitAside(new TVectorD(5)), // drift fit - A side
313 fFitCside(new TVectorD(5)), // drift fit - C- side
314 fFitACside(new TVectorD(6)), // drift fit - AC- side
315 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
316 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
317 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
318 fNcuts(0), // number of cuts
319 fBeamSectorOuter(336),
320 fBeamSectorInner(336),
321 fBeamOffsetYOuter(336),
322 fBeamSlopeYOuter(336),
323 fBeamOffsetYInner(336),
324 fBeamSlopeYInner(336),
325 fBeamOffsetZOuter(336),
326 fBeamSlopeZOuter(336),
327 fBeamOffsetZInner(336),
328 fBeamSlopeZInner(336),
329 fInverseSlopeZ(kTRUE),
331 fFixedFitAside0(0.0),
332 fFixedFitAside1(1.0),
333 fFixedFitCside0(0.0),
341 fTracksEsdParam.SetOwner(kTRUE);
344 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
345 AliTPCcalibBase(calibLaser),
351 fTracksEsdParam(336),
353 fFullCalib(calibLaser.fFullCalib),
355 fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
356 fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
357 fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
358 fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
359 fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
360 fSignals(((calibLaser.fSignals))), // array of dedx signals
363 fHisLaser(0), // N dim histogram of laser
364 fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
365 fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
366 fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
367 fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
368 fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
369 fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
370 fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
371 fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
374 fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
375 fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
376 fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
377 fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
378 fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
379 fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
380 fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
381 fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
382 fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
383 fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
386 fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
387 fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
388 fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
389 fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
390 fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
391 fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
392 fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
393 fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
394 fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
397 fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
398 fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
399 fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
400 fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
401 fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
402 fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
403 fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
404 fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
405 fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
408 fDeltaYres(((calibLaser.fDeltaYres))),
409 fDeltaZres(((calibLaser.fDeltaZres))),
410 fDeltaYres2(((calibLaser.fDeltaYres))),
411 fDeltaZres2(((calibLaser.fDeltaZres))),
412 fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
413 fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
414 fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
415 fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
416 // fDeltaYres3(((calibLaser.fDeltaYres))),
417 //fDeltaZres3(((calibLaser.fDeltaZres))),
418 fFitAside(new TVectorD(5)), // drift fit - A side
419 fFitCside(new TVectorD(5)), // drift fit - C- side
420 fFitACside(new TVectorD(6)), // drift fit - C- side
421 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
422 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
423 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
424 fNcuts(0), // number of cuts
425 fBeamSectorOuter(336),
426 fBeamSectorInner(336),
427 fBeamOffsetYOuter(336),
428 fBeamSlopeYOuter(336),
429 fBeamOffsetYInner(336),
430 fBeamSlopeYInner(336),
431 fBeamOffsetZOuter(336),
432 fBeamSlopeZOuter(336),
433 fBeamOffsetZInner(336),
434 fBeamSlopeZInner(336),
435 fInverseSlopeZ(calibLaser.fInverseSlopeZ),
436 fUseFixedDriftV(calibLaser.fUseFixedDriftV),
437 fFixedFitAside0(calibLaser.fFixedFitAside0),
438 fFixedFitAside1(calibLaser.fFixedFitAside1),
439 fFixedFitCside0(calibLaser.fFixedFitCside0),
440 fFixedFitCside1(calibLaser.fFixedFitCside1)
449 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
451 // assgnment operator
453 if (this != &calibLaser) {
454 new (this) AliTPCcalibLaser(calibLaser);
463 AliTPCcalibLaser::~AliTPCcalibLaser() {
468 delete fHisLaser; //->
469 delete fHisNclIn; //->Number of clusters inner
470 delete fHisNclOut; //->Number of clusters outer
471 delete fHisNclIO; //->Number of cluster inner outer
472 delete fHisLclIn; //->Level arm inner
473 delete fHisLclOut; //->Level arm outer
474 delete fHisLclIO; //->Number of cluster inner outer
479 delete fHisChi2YIn1; //->chi2 y inner - line
480 delete fHisChi2YOut1; //->chi2 y inner - line
481 delete fHisChi2YIn2; //->chi2 y inner - parabola
482 delete fHisChi2YOut2; //->chi2 y inner - parabola
483 delete fHisChi2YIO1; //->chi2 y IO - common
484 delete fHisChi2ZIn1; //->chi2 z inner - line
485 delete fHisChi2ZOut1; //->chi2 z inner - line
486 delete fHisChi2ZIn2; //->chi2 z inner - parabola
487 delete fHisChi2ZOut2; //->chi2 z inner - parabola
488 delete fHisChi2ZIO1; //->chi2 z IO - common
491 delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
492 delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
493 delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
494 delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
495 delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
496 delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
497 delete fHisPy2vP2In; //-> Curv P2inner - parabola
498 delete fHisPy2vP2Out; //-> Curv P2outer - parabola
499 delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
502 delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
503 delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
504 delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
505 delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
506 delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
507 delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
508 delete fHisPz2vP2In; //-> Curv P2inner - parabola
509 delete fHisPz2vP2Out; //-> Curv P2outer - parabola
510 delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
515 fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
516 fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
517 fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
518 fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
519 fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
520 fSignals.SetOwner(); //->Array of dedx signals
522 fDeltaZ.Delete(); //-> array of histograms of delta z for each track
523 fDeltaP3.Delete(); //-> array of histograms of P3 for each track
524 fDeltaP4.Delete(); //-> array of histograms of P4 for each track
525 fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
526 fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
527 fSignals.Delete(); //->Array of dedx signals
529 fDeltaYres.SetOwner();
531 delete fHisYAbsErrors;
532 fDeltaZres.SetOwner();
534 delete fHisZAbsErrors;
535 fDeltaYres2.SetOwner();
536 fDeltaYres2.Delete();
537 fDeltaZres2.SetOwner();
538 fDeltaZres2.Delete();
540 fDeltaYresAbs.SetOwner();
541 fDeltaYresAbs.Delete();
542 fDeltaZresAbs.SetOwner();
543 fDeltaZresAbs.Delete();
548 void AliTPCcalibLaser::Process(AliESDEvent * event) {
551 // Loop over tracks and call Process function
553 const Int_t kMinTracks=20;
554 const Int_t kMinClusters=40;
560 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
564 if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
565 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
567 // find CE background if present
569 if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
570 TH1D hisCE("hhisCE","hhisCE",100,-100,100);
571 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
572 AliESDtrack *track=fESD->GetTrack(i);
573 if (!track) continue;
574 hisCE.Fill(track->GetZ());
575 hisCE.Fill(track->GetZ()+2);
576 hisCE.Fill(track->GetZ()-2);
584 fTracksEsdParam.Delete();
585 for (Int_t id=0; id<336;id++) {
587 fClusterCounter[id]=0;
591 Int_t n=fESD->GetNumberOfTracks();
593 for (Int_t i=0;i<n;++i) {
594 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
595 if (!friendTrack) continue;
596 AliESDtrack *track=fESD->GetTrack(i);
597 Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
598 if (binC>336) continue; //remove CE background
599 TObject *calibObject=0;
601 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
602 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
604 if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
606 Int_t id = FindMirror(track,seed);
607 if (id>=0) counter++;
612 if (counter<kMinTracks) return;
616 if (!fFullCalib) return;
617 static Bool_t init=kFALSE;
619 init = kTRUE; // way around for PROOF - to be investigated
623 for (Int_t id=0; id<336; id++){
625 if (!fTracksEsdParam.At(id)) continue;
626 if (fClusterSatur[id]>0.3) continue; // tracks in saturation
628 if ( AcceptLaser(id) && fFitZ[id]<0.5){
636 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
640 // for (Int_t id=0; id<336; id++){
643 if (!fTracksEsdParam.At(id)) return;
644 if (!AcceptLaser(id)) return;
645 Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
648 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
649 if (!hisdz) MakeFitHistos();
650 hisdz = (TH1F*)fDeltaZ.At(id);
651 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
652 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
653 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
654 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
655 TH1F * hisSignal = (TH1F*)fSignals.At(id);
658 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
659 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
660 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
669 param->GetPxPyPz(pxyz);
671 ltrp->GetPxPyPz(lpxyz);
673 Float_t dz = param->GetZ()-ltrp->GetZ();
674 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
675 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
676 if (hisdz) hisdz->Fill(dz);
677 if (hisP3) hisP3->Fill(param->GetParameter()[3]);
678 if (hisP4) hisP4->Fill(param->GetParameter()[4]);
679 if (hisdphi) hisdphi->Fill(dphi);
680 if (hisdphiP) hisdphiP->Fill(dphiP);
681 if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
683 xhis[0] = ltrp->GetId();
684 xhis[1] = ltrp->GetSide();
685 xhis[2] = ltrp->GetRod();
686 xhis[3] = ltrp->GetBundle();
687 xhis[4] = ltrp->GetBeam();
690 xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
691 xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
692 xhis[9] = param->GetParameter()[4];
693 xhis[10]= track->GetTPCNcls();
694 xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
696 fHisLaser->Fill(xhis);
699 void AliTPCcalibLaser::FitDriftV(){
701 // Fit corrections to the drift velocity - linear approximation in the z and global y
702 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
707 z = s* (z0 - vd*(t-t0))
711 vd - nominal drift velocity
712 zs - miscalibrated position
714 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
715 vr - relative change of the drift velocity
720 zs ~ z - s*vr*(z0-s*z)+s*dzt
721 --------------------------------
722 1. Correction function vr constant:
725 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
726 dzs/dl = dz/dl +s*s*vr*dz/dl
730 const Float_t kZCut = 240; // remove the closest laser beam
731 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
732 const Float_t kDistCut = 3; // distance sigma cut
733 const Float_t kDistCutAbs = 0.25;
734 const Float_t kMinClusters = 60; // minimal amount of the clusters
735 const Float_t kMinSignal = 16; // minimal mean height of the signal
736 const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
737 static TLinearFitter fdriftA(3,"hyp2");
738 static TLinearFitter fdriftC(3,"hyp2");
739 static TLinearFitter fdriftAC(4,"hyp3");
740 TVectorD fitA(3),fitC(3),fitAC(4);
742 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
743 AliTPCParam * tpcparam = calib->GetParameters();
745 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
755 for (Int_t iter=0; iter<3; iter++){
756 fdriftA.ClearPoints();
757 fdriftC.ClearPoints();
758 fdriftAC.ClearPoints();
760 for (Int_t id=0; id<336; id++){
761 if (!fTracksEsdParam.At(id)) continue;
762 if (!AcceptLaser(id)) continue;
763 if ( fClusterSatur[id]>kSaturCut) continue;
764 if ( fClusterCounter[id]<kMinClusters) continue;
765 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
766 if (track->GetTPCsignal()<kMinSignal) continue;
767 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
768 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
775 param->GetPxPyPz(pxyz);
777 ltrp->GetPxPyPz(lpxyz);
778 if (TMath::Abs(lxyz[2])>kZCut) continue;
779 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
780 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
781 if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
782 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
785 Double_t zlength = tpcparam->GetZLength(0);
786 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
787 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
788 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
789 if (ltrp->GetSide()==0){
790 fdriftA.AddPoint(xxx,mdrift,1);
792 fdriftC.AddPoint(xxx,mdrift,1);
794 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
795 fdriftAC.AddPoint(xxx2,mdrift,1);
798 if (fdriftA.GetNpoints()>10){
801 if (iter==0) fdriftA.EvalRobust(0.7);
802 else fdriftA.EvalRobust(0.8);
803 fdriftA.GetParameters(fitA);
804 npointsA= fdriftA.GetNpoints();
805 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
806 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
807 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
808 (*fFitAside)[0] = fitA[0];
809 (*fFitAside)[1] = fitA[1];
810 (*fFitAside)[2] = fitA[2];
811 (*fFitAside)[3] = fdriftA.GetNpoints();
812 (*fFitAside)[4] = chi2A;
815 if (fdriftC.GetNpoints()>10){
817 if (iter==0) fdriftC.EvalRobust(0.7);
818 else fdriftC.EvalRobust(0.8);
820 fdriftC.GetParameters(fitC);
821 npointsC= fdriftC.GetNpoints();
822 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
823 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
824 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
825 (*fFitCside)[0] = fitC[0];
826 (*fFitCside)[1] = fitC[1];
827 (*fFitCside)[2] = fitC[2];
828 (*fFitCside)[3] = fdriftC.GetNpoints();
829 (*fFitCside)[4] = chi2C;
833 if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
835 if (iter==0) fdriftAC.EvalRobust(0.7);
836 else fdriftAC.EvalRobust(0.8);
838 fdriftAC.GetParameters(fitAC);
839 npointsAC= fdriftAC.GetNpoints();
840 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
841 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
844 for (Int_t id=0; id<336; id++){
845 if (!fTracksEsdParam.At(id)) continue;
847 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
848 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
854 param->GetPxPyPz(pxyz);
856 ltrp->GetPxPyPz(lpxyz);
857 Double_t zlength = tpcparam->GetZLength(0);
858 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
859 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
862 if (ltrp->GetSide()==0){
863 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
865 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
868 fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
873 TTreeSRedirector *cstream = GetDebugStreamer();
874 TTimeStamp tstamp(fTime);
875 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
876 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
877 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
878 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
879 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
880 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
881 TVectorD vecGoofie(20);
882 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
884 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
885 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
886 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
890 (*cstream)<<"driftv"<<
891 "run="<<fRun<< // run number
892 "event="<<fEvent<< // event number
893 "time="<<fTime<< // time stamp of event
894 "trigger="<<fTrigger<< // trigger
895 "mag="<<fMagF<< // magnetic field
896 // Environment values
897 "press0="<<valuePressure0<<
898 "press1="<<valuePressure1<<
899 "pt0="<<ptrelative0<<
900 "pt1="<<ptrelative1<<
903 "vecGoofie.="<<&vecGoofie<<
907 "driftA.="<<fFitAside<<
908 "driftC.="<<fFitCside<<
909 "driftAC.="<<fFitACside<<
921 Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
923 // Fit corrections to the drift velocity - linear approximation in the z and global y
924 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
926 // Source of outlyers :
927 // 0. Track in the saturation - postpeak
928 // 1. gating grid close the part of the signal for first bundle
930 // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
931 // 1. Robust fit is used in the itteration number 0
932 // only fraction of laser uted
933 // 2. Only the tracks close to the fit used in the second itteration
937 z = s* (z0 - vd*(t-t0))
941 vd - nominal drift velocity
942 zs - miscalibrated position
944 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
945 vr - relative change of the drift velocity
950 zs ~ z - s*vr*(z0-s*z)+s*dzt
951 --------------------------------
952 1. Correction function vr constant:
955 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
956 dzs/dl = dz/dl +s*s*vr*dz/dl
959 const Int_t knLaser = 336; //n laser tracks
960 const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
962 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
963 const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
964 const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
965 const Float_t kMinClusters = 40.; // minimal amount of the clusters
966 const Float_t kMinSignal = 2.5; // minimal mean height of the signal
967 const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
969 static TLinearFitter fdriftA(3,"hyp2");
970 static TLinearFitter fdriftC(3,"hyp2");
971 static TLinearFitter fdriftAC(4,"hyp3");
972 TVectorD fitA(3),fitC(3),fitAC(4);
974 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
975 AliTPCParam * tpcparam = calib->GetParameters();
979 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
980 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
981 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
982 for (Int_t i=0;i<5; i++){
994 Int_t nbA[4]={0,0,0,0};
995 Int_t nbC[4]={0,0,0,0};
996 TVectorD vecZM(336); // measured z potion of laser
997 TVectorD vecA(336); // accepted laser
998 TVectorD vecZF(336); // fitted position
999 TVectorD vecDz(336); // deltaZ
1000 TVectorD vecZS(336); // surveyed position of laser
1001 // additional variable to cut
1002 TVectorD vecdEdx(336); // dEdx
1003 TVectorD vecSy(336); // shape y
1004 TVectorD vecSz(336); // shape z
1007 for (Int_t id=0; id<336; id++){
1009 AliTPCLaserTrack *ltrp =
1010 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1011 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1012 AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1013 vecZM(id)= (param==0) ? 0:param->GetZ();
1014 vecZS(id)= ltrp->GetZ();
1017 vecdEdx(id)=(seed)?seed->GetdEdx():0;
1018 vecSy(id) =(seed)?seed->CookShape(1):0;
1019 vecSz(id) =(seed)?seed->CookShape(2):0;
1022 if (!fTracksEsdParam.At(id)) reject|=1;
1023 if (!param) continue;
1024 if (!AcceptLaser(id)) reject|=2;
1025 if ( fClusterSatur[id]>kSaturCut) reject|=4;
1026 if ( fClusterCounter[id]<kMinClusters) reject|=8;
1028 if (reject>0) continue;
1029 if (ltrp->GetSide()==0){
1031 nbA[ltrp->GetBundle()]++;
1033 if (ltrp->GetSide()==1){
1035 nbC[ltrp->GetBundle()]++;
1039 // reject "bad events"
1044 if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1046 for (Int_t i=0;i<4;i++){
1047 //count accepted for all layers
1048 if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1049 if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1051 if (naA<3 &&naC<3) isOK=kFALSE;
1052 if (isOK==kFALSE) return kFALSE;
1057 for (Int_t iter=0; iter<2; iter++){
1058 fdriftA.ClearPoints();
1059 fdriftC.ClearPoints();
1060 fdriftAC.ClearPoints();
1061 npointsA=0; npointsC=0; npointsAC=0;
1063 for (Int_t id=0; id<336; id++){
1064 if (!fTracksEsdParam.At(id)) continue;
1065 if (!AcceptLaser(id)) continue;
1066 if ( fClusterSatur[id]>kSaturCut) continue;
1067 if ( fClusterCounter[id]<kMinClusters) continue;
1068 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1069 if (track->GetTPCsignal()<kMinSignal) continue;
1070 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1071 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1077 param->GetPxPyPz(pxyz);
1079 ltrp->GetPxPyPz(lpxyz);
1080 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1081 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1082 if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1083 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1085 // // drift distance
1086 // Double_t zlength = tpcparam->GetZLength(0);
1087 // Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1088 // Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1090 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1091 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1092 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1095 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1096 if (iter==0 &<rp->GetBundle()==0) continue;
1097 // skip bundle 0 - can be wrong in the 0 iteration
1098 if (ltrp->GetSide()==0){
1099 fdriftA.AddPoint(xxx,mdrift,1);
1101 fdriftC.AddPoint(xxx,mdrift,1);
1103 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1104 fdriftAC.AddPoint(xxx2,mdrift,1);
1107 if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1110 npointsA= fdriftA.GetNpoints();
1111 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1112 fdriftA.EvalRobust(kFraction[iter]);
1113 fdriftA.GetParameters(fitA);
1114 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1115 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1116 (*fFitAside)[0] = fitA[0];
1117 (*fFitAside)[1] = fitA[1];
1118 (*fFitAside)[2] = fitA[2];
1119 (*fFitAside)[3] = fdriftA.GetNpoints();
1120 (*fFitAside)[4] = chi2A;
1123 if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1125 npointsC= fdriftC.GetNpoints();
1126 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1127 fdriftC.EvalRobust(kFraction[iter]);
1128 fdriftC.GetParameters(fitC);
1129 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1130 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1131 (*fFitCside)[0] = fitC[0];
1132 (*fFitCside)[1] = fitC[1];
1133 (*fFitCside)[2] = fitC[2];
1134 (*fFitCside)[3] = fdriftC.GetNpoints();
1135 (*fFitCside)[4] = chi2C;
1139 if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1141 npointsAC= fdriftAC.GetNpoints();
1142 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1143 fdriftAC.EvalRobust(kFraction[iter]);
1144 fdriftAC.GetParameters(fitAC);
1145 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1146 (*fFitACside)[0] = fitAC[0];
1147 (*fFitACside)[1] = fitAC[1];
1148 (*fFitACside)[2] = fitAC[2];
1149 (*fFitACside)[3] = fdriftAC.GetNpoints();
1150 (*fFitACside)[4] = chi2AC;
1153 for (Int_t id=0; id<336; id++){
1154 if (!fTracksEsdParam.At(id)) continue;
1156 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1157 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1163 param->GetPxPyPz(pxyz);
1165 ltrp->GetPxPyPz(lpxyz);
1166 //Double_t zlength = tpcparam->GetZLength(0);
1167 //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1168 //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1169 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1170 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1171 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1175 if (ltrp->GetSide()==0){
1176 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1178 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1181 //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1183 fFitZ[id]=mdrift-fz;
1185 vecDz[id]=mdrift-fz;
1187 if (fStreamLevel>0){
1188 TTreeSRedirector *cstream = GetDebugStreamer();
1189 TTimeStamp tstamp(fTime);
1190 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1191 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1192 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1193 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1194 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1195 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1196 TVectorD vecGoofie(20);
1197 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1199 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1200 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1201 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1205 (*cstream)<<"driftvN"<<
1206 "run="<<fRun<< // run number
1207 "event="<<fEvent<< // event number
1208 "time="<<fTime<< // time stamp of event
1209 "trigger="<<fTrigger<< // trigger
1210 "mag="<<fMagF<< // magnetic field
1211 // Environment values
1212 "press0="<<valuePressure0<<
1213 "press1="<<valuePressure1<<
1214 "pt0="<<ptrelative0<<
1215 "pt1="<<ptrelative1<<
1218 "vecGoofie.="<<&vecGoofie<<
1221 "vecZM.="<<&vecZM<< // measured z position
1222 "vecZS.="<<&vecZS<< // surveyed z position
1223 "vecZF.="<<&vecZF<< // fitted z position
1224 "vecDz.="<<&vecDz<< // fitted z position
1225 "vecA.="<<&vecA<< // accept laser flag
1226 "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
1227 "vecSy.="<<&vecSy<< // shape y - to cut on
1228 "vecSz.="<<&vecSz<< // shape z - to cut on
1231 "driftA.="<<fFitAside<<
1232 "driftC.="<<fFitCside<<
1233 "driftAC.="<<fFitACside<<
1243 variables to check in debug mode:
1245 chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1246 chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1247 chainDriftN->SetAlias("driftF","vecZF.fElements");
1248 chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
1249 TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1250 TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1259 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1261 // get distance between mirror and laser track
1272 dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1275 // apply drift correction if already exist
1278 if (ltrp->GetSide()==0){
1279 if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1281 if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1283 if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1284 dist+=TMath::Abs(dz);
1286 // phi dist - divergence on 50 cm
1288 dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1293 Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
1298 TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1299 TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1300 TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1301 TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1303 TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1305 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1306 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1307 //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1312 if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1313 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
1314 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
1315 if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
1322 Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1324 // Find corresponding mirror
1325 // add the corresponding tracks
1327 if (!track->GetOuterParam()) return -1;
1329 Float_t kRadius0 = 252;
1330 Float_t kRadius = 254.2;
1332 Float_t counterSatur=0;
1335 for (Int_t irow=158;irow>-1;--irow) {
1336 AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1338 Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1339 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1340 if (pedgeY<3) continue;
1341 if (pedgeX<3) continue;
1343 if (c->GetDetector()%36<18) csideA++;
1344 if (c->GetDetector()%36>=18) csideC++;
1345 if (c->GetMax()>900) counterSatur++;
1347 counterSatur/=(countercl+1);
1351 if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
1354 if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1357 AliExternalTrackParam param(*(track->GetOuterParam()));
1358 AliTracker::PropagateTrackTo(¶m,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
1359 AliTracker::PropagateTrackTo(¶m,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
1360 AliTPCLaserTrack ltr;
1361 AliTPCLaserTrack *ltrp=0x0;
1362 // AliTPCLaserTrack *ltrpjw=0x0;
1364 Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
1365 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
1366 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1368 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1369 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1373 if (id<0) return -1;
1374 if (ltrp->GetSide()!=side) return -1;
1380 if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1383 Float_t radius=TMath::Abs(ltrp->GetX());
1384 param.Rotate(ltrp->GetAlpha());
1385 AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
1387 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1388 Bool_t accept=kTRUE;
1390 // choose closer track
1392 AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1394 Float_t dist0=GetDistance(param0,ltrp);
1395 Float_t dist1=GetDistance(¶m,ltrp);
1396 if (dist0<dist1) accept=kFALSE;
1400 fClusterCounter[id]=countercl;
1401 fTracksEsdParam.AddAt(param.Clone(),id);
1402 fTracksEsd.AddAt(track,id);
1403 fTracksTPC.AddAt(seed,id);
1410 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1412 // Dump Laser info to the tree
1414 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1415 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1416 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1425 param->GetPxPyPz(pxyz);
1427 ltrp->GetPxPyPz(lpxyz);
1428 Float_t dist3D = GetDistance(param,ltrp);
1429 Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1430 Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1433 if (fStreamLevel>0){
1434 TTreeSRedirector *cstream = GetDebugStreamer();
1435 Int_t time = fESD->GetTimeStamp();
1436 Bool_t accept = AcceptLaser(id);
1438 (*cstream)<<"Track"<<
1440 "run="<<fRun<< // run number
1441 "event="<<fEvent<< // event number
1442 "time="<<fTime<< // time stamp of event
1443 "trigger="<<fTrigger<< // trigger
1444 "mag="<<fMagF<< // magnetic field
1448 "driftA.="<<fFitAside<<
1449 "driftC.="<<fFitCside<<
1453 "distphi="<<distphi<<
1456 "counter="<<fCounter[id]<<
1457 "clcounter="<<fClusterCounter[id]<<
1458 "satur="<<fClusterSatur[id]<<
1482 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1484 // Refit the track with different tracklet models:
1485 // 1. Per ROC using the kalman filter, different edge cuts
1486 // 2. Per ROC linear in y and z
1487 // 3. Per ROC quadratic in y and z
1488 // 4. Per track offset for each sector, linear for each sector, common quadratic
1489 // store x, y, z information for all models and the cluster to calculate the residuals
1492 // The clusters which do not fulfill given criteria are skipped
1494 // Cluters removed from fit
1495 // 1. Extended shape > kShapeCut
1496 // 2. In saturation Max > kMax
1497 // 3. Distance to edge < cutEdge
1499 // Clusters not used for the calibration:
1501 // 1. Extended shape > kShapeCut
1502 // 2. In saturation Max > kMax
1505 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1506 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1507 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1509 AliTPCclusterMI dummyCl;
1512 Int_t kMaxTracklets=2;
1513 Float_t kShapeCut = 1.3;
1514 Float_t kRatioCut = 2.;
1516 Float_t kMax = 900.;
1519 //=============================================//
1520 // Linear Fitters for the Different Approaches //
1521 //=============================================//
1522 //linear fit model in y and z; inner - outer sector, combined with offset
1523 static TLinearFitter fy1I(2,"hyp1");
1524 static TLinearFitter fy1O(2,"hyp1");
1525 static TLinearFitter fz1I(2,"hyp1");
1526 static TLinearFitter fz1O(2,"hyp1");
1527 static TLinearFitter fy1IO(3,"hyp2");
1528 static TLinearFitter fz1IO(3,"hyp2");
1529 //quadratic fit model in y and z; inner - sector
1530 static TLinearFitter fy2I(3,"hyp2");
1531 static TLinearFitter fy2O(3,"hyp2");
1532 static TLinearFitter fz2I(3,"hyp2");
1533 static TLinearFitter fz2O(3,"hyp2");
1534 //common quadratic fit for IROC and OROC in y and z
1535 static TLinearFitter fy4(5,"hyp4");
1536 static TLinearFitter fz4(5,"hyp4");
1546 //=============================//
1547 // Loop over all Tracklet Cuts //
1548 //=============================//
1549 for (Int_t icut=0; icut<fNcuts; icut++){
1550 Float_t xinMin = 2500, xinMax=-90;
1551 Float_t xoutMin=2500, xoutMax=-90;
1552 Float_t msigmaYIn=0;
1553 Float_t msigmaYOut=0;
1554 Float_t mqratioIn=0;
1555 Float_t mqratioOut=0;
1558 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1560 Double_t edgeCutX = fEdgeXcuts[icut];
1561 Double_t edgeCutY = fEdgeYcuts[icut];
1562 Int_t nclCut = (Int_t)fNClCuts[icut];
1563 //=========================//
1564 // Parameters to calculate //
1565 //=========================//
1566 //fit parameter inner, outer tracklet and combined fit
1567 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1568 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1570 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1571 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1572 TVectorD vecy4res(5),vecz4res(5);
1573 TVectorD vecy1resIO(3),vecz1resIO(3);
1574 // cluster and track positions for each row - used for residuals
1575 TVectorD vecgX(159); // global X
1576 TVectorD vecgY(159); // global Y
1577 TVectorD vecgZ(159); // global Z
1579 TVectorD vecX(159); // x is the same for all (row center)
1580 TVectorD vecYkalman(159); // y from kalman fit
1581 TVectorD vecZkalman(159); // z from kalman fit
1582 TVectorD vecY1(159); // y from pol1 fit per ROC
1583 TVectorD vecZ1(159); // z from pol1 fit per ROC
1584 TVectorD vecY1IO(159); // y from pol1 fit per ROC
1585 TVectorD vecZ1IO(159); // z from pol1 fit per ROC
1586 TVectorD vecY2(159); // y from pol2 fit per ROC
1587 TVectorD vecZ2(159); // z from pol2 fit per ROC
1588 TVectorD vecY4(159); // y from sector fit
1589 TVectorD vecZ4(159); // z from sector fit
1590 TVectorD vecClY(159); // y cluster position
1591 TVectorD vecClZ(159); // z cluster position
1592 TVectorD vecSec(159); // sector for each row
1593 TVectorD isReject(159); // flag - cluster to be rejected
1595 Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
1596 Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
1597 Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
1598 Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
1599 Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
1600 Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
1601 Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
1602 Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
1603 Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
1604 Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
1605 Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
1606 Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
1608 Int_t innerSector = -1; // number of inner sector
1609 Int_t outerSector = -1; // number of outer sector
1610 Int_t nclI=0; // number of clusters (inner)
1611 Int_t nclO=0; // number of clusters (outer)
1612 //=================================================//
1613 // Perform the Kalman Fit using the Tracklet Class //
1614 //=================================================//
1615 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1616 TObjArray tracklets=
1617 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1618 kFALSE,nclCut,kMaxTracklets);
1619 // tracklet pointers
1620 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1621 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1622 AliTPCTracklet *tr=0x0;
1623 AliTPCTracklet dummy;
1624 //continue if we didn't find a tracklet
1625 if ( !trInner && !trOuter ) continue;
1626 //================================================================================//
1627 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1628 //================================================================================//
1629 if ( trInner && trOuter ){
1630 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1631 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1635 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1639 if ( !trInner->GetInner() ) continue;
1641 if ( trInner->GetSector()>35 ){
1646 if ( !trOuter->GetInner() ) continue;
1648 if ( trOuter->GetSector()<36 ){
1654 innerSector = trInner->GetSector();
1655 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1656 outerSector = trOuter->GetSector();
1657 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1659 // array of clusters
1660 TClonesArray arrCl("AliTPCclusterMI",159);
1661 arrCl.ExpandCreateFast(159);
1662 //=======================================//
1663 // fill fitters with cluster information //
1664 //=======================================//
1665 AliDebug(3,"Filing Cluster Information");
1670 for (Int_t irow=158;irow>-1;--irow) {
1671 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1672 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1677 Float_t meanY=0, sumY=0;
1678 for (Int_t drow=-1;drow<=1;drow++) {
1679 if (irow+drow<0) continue;
1680 if (irow+drow>158) continue;
1681 AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1682 if (!ccurrent) continue;
1683 Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1684 if ( roc!=innerSector && roc!=outerSector ) continue;
1685 if (ccurrent->GetX()<10) continue;
1686 if (ccurrent->GetY()==0) continue;
1687 meanY+=ccurrent->GetY();
1690 if (sumY>0) meanY/=sumY;
1695 Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1696 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1699 Int_t roc = static_cast<Int_t>(c->GetDetector());
1700 if ( roc!=innerSector && roc!=outerSector ) continue;
1702 //store clusters in clones array
1705 if (c->GetMax()<4) continue; // noise cluster?
1706 if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
1708 vecX[irow] = c->GetX();
1709 vecClY[irow] = c->GetY();
1710 vecClZ[irow] = c->GetZ();
1713 // c->GetGlobalXYZ(gxyz);
1714 // vecgX[irow] = gxyz[0];
1715 // vecgY[irow] = gxyz[1];
1716 // vecgZ[irow] = gxyz[2];
1718 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1719 Double_t y = vecClY[irow];
1720 Double_t z = vecClZ[irow];
1722 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
1723 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1724 Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
1725 if ( roc == innerSector ) {
1726 x4[0]=1; //offset inner - outer sector
1727 x4[1]=x; //slope parameter inner sector
1730 x4[2]=x; //slope parameter outer sector
1732 x4[3]=x*x; //common parabolic shape
1733 if (pedgeX < fEdgeXcuts[icut]) continue;
1734 if (pedgeY < fEdgeYcuts[icut]) continue;
1735 if (c->GetMax()>900) continue; // cluster in saturation
1737 if ( roc==innerSector ){
1738 fy1I.AddPoint(x2,y);
1739 fz1I.AddPoint(x2,z);
1740 fy2I.AddPoint(x2,y);
1741 fz2I.AddPoint(x2,z);
1743 if (c->GetX()<xinMin) xinMin=c->GetX();
1744 if (c->GetX()>xinMax) xinMax=c->GetX();
1746 msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1747 mqratioIn +=c->GetMax()/c->GetQ();
1750 if ( roc==outerSector ){
1751 fy1O.AddPoint(x2,y);
1752 fz1O.AddPoint(x2,z);
1753 fy2O.AddPoint(x2,y);
1754 fz2O.AddPoint(x2,z);
1756 if (c->GetX()<xoutMin) xoutMin=c->GetX();
1757 if (c->GetX()>xoutMax) xoutMax=c->GetX();
1758 msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1759 mqratioOut +=c->GetMax()/c->GetQ();
1764 fy1IO.AddPoint(xIO,y);
1765 fz1IO.AddPoint(xIO,z);
1775 //======================================//
1776 // Evaluate and retrieve fit parameters //
1777 //======================================//
1778 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1780 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1785 fy1I.GetParameters(vecy1resInner);
1786 fz1I.GetParameters(vecz1resInner);
1787 fy2I.GetParameters(vecy2resInner);
1788 fz2I.GetParameters(vecz2resInner);
1789 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1790 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1791 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1792 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1795 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1800 fy1O.GetParameters(vecy1resOuter);
1801 fz1O.GetParameters(vecz1resOuter);
1802 fy2O.GetParameters(vecy2resOuter);
1803 fz2O.GetParameters(vecz2resOuter);
1804 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1805 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1806 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1807 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1810 if ( innerSector>0 && outerSector>0 ){
1811 if (fy4.GetNpoints()>0) {
1813 fy4.GetParameters(vecy4res);
1814 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1816 if (fz4.GetNpoints()>0) {
1818 fz4.GetParameters(vecz4res);
1819 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1821 if (fy1IO.GetNpoints()>0) {
1823 fy1IO.GetParameters(vecy1resIO);
1824 chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1826 if (fz1IO.GetNpoints()>0) {
1828 fz1IO.GetParameters(vecz1resIO);
1829 chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1833 fy4.ClearPoints(); fz4.ClearPoints();
1834 fy1I.ClearPoints(); fy1O.ClearPoints();
1835 fz1I.ClearPoints(); fz1O.ClearPoints();
1836 fy2I.ClearPoints(); fy2O.ClearPoints();
1837 fz2I.ClearPoints(); fz2O.ClearPoints();
1838 fy1IO.ClearPoints(); fz1IO.ClearPoints();
1839 //==============================//
1840 // calculate tracklet positions //
1841 //==============================//
1842 AliDebug(4,"Calculate tracklet positions");
1843 for (Int_t irow=158;irow>-1;--irow) {
1845 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1846 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1850 if (!c) { //no cluster
1853 if (c->GetMax()>kMax){ //saturation
1856 if ( vecSec[irow] == outerSector ) { //extended shape
1857 if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
1858 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
1860 if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
1861 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
1867 if ( vecSec[irow]==-1 ) continue; //no cluster info
1868 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1870 Double_t x=vecX[irow];
1871 Double_t xref=x-133.4;
1873 Double_t yoffInner=0;
1874 Double_t zoffInner=0;
1875 Double_t yoffInner1=0;
1876 Double_t zoffInner1=0;
1877 Double_t yslopeInner=0;
1878 Double_t yslopeOuter=0;
1879 Double_t zslopeInner=0;
1880 Double_t zslopeOuter=0;
1881 //positions of hyperplane fits
1882 if ( vecSec[irow] == outerSector ) {
1884 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1885 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1886 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1887 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1888 yslopeOuter=vecy4res[3];
1889 zslopeOuter=vecz4res[3];
1892 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1893 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1894 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1895 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1896 yoffInner=vecy4res[1];
1897 zoffInner=vecz4res[1];
1898 yoffInner1=vecy1resIO[1];
1899 zoffInner1=vecz1resIO[1];
1900 yslopeInner=vecy4res[2];
1901 zslopeInner=vecz4res[2];
1903 vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1904 vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1905 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1906 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1907 //positions of kalman fits
1908 Double_t gxyz[3],xyz[3];
1909 AliExternalTrackParam *param = 0x0;
1911 param=tr->GetInner();
1913 param->GetXYZ(gxyz);
1914 Float_t bz = AliTracker::GetBz(gxyz);
1915 param->GetYAt(x, bz, xyz[1]);
1916 param->GetZAt(x, bz, xyz[2]);
1917 vecYkalman[irow]=xyz[1];
1918 vecZkalman[irow]=xyz[2];
1925 //=====================================================================//
1926 // write results from the different tracklet fits with debug streamers //
1927 //=====================================================================//
1928 if (fStreamLevel>4){
1929 TTreeSRedirector *cstream = GetDebugStreamer();
1931 Float_t dedx = track->GetdEdx();
1932 (*cstream)<<"FitModels"<<
1933 "run="<<fRun<< // run number
1934 "event="<<fEvent<< // event number
1935 "time="<<fTime<< // time stamp of event
1936 "trigger="<<fTrigger<< // trigger
1937 "mag="<<fMagF<< // magnetic field
1940 "edgeCutX=" << edgeCutX <<
1941 "edgeCutY=" << edgeCutY <<
1942 "nclCut=" << nclCut <<
1943 "innerSector="<< innerSector <<
1944 "outerSector="<< outerSector <<
1947 "Tr.=" << extparam <<
1948 "yPol1In.=" << &vecy1resInner <<
1949 "zPol1In.=" << &vecz1resInner <<
1950 "yPol1InOut.="<< &vecy1resIO <<
1951 "zPol1InOut.="<< &vecz1resIO <<
1952 "yPol2In.=" << &vecy2resInner <<
1953 "zPol2In.=" << &vecz2resInner <<
1954 "yPol1Out.=" << &vecy1resOuter <<
1955 "zPol1Out.=" << &vecz1resOuter <<
1956 "yPol2Out.=" << &vecy2resOuter <<
1957 "zPol2Out.=" << &vecz2resOuter <<
1958 "yInOut.=" << &vecy4res <<
1959 "zInOut.=" << &vecz4res <<
1960 "chi2y1In=" << chi2I1y <<
1961 "chi2z1In=" << chi2I1z <<
1962 "chi2y1InOut="<< chi2IO1y <<
1963 "chi2z1InOut="<< chi2IO1z <<
1964 "chi2y1Out=" << chi2O1y <<
1965 "chi2z1Out=" << chi2O1z <<
1966 "chi2y2In=" << chi2I2y <<
1967 "chi2z2In=" << chi2I2z <<
1968 "chi2y2Out=" << chi2O2y <<
1969 "chi2z2Out=" << chi2O2z <<
1970 "chi2yInOut=" << chi2IOy <<
1971 "chi2zInOut=" << chi2IOz <<
1972 "trletIn.=" << trInner <<
1973 "trletOut.=" << trOuter <<
1976 "xinMin=" << xinMin<<
1977 "xinMax=" << xinMax<<
1978 "xoutMin=" << xoutMin<<
1979 "xoutMax=" << xoutMax<<
1980 "msigmaYIn=" <<msigmaYIn<<
1981 "msigmaYOut=" <<msigmaYOut<<
1982 "mqratioIn=" <<mqratioIn<<
1983 "mqratioOut=" << mqratioOut <<
1988 // wirte residuals information
1989 if (fStreamLevel>5){
1990 TTreeSRedirector *cstream = GetDebugStreamer();
1992 Float_t dedx = track->GetdEdx();
1993 (*cstream)<<"Residuals"<<
1994 "run="<<fRun<< // run number
1995 "event="<<fEvent<< // event number
1996 "time="<<fTime<< // time stamp of event
1997 "trigger="<<fTrigger<< // trigger
1998 "mag="<<fMagF<< // magnetic field
2001 "edgeCutX=" << edgeCutX <<
2002 "edgeCutY=" << edgeCutY <<
2003 "nclCut=" << nclCut <<
2005 "Tr.=" << extparam<<
2008 "vX.=" << &vecgX<< // global x
2009 "vY.=" << &vecgY<< // global y
2010 "vZ.=" << &vecgZ<< // global z
2012 "TrYpol1.=" << &vecY1 <<
2013 "TrZpol1.=" << &vecZ1 <<
2014 "TrYpol2.=" << &vecY2 <<
2015 "TrZpol2.=" << &vecZ2 <<
2016 "TrYpol1InOut.="<< &vecY1IO <<
2017 "TrZpol1InOut.="<< &vecZ1IO <<
2018 "TrYInOut.=" << &vecY4 <<
2019 "TrZInOut.=" << &vecZ4 <<
2020 "ClY.=" << &vecClY <<
2021 "ClZ.=" << &vecClZ <<
2022 "isReject.=" << &isReject<<
2023 "sec.=" << &vecSec <<
2026 "xinMin=" << xinMin<<
2027 "xinMax=" << xinMax<<
2028 "xoutMin=" << xoutMin<<
2029 "xoutMax=" << xoutMax<<
2030 "msigmaYIn=" <<msigmaYIn<<
2031 "msigmaYOut=" <<msigmaYOut<<
2032 "mqratioIn=" <<mqratioIn<<
2033 "mqratioOut=" << mqratioOut <<
2034 "yInOut.=" << &vecy4res <<
2035 "zInOut.=" << &vecz4res <<
2037 "chi2y1In=" << chi2I1y << //
2038 "chi2z1In=" << chi2I1z <<
2039 "chi2y1Out=" << chi2O1y <<
2040 "chi2z1Out=" << chi2O1z <<
2041 "chi2y1InOut="<< chi2IO1y <<
2042 "chi2z1InOut="<< chi2IO1z <<
2043 "chi2y2In=" << chi2I2y <<
2044 "chi2z2In=" << chi2I2z <<
2045 "chi2y2Out=" << chi2O2y <<
2046 "chi2z2Out=" << chi2O2z <<
2047 "chi2yInOut=" << chi2IOy <<
2048 "chi2zInOut=" << chi2IOz <<
2050 "yPol1In.=" << &vecy1resInner <<
2051 "zPol1In.=" << &vecz1resInner <<
2052 "yPol2In.=" << &vecy2resInner <<
2053 "zPol2In.=" << &vecz2resInner <<
2054 "yPol1Out.=" << &vecy1resOuter <<
2055 "zPol1Out.=" << &vecz1resOuter <<
2056 "yPol1InOut.="<< &vecy1resIO <<
2057 "zPol1InOut.="<< &vecz1resIO <<
2058 "yPol2Out.=" << &vecy2resOuter <<
2059 "zPol2Out.=" << &vecz2resOuter <<
2065 //==========================//
2066 // Fill Residual Histograms //
2067 //==========================//
2068 if (!fHisNclIn) MakeFitHistos();
2070 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2071 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2072 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2073 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2074 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2075 TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2076 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2077 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2079 for (Int_t irow=158;irow>-1;--irow) {
2080 if (vecSec[irow]==-1)continue; // no cluster info
2081 if (isReject[irow]>0.5) continue; //
2082 Double_t ycl = vecClY[irow];
2083 Double_t yfit = vecY1[irow];
2084 Double_t yfit2 = vecY2[irow];
2085 Double_t x = vecX[irow];
2086 Double_t yabsbeam = -1000;
2087 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2088 yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2089 else if(innerSector==fBeamSectorInner[id])
2090 yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2092 // Double_t yfit3 = vecY2[irow];
2093 Double_t zcl = vecClZ[irow];
2094 Double_t zfit = vecZ1[irow];
2095 Double_t zfit2 = vecZ2[irow];
2096 //Double_t zfit3 = vecZ2[irow];
2099 // The expressions for zcorrected has been obtained by
2100 // inverting the fits in the FitDriftV() method (ignoring the
2101 // global y dependence for now):
2103 // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2105 // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2108 // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
2110 // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2112 Double_t dzabs = -1000;
2113 Double_t zcorrected = -1000;
2114 if (ltrp->GetSide()==0){
2115 if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2116 // ignore global y dependence for now
2118 if(!fUseFixedDriftV)
2119 zcorrected = (zcl + (*fFitAside)[0] -
2120 (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2122 zcorrected = (zcl + fFixedFitAside0 -
2123 (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2124 // zcorrected = zcl;
2125 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2126 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2127 else if(innerSector==fBeamSectorInner[id])
2128 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2131 if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2133 if(!fUseFixedDriftV)
2134 zcorrected = (zcl - (*fFitCside)[0] +
2135 (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2137 zcorrected = (zcl - fFixedFitCside0 +
2138 (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2140 // zcorrected = zcl;
2141 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2142 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2143 else if(innerSector==fBeamSectorInner[id])
2144 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2148 if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2150 profy->Fill(irow,ycl-yfit);
2151 profy2->Fill(irow,ycl-yfit2);
2153 fHisYAbsErrors->Fill(id);
2154 // profyabs->Fill(irow,-0.99);
2156 profyabs->Fill(irow,ycl-yabsbeam);
2158 // profy3->Fill(irow,ycl-yfit3);
2161 profz->Fill(irow,zcl-zfit);
2162 profz2->Fill(irow,zcl-zfit2);
2163 //profz3->Fill(irow,zcl-zfit3);
2166 fHisZAbsErrors->Fill(id);
2168 profzabs->Fill(irow,dzabs);
2174 // Fill laser fit histograms
2176 Float_t dedx = track->GetdEdx();
2177 if (nclI>20&&nclO>20){
2178 fHisNclIn->Fill(id,nclI); //->Number of clusters inner
2179 fHisNclOut->Fill(id,nclO); //->Number of clusters outer
2180 fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
2182 fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
2183 fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
2184 fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
2186 fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2187 fHisdZfit->Fill(id,fFitZ[id]);
2190 fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
2191 fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
2192 fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
2193 fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
2194 fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
2197 fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
2198 fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
2199 fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
2200 fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
2201 fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
2204 fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
2205 fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
2206 fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
2208 fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
2209 fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
2210 fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
2212 fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
2213 fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
2214 fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
2215 fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
2217 fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
2218 fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
2219 fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
2220 fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
2223 fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
2224 fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
2228 fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
2229 fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
2237 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2239 // Dump information about laser beams
2240 // isOK variable indicates usability of the beam
2241 // Beam is not usable if:
2242 // a. No entries in range (krmsCut0)
2243 // b. Big sperad (krmscut1)
2244 // c. RMSto fit sigma bigger then (kmultiCut)
2245 // d. Too big angular spread
2248 const Float_t krmsCut0=0.001;
2249 const Float_t krmsCut1=0.16;
2250 const Float_t kmultiCut=2;
2251 const Float_t kcutP0=0.002;
2252 AliMagF* magF= dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
2253 Double_t xyz[3]={90,0,10}; // tmp. global position
2254 Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical
2255 Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local
2257 AliTPCcalibLaser *laser = this;
2258 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2259 TF1 fg("fg","gaus");
2260 AliTPCParam * tpcparam = 0;
2261 // start set up for absolute residuals analysis
2263 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2264 tpcparam = calib->GetParameters();
2265 if (!tpcparam) tpcparam = new AliTPCParamSR;
2267 AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2269 Float_t bfield = 0, bz=0;
2272 Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
2273 current = grp->GetL3Current((AliGRPObject::Stats)0);
2274 bfield = polarity*5*current/30000.;
2275 bz = polarity*5*current/30000.;
2276 printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2279 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2280 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2281 TLinearFitter lfabsyInner(2);
2282 lfabsyInner.SetFormula("1 ++ x");
2283 TLinearFitter lfabszInner(2);
2284 lfabszInner.SetFormula("1 ++ x");
2286 TLinearFitter lfabsyOuter(2);
2287 lfabsyOuter.SetFormula("1 ++ x");
2288 TLinearFitter lfabszOuter(2);
2289 lfabszOuter.SetFormula("1 ++ x");
2290 // end set up for absolute residuals analysis
2294 for (Int_t id=0; id<336; id++){
2296 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
2297 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2298 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
2299 TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
2300 TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
2301 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2302 //if (!hisphi) continue;
2303 Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2304 //if (entries<minEntries) continue;
2306 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2308 AliTPCLaserTrack::LoadTracks();
2309 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2311 ltrp->UpdatePoints();
2312 pcstream->GetFile()->cd();
2313 if (hisphi) hisphi->Write();
2314 if (hisphiP) hisphiP->Write();
2315 if (hisZ) hisZ->Write();
2316 if (hisP3) hisP3->Write();
2317 if (hisP4) hisP4->Write();
2319 Float_t meanphi = hisphi->GetMean();
2320 Float_t rmsphi = hisphi->GetRMS();
2322 Float_t meanphiP = hisphiP->GetMean();
2323 Float_t rmsphiP = hisphiP->GetRMS();
2324 Float_t meanZ = hisZ->GetMean();
2325 Float_t rmsZ = hisZ->GetRMS();
2326 if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2327 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2328 Double_t gphi1 = fg.GetParameter(1);
2329 Double_t gphi2 = fg.GetParameter(2);
2330 if (hisphiP->GetRMS()>0)
2331 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2332 Double_t gphiP1 = fg.GetParameter(1);
2333 Double_t gphiP2 = fg.GetParameter(2);
2335 if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2336 hisZ->Fit(&fg,"","");
2337 Double_t gz1 = fg.GetParameter(1);
2338 Double_t gz2 = fg.GetParameter(2);
2340 if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2341 hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2342 Double_t gp31 = fg.GetParameter(1);
2343 Double_t gp32 = fg.GetParameter(2);
2344 Double_t meanp3 = hisP3->GetMean();
2345 Double_t rmsp3 = hisP3->GetRMS();
2347 if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2348 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2349 Double_t gp41 = fg.GetParameter(1);
2350 Double_t gp42 = fg.GetParameter(2);
2351 Double_t meanp4 = hisP4->GetMean();
2352 Double_t rmsp4 = hisP4->GetRMS();
2354 Float_t meanS=hisS->GetMean();
2359 ltrp->GetPxPyPz(lpxyz);
2361 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2362 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2363 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2364 if (gphiP2>kcutP0) isOK=kFALSE;
2370 his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2371 Float_t mnclIn = his->GetMean();
2373 his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2374 Float_t mnclOut = his->GetMean();
2376 his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2377 Float_t mnclIO = his->GetMean();
2379 his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2380 Float_t mLclIn = his->GetMean();
2382 his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2383 Float_t mLclOut = his->GetMean();
2385 his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2386 Float_t mLclIO = his->GetMean();
2389 his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2390 Float_t mdEdx = his->GetMean();
2396 his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2397 Float_t mChi2YIn1= his->GetMean();
2399 his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2400 Float_t mChi2YOut1 = his->GetMean();
2402 his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2403 Float_t mChi2YIn2 = his->GetMean();
2405 his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2406 Float_t mChi2YOut2 = his->GetMean();
2408 his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2409 Float_t mChi2YIO1 = his->GetMean();
2411 his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2412 Float_t mChi2ZIn1 = his->GetMean();
2414 his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2415 Float_t mChi2ZOut1 = his->GetMean();
2417 his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2418 Float_t mChi2ZIn2 = his->GetMean();
2420 his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2421 Float_t mChi2ZOut2 = his->GetMean();
2423 his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2424 Float_t mChi2ZIO1 = his->GetMean();
2429 his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2430 Float_t edZfit = his->GetEntries();
2431 Float_t mdZfit = his->GetMean();
2432 Float_t rdZfit = his->GetRMS();
2435 his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2436 Float_t ePy1vP0 = his->GetEntries();
2437 Float_t mPy1vP0 = his->GetMean();
2438 Float_t rPy1vP0 = his->GetRMS();
2441 his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2442 Float_t ePy2vP0 = his->GetEntries();
2443 Float_t mPy2vP0 = his->GetMean();
2444 Float_t rPy2vP0 = his->GetRMS();
2447 his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2448 Float_t ePy3vP0 = his->GetEntries();
2449 Float_t mPy3vP0 = his->GetMean();
2450 Float_t rPy3vP0 = his->GetRMS();
2453 his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2454 Float_t ePy1vP1 = his->GetEntries();
2455 Float_t mPy1vP1 = his->GetMean();
2456 Float_t rPy1vP1 = his->GetRMS();
2459 his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2460 Float_t ePy2vP1 = his->GetEntries();
2461 Float_t mPy2vP1 = his->GetMean();
2462 Float_t rPy2vP1 = his->GetRMS();
2465 his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2466 Float_t ePy3vP1 = his->GetEntries();
2467 Float_t mPy3vP1 = his->GetMean();
2468 Float_t rPy3vP1 = his->GetRMS();
2471 his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2472 Float_t ePy2vP2In = his->GetEntries();
2473 Float_t mPy2vP2In = his->GetMean();
2474 Float_t rPy2vP2In = his->GetRMS();
2477 his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2478 Float_t ePy2vP2Out = his->GetEntries();
2479 Float_t mPy2vP2Out = his->GetMean();
2480 Float_t rPy2vP2Out = his->GetRMS();
2483 his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2484 Float_t ePy3vP2IO = his->GetEntries();
2485 Float_t mPy3vP2IO = his->GetMean();
2486 Float_t rPy3vP2IO = his->GetRMS();
2491 his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2492 Float_t ePz1vP0 = his->GetEntries();
2493 Float_t mPz1vP0 = his->GetMean();
2494 Float_t rPz1vP0 = his->GetRMS();
2497 his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2498 Float_t ePz2vP0 = his->GetEntries();
2499 Float_t mPz2vP0 = his->GetMean();
2500 Float_t rPz2vP0 = his->GetRMS();
2503 his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2504 Float_t ePz3vP0 = his->GetEntries();
2505 Float_t mPz3vP0 = his->GetMean();
2506 Float_t rPz3vP0 = his->GetRMS();
2509 his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2510 Float_t ePz1vP1 = his->GetEntries();
2511 Float_t mPz1vP1 = his->GetMean();
2512 Float_t rPz1vP1 = his->GetRMS();
2515 his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2516 Float_t ePz2vP1 = his->GetEntries();
2517 Float_t mPz2vP1 = his->GetMean();
2518 Float_t rPz2vP1 = his->GetRMS();
2521 his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2522 Float_t ePz3vP1 = his->GetEntries();
2523 Float_t mPz3vP1 = his->GetMean();
2524 Float_t rPz3vP1 = his->GetRMS();
2527 his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2528 Float_t ePz2vP2In = his->GetEntries();
2529 Float_t mPz2vP2In = his->GetMean();
2530 Float_t rPz2vP2In = his->GetRMS();
2533 his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2534 Float_t ePz2vP2Out = his->GetEntries();
2535 Float_t mPz2vP2Out = his->GetMean();
2536 Float_t rPz2vP2Out = his->GetRMS();
2539 his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2540 Float_t ePz3vP2IO = his->GetEntries();
2541 Float_t mPz3vP2IO = his->GetMean();
2542 Float_t rPz3vP2IO = his->GetRMS();
2545 // Fit absolute laser residuals
2546 TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2547 TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2549 Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2550 Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2552 TVectorD vecX(159); // X
2553 TVectorD vecY(159); // Y
2554 TVectorD vecR(159); // R
2555 TVectorD vecDY(159); // absolute residuals in Y
2556 TVectorD vecDZ(159); // absolute residuals in Z
2557 TVectorD vecN(159); // number of clusters
2558 TVectorD vecEy(159); //error y
2559 TVectorD vecEz(159); //error z
2560 TVectorD vecPhi(159); // local tangent
2561 TVectorD vecPhiR(159); // local tangent
2562 // magnetic field integrals
2563 TVectorD vecIBR(159); // radial
2564 TVectorD vecIBRPhi(159); // r-phi
2565 TVectorD vecIBLX(159); // local x
2566 TVectorD vecIBLY(159); // local y
2567 TVectorD vecIBGX(159); // local x
2568 TVectorD vecIBGY(159); // local y
2569 TVectorD vecIBZ(159); // z
2571 for (Int_t irow=0;irow<159;irow++){
2579 Double_t gx =(*(ltrp->fVecGX))[irow];
2580 Double_t gy =(*(ltrp->fVecGY))[irow];
2581 Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
2582 Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2583 Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2584 xyz[2]=(*(ltrp->fVecGZ))[irow];
2585 xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2586 xyz[1]=TMath::ATan2(gy,gx);
2587 Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2589 magF->GetTPCIntCyl(xyz,bxyz);
2590 magF->GetTPCInt(gxyz,bgxyz);
2591 vecIBR[irow]=bxyz[0];
2592 vecIBRPhi[irow]=bxyz[1];
2594 vecIBGX[irow]=bgxyz[0];
2595 vecIBGY[irow]=bgxyz[1];
2597 vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
2598 vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2601 vecIBZ[irow]=bxyz[2];
2606 lfabsyInner.ClearPoints();
2607 lfabszInner.ClearPoints();
2608 lfabsyOuter.ClearPoints();
2609 lfabszOuter.ClearPoints();
2612 Float_t yAbsInOffset = -100;
2613 Float_t yAbsInSlope = -100;
2614 Float_t yAbsInDeltay = -100;
2616 Float_t zAbsInOffset = -100;
2617 Float_t zAbsInSlope = -100;
2618 Float_t zAbsInDeltaz = -100;
2620 Float_t yAbsOutOffset = -100;
2621 Float_t yAbsOutSlope = -100;
2622 Float_t yAbsOutDeltay = -100;
2624 Float_t zAbsOutOffset = -100;
2625 Float_t zAbsOutSlope = -100;
2626 Float_t zAbsOutDeltaz = -100;
2628 Float_t lasTanPhiLocIn = -100;
2629 Float_t lasTanPhiLocOut = -100;
2631 if(histAbsY->GetEntries()>0) {
2633 Double_t rotAngOut = 10;
2634 Double_t rotAngIn = 10;
2635 if((secInner%36)!=(secOuter%36))
2636 rotAngIn += 20; // 30 degrees
2638 // Calculate laser mirror X position in local frame
2639 Double_t laserposOut =
2640 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2641 Double_t laserposIn =
2642 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2644 // Calculate laser tan phi in local frame
2645 lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2646 if(lasTanPhiLocOut<0) {
2647 lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2648 lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2651 lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2652 lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2655 lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2656 lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2658 TProfile* yprof = histAbsY->ProfileX("yprof");
2659 TProfile* zprof = histAbsZ->ProfileX("zprof");
2661 for(Int_t bin = 1; bin<=159; bin++) {
2663 if(yprof->GetBinEntries(bin)<10&&
2664 zprof->GetBinEntries(bin)<10) {
2668 // There is a problem in defining inner and outer sectors for
2669 // the outer beams (0 and 6) where both sectors are OROCs. To
2670 // make sure there is no overlap row 94 to 99 are cutted.
2671 if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2674 Int_t row = (bin-1);
2678 Bool_t isOuter = kTRUE;
2679 Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2682 (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2685 sector = TMath::Nint(fBeamSectorInner[id]);
2689 Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2690 vecN[bin-1] =yprof->GetBinEntries(bin);
2691 vecEy[bin-1]=yprof->GetBinError(bin);
2692 vecEz[bin-1]=zprof->GetBinError(bin);
2694 vecDY[bin-1] = yprof->GetBinContent(bin);
2695 vecDZ[bin-1] = zprof->GetBinContent(bin);
2696 if(!isOuter) { // inner
2697 vecPhi[bin-1]=lasTanPhiLocIn;
2698 // Calculate local y from residual and database
2699 Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2702 Double_t r = TMath::Sqrt(x*x + y*y);
2704 // Find angle between laser vector and R vector
2705 // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2706 Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2707 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2708 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2709 if(lasTanPhiLocIn<0)
2710 vecPhiR[bin-1]*=-1; // to have the same sign
2712 if(yprof->GetBinEntries(bin)>=10) {
2713 lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2714 TMath::Max(yprof->GetBinError(bin), 0.001));
2716 if(zprof->GetBinEntries(bin)>=10) {
2717 lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2718 TMath::Max(zprof->GetBinError(bin), 0.001));
2721 vecPhi[bin-1]=lasTanPhiLocOut;
2722 // Calculate local y from residual and database
2723 Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2726 Double_t r = TMath::Sqrt(x*x + y*y);
2729 Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2730 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2731 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2732 if(lasTanPhiLocOut<0)
2733 vecPhiR[bin-1]*=-1; // to have the same sign
2735 if(yprof->GetBinEntries(bin)>=10) {
2736 lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2737 TMath::Max(yprof->GetBinError(bin), 0.001));
2739 if(zprof->GetBinEntries(bin)>=10) {
2740 lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2741 TMath::Max(zprof->GetBinError(bin), 0.001));
2748 delete yprof; delete zprof;
2751 // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2752 nClInY = lfabsyInner.GetNpoints();
2753 if(lfabsyInner.GetNpoints()>10) {
2754 lfabsyInner.EvalRobust(0.95);
2757 lfabsyInner.GetParameters(result);
2758 yAbsInOffset = result[0];
2759 yAbsInSlope = result[1];
2760 yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2762 nClInZ = lfabszInner.GetNpoints();
2763 if(lfabszInner.GetNpoints()>10) {
2764 lfabszInner.EvalRobust(0.95);
2767 lfabszInner.GetParameters(result);
2768 zAbsInOffset = result[0];
2769 zAbsInSlope = result[1];
2770 zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2772 nClOutY = lfabsyOuter.GetNpoints();
2773 if(lfabsyOuter.GetNpoints()>10) {
2774 lfabsyOuter.EvalRobust(0.95);
2777 lfabsyOuter.GetParameters(result);
2778 yAbsOutOffset = result[0];
2779 yAbsOutSlope = result[1];
2780 yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2782 nClOutZ = lfabszOuter.GetNpoints();
2783 if(lfabszOuter.GetNpoints()>10) {
2784 lfabszOuter.EvalRobust(0.95);
2787 lfabszOuter.GetParameters(result);
2788 zAbsOutOffset = result[0];
2789 zAbsOutSlope = result[1];
2790 zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2796 Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2797 Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2798 Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2799 Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2801 Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2802 Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2804 Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2805 Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2808 (*pcstream)<<"Mean"<<
2821 "id="<<id<< // track id
2822 "entries="<<entries<< // number of entries
2823 "bz="<<bfield<< // bfield
2824 "LTr.="<<ltrp<< // refernece track
2826 "mnclIn="<<mnclIn<< // mean number of clusters in inner
2827 "mnclOut="<<mnclOut<< // mean number of clusters in outer
2828 "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2829 "mLclIn="<<mLclIn<< // mean number of clusters in inner
2830 "mLclOut="<<mLclOut<< // mean number of clusters in outer
2831 "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2832 "mdEdx="<<mdEdx<< // mean dEdx
2833 "edZfit="<<edZfit<< // entries z fit
2834 "mdZfit="<<mdZfit<< // mean z fit
2835 "rdZfit="<<rdZfit<< // RMS z fit
2838 "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2839 "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2840 "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2841 "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2842 "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2843 "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2844 "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2845 "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2846 "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2847 "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2850 "ePy1vP0="<<ePy1vP0<<
2851 "mPy1vP0="<<mPy1vP0<<
2852 "rPy1vP0="<<rPy1vP0<<
2853 "ePy2vP0="<<ePy2vP0<<
2854 "mPy2vP0="<<mPy2vP0<<
2855 "rPy2vP0="<<rPy2vP0<<
2856 "ePy3vP0="<<ePy3vP0<<
2857 "mPy3vP0="<<mPy3vP0<<
2858 "rPy3vP0="<<rPy3vP0<<
2859 "ePy1vP1="<<ePy1vP1<<
2860 "mPy1vP1="<<mPy1vP1<<
2861 "rPy1vP1="<<rPy1vP1<<
2862 "ePy2vP1="<<ePy2vP1<<
2863 "mPy2vP1="<<mPy2vP1<<
2864 "rPy2vP1="<<rPy2vP1<<
2865 "ePy3vP1="<<ePy3vP1<<
2866 "mPy3vP1="<<mPy3vP1<<
2867 "rPy3vP1="<<rPy3vP1<<
2868 "ePy2vP2In="<<ePy2vP2In<<
2869 "mPy2vP2In="<<mPy2vP2In<<
2870 "rPy2vP2In="<<rPy2vP2In<<
2871 "ePy2vP2Out="<<ePy2vP2Out<<
2872 "mPy2vP2Out="<<mPy2vP2Out<<
2873 "rPy2vP2Out="<<rPy2vP2Out<<
2874 "ePy3vP2IO="<<ePy3vP2IO<<
2875 "mPy3vP2IO="<<mPy3vP2IO<<
2876 "rPy3vP2IO="<<rPy3vP2IO<<
2879 "ePz1vP0="<<ePz1vP0<<
2880 "mPz1vP0="<<mPz1vP0<<
2881 "rPz1vP0="<<rPz1vP0<<
2882 "ePz2vP0="<<ePz2vP0<<
2883 "mPz2vP0="<<mPz2vP0<<
2884 "rPz2vP0="<<rPz2vP0<<
2885 "ePz3vP0="<<ePz3vP0<<
2886 "mPz3vP0="<<mPz3vP0<<
2887 "rPz3vP0="<<rPz3vP0<<
2888 "ePz1vP1="<<ePz1vP1<<
2889 "mPz1vP1="<<mPz1vP1<<
2890 "rPz1vP1="<<rPz1vP1<<
2891 "ePz2vP1="<<ePz2vP1<<
2892 "mPz2vP1="<<mPz2vP1<<
2893 "rPz2vP1="<<rPz2vP1<<
2894 "ePz3vP1="<<ePz3vP1<<
2895 "mPz3vP1="<<mPz3vP1<<
2896 "rPz3vP1="<<rPz3vP1<<
2897 "ePz2vP2In="<<ePz2vP2In<<
2898 "mPz2vP2In="<<mPz2vP2In<<
2899 "rPz2vP2In="<<rPz2vP2In<<
2900 "ePz2vP2Out="<<ePz2vP2Out<<
2901 "mPz2vP2Out="<<mPz2vP2Out<<
2902 "rPz2vP2Out="<<rPz2vP2Out<<
2903 "ePz3vP2IO="<<ePz3vP2IO<<
2904 "mPz3vP2IO="<<mPz3vP2IO<<
2905 "rPz3vP2IO="<<rPz3vP2IO<<
2909 "lx0="<<lxyz[0]<< // reference x
2910 "lx1="<<lxyz[1]<< // reference y
2911 "lx2="<<lxyz[2]<< // refernece z
2912 "lpx0="<<lpxyz[0]<< // reference x
2913 "lpx1="<<lpxyz[1]<< // reference y
2914 "lpx2="<<lpxyz[2]<< // refernece z
2918 "mphi="<<meanphi<< //
2919 "rmsphi="<<rmsphi<< //
2923 "mphiP="<<meanphiP<< //
2924 "rmsphiP="<<rmsphiP<< //
2933 "gp31="<<gp31<< //gaus mean - tgl
2934 "gp32="<<gp32<< //gaus rms -tgl
2935 "meanp3="<<meanp3<< //mean - tgl
2936 "rmsp3="<<rmsp3<< //rms -tgl
2937 "gp41="<<gp41<< //gaus mean - P4
2938 "gp42="<<gp42<< //gaus rms - P4
2939 "meanp4="<<meanp4<< //mean - P4
2940 "rmsp4="<<rmsp4<< //rms - P4
2941 // Parameters from abs res analysis
2942 "SecIn="<<secInner<< // inner sector
2943 "SecOut="<<secOuter<< // outer sector
2944 "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2945 "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
2946 "ibr.="<<&vecIBR<< // radial filed integral
2947 "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
2948 "ibr.="<<&vecIBR<< // radial filed integral
2949 "ibz.="<<&vecIBZ<< // z filed integral
2951 "iblx.="<<&vecIBLX<< // local bx integral
2952 "ibly.="<<&vecIBLY<< // local by integral
2953 "ibgx.="<<&vecIBGX<< // global bx integral
2954 "ibgy.="<<&vecIBGY<< // global by integral
2956 "X.="<<&vecX<< // local x
2957 "Y.="<<&vecY<< // local y
2958 "R.="<<&vecR<< // radius
2959 "dY.="<<&vecDY<< // abs y residuals
2960 "dZ.="<<&vecDZ<< // abs z residuals
2961 "eY.="<<&vecEy<< // error of y residuals
2962 "eZ.="<<&vecEz<< // error of z residuals
2963 "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
2964 "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
2965 "nCl.="<<&vecN<< // number of clusters
2967 "nClInY="<<nClInY<< // Number of clusters for inner y
2968 "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2969 "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
2970 "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2971 "nClInZ="<<nClInZ<< // Number of clusters for inner z
2972 "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2973 "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
2974 "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2976 "nClOutY="<<nClOutY<< // Number of clusters for outer y
2977 "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
2978 "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
2979 "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
2980 "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
2981 "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
2982 "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
2983 "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
2990 TFile fmean("laserMean.root");
3000 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3004 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3005 TFile * f = pcstream->GetFile();
3007 f->mkdir("dirphiP");
3009 TF1 fp("p1","pol1");
3014 char grnamefull[1000];
3017 Double_t mphiP[100];
3018 Double_t smphi[100];
3019 Double_t smphiP[100];
3030 for (Int_t id=0; id<336; id++){
3032 sprintf(cut,"fId==%d&&%s",id,cutUser);
3033 Int_t entries = chain->Draw("bz",cut,"goff");
3034 if (entries<3) continue;
3035 AliTPCLaserTrack *ltrp = 0;
3036 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3037 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3041 ltrp->GetPxPyPz(lpxyz);
3043 chain->Draw("bz",cut,"goff");
3044 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3045 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3046 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3048 chain->Draw("gphi1",cut,"goff");
3049 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3050 chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3051 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3053 chain->Draw("gphiP1",cut,"goff");
3054 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3055 chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3056 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3058 chain->Draw("gz1",cut,"goff");
3059 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3060 chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3061 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3064 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3065 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3072 for (Int_t ientry=0; ientry<entries; ientry++){
3073 if (TMath::Abs(bz[ientry])<0.05){
3074 phiB0 = mphi[ientry];
3075 phiPB0 = mphiP[ientry];
3079 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3082 pphi[0] = fp.GetParameter(0); // offset
3083 pphi[1] = fp.GetParameter(1); // slope
3084 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3085 sprintf(grname,"phi_id%d",id);
3086 grphi->SetName(grname); grphi->SetTitle(grnamefull);
3087 grphi->GetXaxis()->SetTitle("b_{z} (T)");
3088 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3089 grphi->SetMaximum(1.2);
3090 grphi->SetMinimum(-1.2);
3094 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3097 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3100 pphiP[0] = fp.GetParameter(0); // offset
3101 pphiP[1] = fp.GetParameter(1); // slope
3102 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3103 sprintf(grname,"phiP_id%d",id);
3104 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3105 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3106 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3107 grphiP->SetMaximum(pphiP[0]+0.005);
3108 grphiP->SetMinimum(pphiP[0]-0.005);
3110 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3115 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3118 pmZ[0] = fp.GetParameter(0); // offset
3119 pmZ[1] = fp.GetParameter(1); // slope
3120 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3121 sprintf(grname,"mZ_id%d",id);
3122 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3123 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3124 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3126 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3132 for (Int_t ientry=0; ientry<entries; ientry++){
3133 (*pcstream)<<"Mean"<<
3136 "entries="<<entries<<
3138 "lx0="<<lxyz[0]<< // reference x
3139 "lx1="<<lxyz[1]<< // reference y
3140 "lx2="<<lxyz[2]<< // refernece z
3141 "lpx0="<<lpxyz[0]<< // reference x
3142 "lpx1="<<lpxyz[1]<< // reference y
3143 "lpx2="<<lpxyz[2]<< // refernece z
3145 "phiB0="<<phiB0<< // position shift at 0 field
3146 "phiPB0="<<phiPB0<< // angular shift at 0 field
3147 "zB0="<<zB0<< // z shift for 0 field
3149 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3150 "pphi0="<<pphi[0]<< // offset
3151 "pphi1="<<pphi[1]<< // slope
3152 "pphi2="<<pphi[2]<< // norm chi2
3154 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3155 "pphiP0="<<pphiP[0]<< // offset
3156 "pphiP1="<<pphiP[1]<< // slope
3157 "pphiP2="<<pphiP[2]<< // norm chi2
3159 "gz1="<<mZ[ientry]<<
3160 "pmZ0="<<pmZ[0]<< // offset
3161 "pmZ1="<<pmZ[1]<< // slope
3162 "pmZ2="<<pmZ[2]<< // norm chi2
3172 void AliTPCcalibLaser::Analyze(){
3179 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3181 TIterator* iter = li->MakeIterator();
3182 AliTPCcalibLaser* cal = 0;
3183 static Int_t counter0=0;
3184 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3185 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3186 // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3189 printf("Marging number %d\n", counter0);
3192 MergeFitHistos(cal);
3197 // TProfile *hp=0x0;
3198 //TProfile *hpm=0x0;
3200 for (Int_t id=0; id<336; id++){
3201 // merge fDeltaZ histograms
3202 hm = (TH1F*)cal->fDeltaZ.At(id);
3203 h = (TH1F*)fDeltaZ.At(id);
3205 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3207 fDeltaZ.AddAt(h,id);
3210 // merge fP3 histograms
3211 hm = (TH1F*)cal->fDeltaP3.At(id);
3212 h = (TH1F*)fDeltaP3.At(id);
3214 h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3216 fDeltaP3.AddAt(h,id);
3219 // merge fP4 histograms
3220 hm = (TH1F*)cal->fDeltaP4.At(id);
3221 h = (TH1F*)fDeltaP4.At(id);
3223 h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3225 fDeltaP4.AddAt(h,id);
3230 // merge fDeltaPhi histograms
3231 hm = (TH1F*)cal->fDeltaPhi.At(id);
3232 h = (TH1F*)fDeltaPhi.At(id);
3234 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3236 fDeltaPhi.AddAt(h,id);
3239 // merge fDeltaPhiP histograms
3240 hm = (TH1F*)cal->fDeltaPhiP.At(id);
3241 h = (TH1F*)fDeltaPhiP.At(id);
3243 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3245 fDeltaPhiP.AddAt(h,id);
3248 // merge fSignals histograms
3249 hm = (TH1F*)cal->fSignals.At(id);
3250 h = (TH1F*)fSignals.At(id);
3252 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3254 fSignals.AddAt(h,id);
3259 // merge ProfileY histograms -0
3260 h2m = (TH2F*)cal->fDeltaYres.At(id);
3261 h2 = (TH2F*)fDeltaYres.At(id);
3262 if (h2m&&h2) h2->Add(h2m);
3264 h2m = (TH2F*)cal->fDeltaZres.At(id);
3265 h2 = (TH2F*)fDeltaZres.At(id);
3266 if (h2m&&h2) h2->Add(h2m);
3267 // merge ProfileY histograms - 2
3268 h2m = (TH2F*)cal->fDeltaYres2.At(id);
3269 h2 = (TH2F*)fDeltaYres2.At(id);
3270 if (h2m&&h2) h2->Add(h2m);
3272 h2m = (TH2F*)cal->fDeltaZres2.At(id);
3273 h2 = (TH2F*)fDeltaZres2.At(id);
3274 if (h2m&&h2) h2->Add(h2m);
3276 // merge ProfileY histograms - abs
3277 h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3278 h2 = (TH2F*)fDeltaYresAbs.At(id);
3279 if (h2m&&h2) h2->Add(h2m);
3281 h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3282 h2 = (TH2F*)fDeltaZresAbs.At(id);
3283 if (h2m&&h2) h2->Add(h2m);
3284 // merge ProfileY histograms - 3
3285 //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3286 //h2 = (TH2F*)fDeltaYres3.At(id);
3287 //if (h2m) h2->Add(h2m);
3289 //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3290 //h2 = (TH2F*)fDeltaZres3.At(id);
3291 //if (h2m) h->Add(h2m);
3299 void AliTPCcalibLaser::MakeFitHistos(){
3301 // Make a fit histograms
3303 // Number of clusters
3305 //TH2F *fHisNclIn; //->Number of clusters inner
3306 //TH2F *fHisNclOut; //->Number of clusters outer
3307 //TH2F *fHisNclIO; //->Number of cluster inner outer
3308 // TH2F *fHisdEdx; //->dEdx histo
3309 fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3310 fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3311 fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3313 fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3314 fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3315 fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3317 fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3318 fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3323 // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3324 // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3325 // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3326 // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3327 // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3328 fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3329 fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3330 fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3331 fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3332 fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3333 // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3334 // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3335 // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3336 // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3337 // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3338 fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3339 fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3340 fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3341 fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3342 fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3347 // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3348 // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3349 // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3350 // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3351 // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3352 // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3353 // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3354 // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3355 // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3356 fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3357 fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3358 fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3359 fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3360 fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3361 fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3362 fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3363 fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3364 fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3367 // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3368 // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3369 // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3370 // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3371 // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3372 // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3373 // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3374 // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3375 // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3376 fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3377 fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3378 fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3379 fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3380 fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3381 fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3382 fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3383 fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3384 fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3386 fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3387 fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3389 fHisNclIn->SetDirectory(0); //->Number of clusters inner
3390 fHisNclOut->SetDirectory(0); //->Number of clusters outer
3391 fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3392 fHisLclIn->SetDirectory(0); //->Level arm inner
3393 fHisLclOut->SetDirectory(0); //->Level arm outer
3394 fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3395 fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3396 fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3399 fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3400 fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3401 fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3402 fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3403 fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3404 fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3405 fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3406 fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3407 fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3408 fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3411 fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3412 fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3413 fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3414 fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3415 fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3416 fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3417 fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3418 fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3419 fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3422 fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3423 fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3424 fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3425 fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3426 fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3427 fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3428 fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3429 fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3430 fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3432 fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3433 fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3440 for (Int_t id=0; id<336;id++){
3441 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3442 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3443 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3444 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
3445 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3446 TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3447 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3448 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3450 profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3451 profy->SetDirectory(0);
3452 fDeltaYres.AddAt(profy,id);
3453 profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3454 profy2->SetDirectory(0);
3455 fDeltaYres2.AddAt(profy2,id);
3456 if(!fUseFixedDriftV)
3457 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
3459 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
3460 profyabs->SetDirectory(0);
3461 fDeltaYresAbs.AddAt(profyabs,id);
3462 //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3463 //profy3->SetDirectory(0);
3464 //fDeltaYres3.AddAt(profy3,id);
3467 profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3468 profz->SetDirectory(0);
3469 fDeltaZres.AddAt(profz,id);
3470 profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3471 profz2->SetDirectory(0);
3472 fDeltaZres2.AddAt(profz2,id);
3473 if(!fUseFixedDriftV)
3474 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
3476 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
3477 profzabs->SetDirectory(0);
3478 fDeltaZresAbs.AddAt(profzabs,id);
3479 //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3480 //profz3->SetDirectory(0);
3481 //fDeltaZres3.AddAt(profz3,id);
3486 for (Int_t id=0; id<336;id++){
3487 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3488 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3489 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3491 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3492 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3493 TH1F * hisSignal = (TH1F*)fSignals.At(id);
3496 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3497 hisdz->SetDirectory(0);
3498 fDeltaZ.AddAt(hisdz,id);
3500 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3501 hisP3->SetDirectory(0);
3502 fDeltaP3.AddAt(hisP3,id);
3504 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3505 hisP4->SetDirectory(0);
3506 fDeltaP4.AddAt(hisP4,id);
3509 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3510 hisdphi->SetDirectory(0);
3511 fDeltaPhi.AddAt(hisdphi,id);
3513 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3514 hisdphiP->SetDirectory(0);
3515 fDeltaPhiP.AddAt(hisdphiP,id);
3516 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3517 hisSignal->SetDirectory(0);
3518 fSignals.AddAt(hisSignal,id);
3522 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3523 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3524 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3525 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3529 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3530 Int_t binsLaser[12]= {336, //id
3542 Double_t xminLaser[12]= {0, //id
3554 Double_t xmaxLaser[12]= {336, //id
3567 TString nameLaser[12]= {"id",
3579 TString titleLaser[12]= {"id",
3591 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3592 for (Int_t iaxis=1; iaxis<12; iaxis++){
3593 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3594 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3598 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3600 // Merge content of histograms
3602 // Only first histogram is checked - all other should be the same
3603 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3605 if (!laser->fHisNclIn) return; // empty histograms
3606 if (!fHisNclIn) MakeFitHistos();
3608 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3609 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3610 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3611 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3612 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3613 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3614 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3615 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3618 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3619 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3620 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3621 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3622 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3623 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3624 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3625 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3626 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3627 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3630 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3631 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3632 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3633 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3634 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3635 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3636 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3637 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3638 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3641 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3642 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3643 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3644 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3645 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3646 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3647 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3648 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3649 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3650 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3651 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3664 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3666 // Dump fit information - collect information from the streamers
3669 TChain * chainFit=0;
3670 TChain * chainTrack=0;
3673 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3674 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3675 AliXRDPROOFtoolkit tool;
3676 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3677 chainTrack->Lookup();
3678 chainTrack->SetProof(kTRUE);
3679 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3680 chainDrift->Lookup();
3681 chainDrift->SetProof(kTRUE);
3683 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3685 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3687 chainFit->SetProof(kTRUE);
3688 chain->SetProof(kTRUE);
3689 AliTPCLaserTrack::LoadTracks();
3690 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3696 TCut cutP3("abs(Tr.fP[3])<0.1");
3697 TCut cutP4("abs(Tr.fP[4])<0.5");
3698 TCut cutPx = cutP3+cutP4;
3699 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3700 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3701 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3702 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3704 TCut cutdEdx("sqrt(dEdx)>3");
3705 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3706 TCut cutN("nclO>20&&nclI>20");
3707 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3711 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3712 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3713 TCut cutClX("abs(Cl[].fX)>10");
3714 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3715 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3716 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3717 TCut cutQ("sqrt(Cl[].fMax)>4");
3718 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3722 TH1F * phisAccept = 0;
3724 TProfile * pdEdx = 0;
3732 TProfile * pNclI = 0;
3733 TProfile * pNclO = 0;
3735 TProfile * pchi2YIn =0;
3736 TProfile * pchi2ZIn =0;
3737 TProfile * pchi2YOut =0;
3738 TProfile * pchi2ZOut =0;
3739 TProfile * pchi2YInOut =0;
3740 TProfile * pchi2ZInOut =0;;
3742 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3743 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3744 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3745 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3746 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3747 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3749 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3750 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3753 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3754 pP0 = (TProfile*)gROOT->FindObject("hP0");
3755 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3756 pP1 = (TProfile*)gROOT->FindObject("hP1");
3757 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3758 pP2 = (TProfile*)gROOT->FindObject("hP2");
3759 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3760 pP3 = (TProfile*)gROOT->FindObject("hP3");
3761 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3762 pP4 = (TProfile*)gROOT->FindObject("hP4");
3764 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3765 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3766 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3767 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3770 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3771 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3772 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3773 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3774 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3775 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3776 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3777 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3778 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3779 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3780 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3781 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3783 // second derivatives
3785 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3786 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3787 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3788 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3789 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3790 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3792 phisPy2In->FitSlicesY();
3793 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3794 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3795 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3797 phisPy2Out->FitSlicesY();
3798 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3799 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3800 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3802 phisPy2InOut->FitSlicesY();
3803 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3804 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3805 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3807 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3808 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3809 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3810 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3811 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3812 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3814 phisPz2In->FitSlicesY();
3815 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3816 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3817 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3819 phisPz2Out->FitSlicesY();
3820 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3821 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3822 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3824 phisPz2InOut->FitSlicesY();
3825 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3826 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3827 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3834 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3835 for (Int_t ilaser=0; ilaser<336; ilaser++){
3836 Float_t all=phisAl->GetBinContent(ilaser+1);
3837 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3838 Float_t out=phisOut->GetBinContent(ilaser+1);
3839 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3840 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3841 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3842 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3843 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3844 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3847 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3848 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3850 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3851 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3852 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3853 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3854 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3855 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3857 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3858 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3859 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3861 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3862 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3863 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3865 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3866 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3867 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3869 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3870 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3871 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3873 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3874 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3875 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3877 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3878 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3879 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3881 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3882 (*pcstream)<<"Scan"<<
3896 "chi2YIn="<<chi2YIn<<
3897 "chi2YOut="<<chi2YOut<<
3898 "chi2YInOut="<<chi2YInOut<<
3899 "chi2ZIn="<<chi2ZIn<<
3900 "chi2ZOut="<<chi2ZOut<<
3901 "chi2ZInOut="<<chi2ZInOut<<
3903 "nPy2In="<<entriesPy2In<<
3904 "mPy2In="<<meanPy2In<<
3905 "sPy2In="<<sigmaPy2In<<
3907 "nPy2Out="<<entriesPy2Out<<
3908 "mPy2Out="<<meanPy2Out<<
3909 "sPy2Out="<<sigmaPy2Out<<
3911 "nPy2InOut="<<entriesPy2InOut<<
3912 "mPy2InOut="<<meanPy2InOut<<
3913 "sPy2InOut="<<sigmaPy2InOut<<
3915 "nPz2In="<<entriesPz2In<<
3916 "mPz2In="<<meanPz2In<<
3917 "sPz2In="<<sigmaPz2In<<
3919 "nPz2Out="<<entriesPz2Out<<
3920 "mPz2Out="<<meanPz2Out<<
3921 "sPz2Out="<<sigmaPz2Out<<
3923 "nPz2InOut="<<entriesPz2InOut<<
3924 "mPz2InOut="<<meanPz2InOut<<
3925 "sPz2InOut="<<sigmaPz2InOut<<
3933 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3934 TVectorD& meanSlope,
3935 TVectorD& sectorArray,
3938 // This method should ideally go in AliTPCLaser
3939 // option == 0 (pads outer - closest to beam)
3940 // option == 1 (pads inner)
3941 // option == 2 (time outer)
3942 // option == 3 (time inner)
3943 Int_t nFailures = 0;
3945 for(Int_t id = 0; id < 336; id++) {
3947 if (!AliTPCLaserTrack::GetTracks())
3948 AliTPCLaserTrack::LoadTracks();
3949 AliTPCLaserTrack *ltrp =
3950 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3952 AliExternalTrackParam trackParam(*ltrp);
3954 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3955 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3956 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3958 Double_t angle = trackParam.GetAlpha();
3960 angle += 2*TMath::Pi();
3961 if(trackParam.GetSnp()>0) // track points to sector "before"
3962 angle -= deltaangle*TMath::DegToRad();
3963 else // track points to sector "after"
3964 angle += deltaangle*TMath::DegToRad();
3966 Bool_t success = trackParam.Rotate(angle);
3969 // cout << "WARNING: Rotate failed for ID: " << id << endl;
3973 angle *= TMath::RadToDeg();
3975 Int_t sector = TMath::Nint((angle-10.0)/20.0);
3980 if(ltrp->GetSide()==1) // C side
3982 if(option==0 || option==2)
3984 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
3987 sectorArray[id] = sector;
3989 const Double_t x0 = 0;
3991 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3992 Double_t slopez = trackParam.GetTgl();
3993 // One needs a factor sqrt(1+slopey**2) to take into account the
3994 // longer path length
3995 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3996 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3998 // Double_t offsetz = trackParam.GetZ();
3999 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4000 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4001 if(option==2 || option==3) {
4002 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4004 meanOffset[id] = offsety; meanSlope[id] = slopey;
4009 AliWarning(Form("Rotate method failed %d times", nFailures));
4014 TFile f("vscan.root");
4019 chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
4021 chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
4024 chain->Draw("Cl.fY-TrYpol1.fElements-AliTPCClusterParam::SPosCorrection(0,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
4027 chain->Draw("Cl[].fZ-TrZpol1.fElements-0*AliTPCClusterParam::SPosCorrection(1,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fTimeBin-int(Cl[].fTimeBin)",cutA+cutCl+"Cl[].fZ>0","prof",10000)
4036 // check edge effects
4037 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4039 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4041 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
4045 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4046 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")
4057 dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
4060 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):pi/18-abs(Cl.fY/Cl.fX)>>hisYdphi(100,0,0.03)",""+cutA+cutCl,"prof",10000)
4062 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdy(100,0,5)",""+cutA+cutCl,"prof",10000)
4068 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyIROC(100,0,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
4070 chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyOROC(100,0,5)","Cl.fDetector>36"+cutA+cutCl,"prof",100000)
4074 chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>his(100,-5,5)",""+cutA+cutCl,"prof",100000)
4076 chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>hisdyInner(100,-5,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
4085 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutDY,"prof")
4087 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&<r.fP[1]<0"+cutA+cutDY,"prof")
4091 chainFit->Draw("LTr.fId","nclI>10",100000)
4093 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4095 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4097 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4107 gSystem->Load("libSTAT.so")
4108 TStatToolkit toolkit;
4114 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
4119 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
4120 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
4121 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
4122 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
4124 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
4125 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
4126 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
4127 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
4129 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
4130 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
4131 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
4132 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
4137 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4139 treeT->SetAlias("fit",strq0->Data());
4142 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4144 treeT->SetAlias("fitP",strqP->Data());
4147 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4148 treeT->SetAlias("fitD",strqDrift->Data());
4151 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4153 for (Int_t i=0; i<6;i++){
4154 treeT->SetLineColor(i+2);
4155 treeT->SetMarkerSize(1);
4156 treeT->SetMarkerStyle(22+i);
4157 treeT->SetMarkerColor(i+2);
4159 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4167 TTree * tree = (TTree*)f.Get("FitModels");
4169 TEventList listLFit0("listLFit0","listLFit0");
4170 TEventList listLFit1("listLFit1","listLFit1");
4171 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4172 tree->SetEventList(&listLFit0);
4177 gSystem->Load("libSTAT.so")
4178 TStatToolkit toolkit;
4184 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4185 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4189 fstring+="cos(dp)++";
4190 fstring+="sin(dp)++";
4191 fstring+="cos(dt)++";
4192 fstring+="sin(dt)++";
4194 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4207 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4208 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4209 AliXRDPROOFtoolkit tool;
4210 TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4211 chainTrack->Lookup();
4212 chainTrack->SetProof(kTRUE);
4214 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4216 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4218 chainFit->SetProof(kTRUE);
4219 chain->SetProof(kTRUE);
4223 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4224 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4225 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4226 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4228 TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4229 TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4230 TCut cutN("nclO>20&&nclI>20");
4231 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4235 TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4236 TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4237 TCut cutClX("abs(Cl.fX)>10");
4238 TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4239 TCut cutCl=cutClY+cutClZ+cutClX;
4242 // check edge effects
4243 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4245 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4247 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
4251 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4252 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")