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);
235 for (Int_t i=0; i<336; i++) {
237 fCounter[i]=0; //! counter of usage
238 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
239 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
243 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
250 fTracksEsdParam(336),
254 fDeltaZ(336), // array of histograms of delta z for each track
255 fDeltaP3(336), // array of histograms of delta z for each track
256 fDeltaP4(336), // array of histograms of P3 for each track
257 fDeltaPhi(336), // array of histograms of P4 for each track
258 fDeltaPhiP(336), // array of histograms of delta z for each track
259 fSignals(336), // array of dedx signals
262 fHisLaser(0), // N dim histogram of laser
263 fHisNclIn(0), //->Number of clusters inner
264 fHisNclOut(0), //->Number of clusters outer
265 fHisNclIO(0), //->Number of cluster inner outer
266 fHisLclIn(0), //->Level arm inner
267 fHisLclOut(0), //->Level arm outer
268 fHisLclIO(0), //->Number of cluster inner outer
269 fHisdEdx(0), //->dEdx histo
270 fHisdZfit(0), //->distance to the mirror after linear fit
273 fHisChi2YIn1(0), //->chi2 y inner - line
274 fHisChi2YOut1(0), //->chi2 y inner - line
275 fHisChi2YIn2(0), //->chi2 y inner - parabola
276 fHisChi2YOut2(0), //->chi2 y inner - parabola
277 fHisChi2YIO1(0), //->chi2 y IO - common
278 fHisChi2ZIn1(0), //->chi2 z inner - line
279 fHisChi2ZOut1(0), //->chi2 z inner - line
280 fHisChi2ZIn2(0), //->chi2 z inner - parabola
281 fHisChi2ZOut2(0), //->chi2 z inner - parabola
282 fHisChi2ZIO1(0), //->chi2 z IO - common
285 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
286 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
287 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
288 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
289 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
290 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
291 fHisPy2vP2In(0), //-> Curv P2inner - parabola
292 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
293 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
296 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
297 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
298 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
299 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
300 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
301 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
302 fHisPz2vP2In(0), //-> Curv P2inner - parabola
303 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
304 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
318 fFitAside(new TVectorD(5)), // drift fit - A side
319 fFitCside(new TVectorD(5)), // drift fit - C- side
320 fFitACside(new TVectorD(6)), // drift fit - AC- side
321 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
322 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
323 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
324 fNcuts(0), // number of cuts
325 fBeamSectorOuter(336),
326 fBeamSectorInner(336),
327 fBeamOffsetYOuter(336),
328 fBeamSlopeYOuter(336),
329 fBeamOffsetYInner(336),
330 fBeamSlopeYInner(336),
331 fBeamOffsetZOuter(336),
332 fBeamSlopeZOuter(336),
333 fBeamOffsetZInner(336),
334 fBeamSlopeZInner(336),
335 fInverseSlopeZ(kTRUE),
337 fFixedFitAside0(0.0),
338 fFixedFitAside1(1.0),
339 fFixedFitCside0(0.0),
347 fTracksEsdParam.SetOwner(kTRUE);
348 for (Int_t i=0; i<336; i++) {
350 fCounter[i]=0; //! counter of usage
351 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
352 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
356 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
357 AliTPCcalibBase(calibLaser),
363 fTracksEsdParam(336),
365 fFullCalib(calibLaser.fFullCalib),
367 fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
368 fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
369 fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
370 fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
371 fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
372 fSignals(((calibLaser.fSignals))), // array of dedx signals
375 fHisLaser(0), // N dim histogram of laser
376 fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
377 fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
378 fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
379 fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
380 fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
381 fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
382 fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
383 fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
386 fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
387 fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
388 fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
389 fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
390 fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
391 fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
392 fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
393 fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
394 fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
395 fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
398 fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
399 fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
400 fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
401 fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
402 fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
403 fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
404 fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
405 fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
406 fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
409 fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
410 fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
411 fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
412 fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
413 fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
414 fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
415 fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
416 fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
417 fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
420 fDeltaYres(((calibLaser.fDeltaYres))),
421 fDeltaZres(((calibLaser.fDeltaZres))),
422 fDeltaYres2(((calibLaser.fDeltaYres))),
423 fDeltaZres2(((calibLaser.fDeltaZres))),
424 fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
425 fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
426 fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
427 fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
428 // fDeltaYres3(((calibLaser.fDeltaYres))),
429 //fDeltaZres3(((calibLaser.fDeltaZres))),
430 fFitAside(new TVectorD(5)), // drift fit - A side
431 fFitCside(new TVectorD(5)), // drift fit - C- side
432 fFitACside(new TVectorD(6)), // drift fit - C- side
433 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
434 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
435 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
436 fNcuts(0), // number of cuts
437 fBeamSectorOuter(336),
438 fBeamSectorInner(336),
439 fBeamOffsetYOuter(336),
440 fBeamSlopeYOuter(336),
441 fBeamOffsetYInner(336),
442 fBeamSlopeYInner(336),
443 fBeamOffsetZOuter(336),
444 fBeamSlopeZOuter(336),
445 fBeamOffsetZInner(336),
446 fBeamSlopeZInner(336),
447 fInverseSlopeZ(calibLaser.fInverseSlopeZ),
448 fUseFixedDriftV(calibLaser.fUseFixedDriftV),
449 fFixedFitAside0(calibLaser.fFixedFitAside0),
450 fFixedFitAside1(calibLaser.fFixedFitAside1),
451 fFixedFitCside0(calibLaser.fFixedFitCside0),
452 fFixedFitCside1(calibLaser.fFixedFitCside1)
457 for (Int_t i=0; i<336; i++) {
459 fCounter[i]=0; //! counter of usage
460 fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
461 fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
467 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
469 // assgnment operator
471 if (this != &calibLaser) {
472 new (this) AliTPCcalibLaser(calibLaser);
481 AliTPCcalibLaser::~AliTPCcalibLaser() {
486 delete fHisLaser; //->
487 delete fHisNclIn; //->Number of clusters inner
488 delete fHisNclOut; //->Number of clusters outer
489 delete fHisNclIO; //->Number of cluster inner outer
490 delete fHisLclIn; //->Level arm inner
491 delete fHisLclOut; //->Level arm outer
492 delete fHisLclIO; //->Number of cluster inner outer
497 delete fHisChi2YIn1; //->chi2 y inner - line
498 delete fHisChi2YOut1; //->chi2 y inner - line
499 delete fHisChi2YIn2; //->chi2 y inner - parabola
500 delete fHisChi2YOut2; //->chi2 y inner - parabola
501 delete fHisChi2YIO1; //->chi2 y IO - common
502 delete fHisChi2ZIn1; //->chi2 z inner - line
503 delete fHisChi2ZOut1; //->chi2 z inner - line
504 delete fHisChi2ZIn2; //->chi2 z inner - parabola
505 delete fHisChi2ZOut2; //->chi2 z inner - parabola
506 delete fHisChi2ZIO1; //->chi2 z IO - common
509 delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
510 delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
511 delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
512 delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
513 delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
514 delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
515 delete fHisPy2vP2In; //-> Curv P2inner - parabola
516 delete fHisPy2vP2Out; //-> Curv P2outer - parabola
517 delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
520 delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
521 delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
522 delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
523 delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
524 delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
525 delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
526 delete fHisPz2vP2In; //-> Curv P2inner - parabola
527 delete fHisPz2vP2Out; //-> Curv P2outer - parabola
528 delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
533 fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
534 fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
535 fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
536 fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
537 fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
538 fSignals.SetOwner(); //->Array of dedx signals
540 fDeltaZ.Delete(); //-> array of histograms of delta z for each track
541 fDeltaP3.Delete(); //-> array of histograms of P3 for each track
542 fDeltaP4.Delete(); //-> array of histograms of P4 for each track
543 fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
544 fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
545 fSignals.Delete(); //->Array of dedx signals
547 fDeltaYres.SetOwner();
549 delete fHisYAbsErrors;
550 fDeltaZres.SetOwner();
552 delete fHisZAbsErrors;
553 fDeltaYres2.SetOwner();
554 fDeltaYres2.Delete();
555 fDeltaZres2.SetOwner();
556 fDeltaZres2.Delete();
558 fDeltaYresAbs.SetOwner();
559 fDeltaYresAbs.Delete();
560 fDeltaZresAbs.SetOwner();
561 fDeltaZresAbs.Delete();
566 void AliTPCcalibLaser::Process(AliESDEvent * event) {
569 // Loop over tracks and call Process function
571 const Int_t kMinTracks=20;
572 const Int_t kMinClusters=40;
578 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
582 if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
583 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
585 // find CE background if present
587 if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
588 TH1D hisCE("hhisCE","hhisCE",100,-100,100);
589 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
590 AliESDtrack *track=fESD->GetTrack(i);
591 if (!track) continue;
592 hisCE.Fill(track->GetZ());
593 hisCE.Fill(track->GetZ()+2);
594 hisCE.Fill(track->GetZ()-2);
602 fTracksEsdParam.Delete();
603 for (Int_t id=0; id<336;id++) {
605 fClusterCounter[id]=0;
609 Int_t n=fESD->GetNumberOfTracks();
611 for (Int_t i=0;i<n;++i) {
612 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
613 if (!friendTrack) continue;
614 AliESDtrack *track=fESD->GetTrack(i);
615 if (!track) continue;
616 Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
617 if (binC>336) continue; //remove CE background
618 TObject *calibObject=0;
620 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
621 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
623 if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
625 Int_t id = FindMirror(track,seed);
626 if (id>=0) counter++;
631 if (counter<kMinTracks) return;
635 if (!fFullCalib) return;
636 static Bool_t init=kFALSE;
638 init = kTRUE; // way around for PROOF - to be investigated
642 for (Int_t id=0; id<336; id++){
644 if (!fTracksEsdParam.At(id)) continue;
645 if (fClusterSatur[id]>0.3) continue; // tracks in saturation
647 if ( AcceptLaser(id) && fFitZ[id]<0.5){
655 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
659 // for (Int_t id=0; id<336; id++){
662 if (!fTracksEsdParam.At(id)) return;
663 if (!AcceptLaser(id)) return;
664 Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
667 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
668 if (!hisdz) MakeFitHistos();
669 hisdz = (TH1F*)fDeltaZ.At(id);
670 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
671 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
672 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
673 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
674 TH1F * hisSignal = (TH1F*)fSignals.At(id);
677 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
678 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
679 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
688 param->GetPxPyPz(pxyz);
690 ltrp->GetPxPyPz(lpxyz);
692 Float_t dz = param->GetZ()-ltrp->GetZ();
693 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
694 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
695 if (hisdz) hisdz->Fill(dz);
696 if (hisP3) hisP3->Fill(param->GetParameter()[3]);
697 if (hisP4) hisP4->Fill(param->GetParameter()[4]);
698 if (hisdphi) hisdphi->Fill(dphi);
699 if (hisdphiP) hisdphiP->Fill(dphiP);
700 if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
702 xhis[0] = ltrp->GetId();
703 xhis[1] = ltrp->GetSide();
704 xhis[2] = ltrp->GetRod();
705 xhis[3] = ltrp->GetBundle();
706 xhis[4] = ltrp->GetBeam();
709 xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
710 xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
711 xhis[9] = param->GetParameter()[4];
712 xhis[10]= track->GetTPCNcls();
713 xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
715 fHisLaser->Fill(xhis);
718 void AliTPCcalibLaser::FitDriftV(){
720 // Fit corrections to the drift velocity - linear approximation in the z and global y
721 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
726 z = s* (z0 - vd*(t-t0))
730 vd - nominal drift velocity
731 zs - miscalibrated position
733 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
734 vr - relative change of the drift velocity
739 zs ~ z - s*vr*(z0-s*z)+s*dzt
740 --------------------------------
741 1. Correction function vr constant:
744 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
745 dzs/dl = dz/dl +s*s*vr*dz/dl
749 const Float_t kZCut = 240; // remove the closest laser beam
750 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
751 const Float_t kDistCut = 3; // distance sigma cut
752 const Float_t kDistCutAbs = 0.25;
753 const Float_t kMinClusters = 60; // minimal amount of the clusters
754 const Float_t kMinSignal = 16; // minimal mean height of the signal
755 const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
756 static TLinearFitter fdriftA(3,"hyp2");
757 static TLinearFitter fdriftC(3,"hyp2");
758 static TLinearFitter fdriftAC(4,"hyp3");
759 TVectorD fitA(3),fitC(3),fitAC(4);
761 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
762 AliTPCParam * tpcparam = calib->GetParameters();
764 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
774 for (Int_t iter=0; iter<3; iter++){
775 fdriftA.ClearPoints();
776 fdriftC.ClearPoints();
777 fdriftAC.ClearPoints();
779 for (Int_t id=0; id<336; id++){
780 if (!fTracksEsdParam.At(id)) continue;
781 if (!AcceptLaser(id)) continue;
782 if ( fClusterSatur[id]>kSaturCut) continue;
783 if ( fClusterCounter[id]<kMinClusters) continue;
784 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
785 if (track->GetTPCsignal()<kMinSignal) continue;
786 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
787 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
794 param->GetPxPyPz(pxyz);
796 ltrp->GetPxPyPz(lpxyz);
797 if (TMath::Abs(lxyz[2])>kZCut) continue;
798 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
799 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
800 if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
801 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
804 Double_t zlength = tpcparam->GetZLength(0);
805 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
806 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
807 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
808 if (ltrp->GetSide()==0){
809 fdriftA.AddPoint(xxx,mdrift,1);
811 fdriftC.AddPoint(xxx,mdrift,1);
813 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
814 fdriftAC.AddPoint(xxx2,mdrift,1);
817 if (fdriftA.GetNpoints()>10){
820 if (iter==0) fdriftA.EvalRobust(0.7);
821 else fdriftA.EvalRobust(0.8);
822 fdriftA.GetParameters(fitA);
823 npointsA= fdriftA.GetNpoints();
824 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
825 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
826 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
827 (*fFitAside)[0] = fitA[0];
828 (*fFitAside)[1] = fitA[1];
829 (*fFitAside)[2] = fitA[2];
830 (*fFitAside)[3] = fdriftA.GetNpoints();
831 (*fFitAside)[4] = chi2A;
834 if (fdriftC.GetNpoints()>10){
836 if (iter==0) fdriftC.EvalRobust(0.7);
837 else fdriftC.EvalRobust(0.8);
839 fdriftC.GetParameters(fitC);
840 npointsC= fdriftC.GetNpoints();
841 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
842 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
843 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
844 (*fFitCside)[0] = fitC[0];
845 (*fFitCside)[1] = fitC[1];
846 (*fFitCside)[2] = fitC[2];
847 (*fFitCside)[3] = fdriftC.GetNpoints();
848 (*fFitCside)[4] = chi2C;
852 if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
854 if (iter==0) fdriftAC.EvalRobust(0.7);
855 else fdriftAC.EvalRobust(0.8);
857 fdriftAC.GetParameters(fitAC);
858 npointsAC= fdriftAC.GetNpoints();
859 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
860 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
863 for (Int_t id=0; id<336; id++){
864 if (!fTracksEsdParam.At(id)) continue;
866 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
867 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
873 param->GetPxPyPz(pxyz);
875 ltrp->GetPxPyPz(lpxyz);
876 Double_t zlength = tpcparam->GetZLength(0);
877 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
878 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
881 if (ltrp->GetSide()==0){
882 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
884 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
887 fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
892 TTreeSRedirector *cstream = GetDebugStreamer();
893 TTimeStamp tstamp(fTime);
894 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
895 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
896 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
897 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
898 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
899 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
900 TVectorD vecGoofie(20);
901 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
903 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
904 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
905 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
909 (*cstream)<<"driftv"<<
910 "run="<<fRun<< // run number
911 "event="<<fEvent<< // event number
912 "time="<<fTime<< // time stamp of event
913 "trigger="<<fTrigger<< // trigger
914 "mag="<<fMagF<< // magnetic field
915 // Environment values
916 "press0="<<valuePressure0<<
917 "press1="<<valuePressure1<<
918 "pt0="<<ptrelative0<<
919 "pt1="<<ptrelative1<<
922 "vecGoofie.="<<&vecGoofie<<
926 "driftA.="<<fFitAside<<
927 "driftC.="<<fFitCside<<
928 "driftAC.="<<fFitACside<<
940 Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
942 // Fit corrections to the drift velocity - linear approximation in the z and global y
943 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
945 // Source of outlyers :
946 // 0. Track in the saturation - postpeak
947 // 1. gating grid close the part of the signal for first bundle
949 // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
950 // 1. Robust fit is used in the itteration number 0
951 // only fraction of laser used
952 // 2. Only the tracks close to the fit used in the second itteration
956 z = s* (z0 - vd*(t-t0))
960 vd - nominal drift velocity
961 zs - miscalibrated position
963 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
964 vr - relative change of the drift velocity
969 zs ~ z - s*vr*(z0-s*z)+s*dzt
970 --------------------------------
971 1. Correction function vr constant:
974 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
975 dzs/dl = dz/dl +s*s*vr*dz/dl
978 const Int_t knLaser = 336; //n laser tracks
979 const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
981 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
982 const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
983 const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
984 const Float_t kMinClusters = 40.; // minimal amount of the clusters
985 const Float_t kMinSignal = 2.5; // minimal mean height of the signal
986 const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
988 static TLinearFitter fdriftA(3,"hyp2");
989 static TLinearFitter fdriftC(3,"hyp2");
990 static TLinearFitter fdriftAC(4,"hyp3");
991 TVectorD fitA(3),fitC(3),fitAC(4);
993 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
994 AliTPCParam * tpcparam = calib->GetParameters();
998 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
999 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1000 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1001 for (Int_t i=0;i<5; i++){
1009 Float_t chi2AC = 10;
1013 Int_t nbA[4]={0,0,0,0};
1014 Int_t nbC[4]={0,0,0,0};
1015 TVectorD vecZM(336); // measured z potion of laser
1016 TVectorD vecA(336); // accepted laser
1017 TVectorD vecZF(336); // fitted position
1018 TVectorD vecDz(336); // deltaZ
1019 TVectorD vecZS(336); // surveyed position of laser
1020 // additional variable to cut
1021 TVectorD vecdEdx(336); // dEdx
1022 TVectorD vecSy(336); // shape y
1023 TVectorD vecSz(336); // shape z
1026 for (Int_t id=0; id<336; id++){
1028 AliTPCLaserTrack *ltrp =
1029 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1030 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1031 AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1032 vecZM(id)= (param==0) ? 0:param->GetZ();
1033 vecZS(id)= ltrp->GetZ();
1036 vecdEdx(id)=(seed)?seed->GetdEdx():0;
1037 vecSy(id) =(seed)?seed->CookShape(1):0;
1038 vecSz(id) =(seed)?seed->CookShape(2):0;
1041 if (!fTracksEsdParam.At(id)) reject|=1;
1042 if (!param) continue;
1043 if (!AcceptLaser(id)) reject|=2;
1044 if ( fClusterSatur[id]>kSaturCut) reject|=4;
1045 if ( fClusterCounter[id]<kMinClusters) reject|=8;
1047 if (reject>0) continue;
1048 if (ltrp->GetSide()==0){
1050 nbA[ltrp->GetBundle()]++;
1052 if (ltrp->GetSide()==1){
1054 nbC[ltrp->GetBundle()]++;
1058 // reject "bad events"
1063 if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1065 for (Int_t i=0;i<4;i++){
1066 //count accepted for all layers
1067 if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1068 if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1070 if (naA<3 &&naC<3) isOK=kFALSE;
1071 if (isOK==kFALSE) return kFALSE;
1076 for (Int_t iter=0; iter<2; iter++){
1077 fdriftA.ClearPoints();
1078 fdriftC.ClearPoints();
1079 fdriftAC.ClearPoints();
1080 npointsA=0; npointsC=0; npointsAC=0;
1082 for (Int_t id=0; id<336; id++){
1083 if (!fTracksEsdParam.At(id)) continue;
1084 if (!AcceptLaser(id)) continue;
1085 if ( fClusterSatur[id]>kSaturCut) continue;
1086 if ( fClusterCounter[id]<kMinClusters) continue;
1087 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1088 if (track->GetTPCsignal()<kMinSignal) continue;
1089 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1090 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1096 param->GetPxPyPz(pxyz);
1098 ltrp->GetPxPyPz(lpxyz);
1099 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1100 //if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1101 if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1102 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1104 // // drift distance
1105 // Double_t zlength = tpcparam->GetZLength(0);
1106 // Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1107 // Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1109 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1110 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1111 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1114 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1115 if (iter==0 &<rp->GetBundle()==0) continue;
1116 // skip bundle 0 - can be wrong in the 0 iteration
1117 if (ltrp->GetSide()==0){
1118 fdriftA.AddPoint(xxx,mdrift,1);
1120 fdriftC.AddPoint(xxx,mdrift,1);
1122 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1123 fdriftAC.AddPoint(xxx2,mdrift,1);
1126 if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1129 //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
1130 npointsA= fdriftA.GetNpoints();
1131 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1132 fdriftA.EvalRobust(kFraction[iter]);
1133 fdriftA.GetParameters(fitA);
1134 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1135 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1136 (*fFitAside)[0] = fitA[0];
1137 (*fFitAside)[1] = fitA[1];
1138 (*fFitAside)[2] = fitA[2];
1139 (*fFitAside)[3] = fdriftA.GetNpoints();
1140 (*fFitAside)[4] = chi2A;
1143 if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1145 //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
1146 npointsC= fdriftC.GetNpoints();
1147 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1148 fdriftC.EvalRobust(kFraction[iter]);
1149 fdriftC.GetParameters(fitC);
1150 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1151 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1152 (*fFitCside)[0] = fitC[0];
1153 (*fFitCside)[1] = fitC[1];
1154 (*fFitCside)[2] = fitC[2];
1155 (*fFitCside)[3] = fdriftC.GetNpoints();
1156 (*fFitCside)[4] = chi2C;
1160 if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1162 //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
1163 npointsAC= fdriftAC.GetNpoints();
1164 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1165 fdriftAC.EvalRobust(kFraction[iter]);
1166 fdriftAC.GetParameters(fitAC);
1167 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1168 (*fFitACside)[0] = fitAC[0];
1169 (*fFitACside)[1] = fitAC[1];
1170 (*fFitACside)[2] = fitAC[2];
1171 (*fFitACside)[3] = fdriftAC.GetNpoints();
1172 (*fFitACside)[4] = chi2AC;
1175 for (Int_t id=0; id<336; id++){
1176 if (!fTracksEsdParam.At(id)) continue;
1178 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1179 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1185 param->GetPxPyPz(pxyz);
1187 ltrp->GetPxPyPz(lpxyz);
1188 //Double_t zlength = tpcparam->GetZLength(0);
1189 //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1190 //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1191 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1192 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1193 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1197 if (ltrp->GetSide()==0){
1198 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1200 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1203 //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1205 fFitZ[id]=mdrift-fz;
1207 vecDz[id]=mdrift-fz;
1209 if (fStreamLevel>0){
1210 TTreeSRedirector *cstream = GetDebugStreamer();
1211 TTimeStamp tstamp(fTime);
1212 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1213 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1214 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1215 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1216 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1217 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1218 TVectorD vecGoofie(20);
1219 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1221 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1222 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1223 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1227 (*cstream)<<"driftvN"<<
1228 "run="<<fRun<< // run number
1229 "event="<<fEvent<< // event number
1230 "time="<<fTime<< // time stamp of event
1231 "trigger="<<fTrigger<< // trigger
1232 "mag="<<fMagF<< // magnetic field
1233 // Environment values
1234 "press0="<<valuePressure0<<
1235 "press1="<<valuePressure1<<
1236 "pt0="<<ptrelative0<<
1237 "pt1="<<ptrelative1<<
1240 "vecGoofie.="<<&vecGoofie<<
1243 "vecZM.="<<&vecZM<< // measured z position
1244 "vecZS.="<<&vecZS<< // surveyed z position
1245 "vecZF.="<<&vecZF<< // fitted z position
1246 "vecDz.="<<&vecDz<< // fitted z position
1247 "vecA.="<<&vecA<< // accept laser flag
1248 "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
1249 "vecSy.="<<&vecSy<< // shape y - to cut on
1250 "vecSz.="<<&vecSz<< // shape z - to cut on
1253 "driftA.="<<fFitAside<<
1254 "driftC.="<<fFitCside<<
1255 "driftAC.="<<fFitACside<<
1265 variables to check in debug mode:
1267 chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1268 chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1269 chainDriftN->SetAlias("driftF","vecZF.fElements");
1270 chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
1271 TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1272 TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1281 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1283 // get distance between mirror and laser track
1294 dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1297 // apply drift correction if already exist
1300 if (ltrp->GetSide()==0){
1301 if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1303 if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1305 if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1306 dist+=TMath::Abs(dz);
1308 // phi dist - divergence on 50 cm
1310 dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1315 Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
1320 TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1321 TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1322 TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1323 TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1325 TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1327 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1328 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1329 //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1334 if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1335 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
1336 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
1337 if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
1344 Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1346 // Find corresponding mirror
1347 // add the corresponding tracks
1349 if (!track->GetOuterParam()) return -1;
1351 Float_t kRadius0 = 252;
1352 Float_t kRadius = 254.2;
1354 Float_t counterSatur=0;
1357 for (Int_t irow=158;irow>-1;--irow) {
1358 AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1360 Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1361 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1362 if (pedgeY<3) continue;
1363 if (pedgeX<3) continue;
1365 if (c->GetDetector()%36<18) csideA++;
1366 if (c->GetDetector()%36>=18) csideC++;
1367 if (c->GetMax()>900) counterSatur++;
1369 counterSatur/=(countercl+1);
1373 if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
1376 if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1379 AliExternalTrackParam param(*(track->GetOuterParam()));
1380 AliTracker::PropagateTrackTo(¶m,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
1381 AliTracker::PropagateTrackTo(¶m,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
1382 AliTPCLaserTrack ltr;
1383 AliTPCLaserTrack *ltrp=0x0;
1384 // AliTPCLaserTrack *ltrpjw=0x0;
1386 Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
1387 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
1388 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1390 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1391 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1395 if (id<0) return -1;
1396 if (ltrp->GetSide()!=side) return -1;
1402 if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1405 Float_t radius=TMath::Abs(ltrp->GetX());
1406 param.Rotate(ltrp->GetAlpha());
1407 AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
1409 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1410 Bool_t accept=kTRUE;
1412 // choose closer track
1414 AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1416 Float_t dist0=GetDistance(param0,ltrp);
1417 Float_t dist1=GetDistance(¶m,ltrp);
1418 if (dist0<dist1) accept=kFALSE;
1422 fClusterCounter[id]=countercl;
1423 fTracksEsdParam.AddAt(param.Clone(),id);
1424 fTracksEsd.AddAt(track,id);
1425 fTracksTPC.AddAt(seed,id);
1432 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1434 // Dump Laser info to the tree
1436 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1437 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1438 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1447 param->GetPxPyPz(pxyz);
1449 ltrp->GetPxPyPz(lpxyz);
1450 Float_t dist3D = GetDistance(param,ltrp);
1451 Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1452 Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1455 if (fStreamLevel>0){
1456 TTreeSRedirector *cstream = GetDebugStreamer();
1457 Int_t time = fESD->GetTimeStamp();
1458 Bool_t accept = AcceptLaser(id);
1460 (*cstream)<<"Track"<<
1462 "run="<<fRun<< // run number
1463 "event="<<fEvent<< // event number
1464 "time="<<fTime<< // time stamp of event
1465 "trigger="<<fTrigger<< // trigger
1466 "mag="<<fMagF<< // magnetic field
1470 "driftA.="<<fFitAside<<
1471 "driftC.="<<fFitCside<<
1475 "distphi="<<distphi<<
1478 "counter="<<fCounter[id]<<
1479 "clcounter="<<fClusterCounter[id]<<
1480 "satur="<<fClusterSatur[id]<<
1504 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1506 // Refit the track with different tracklet models:
1507 // 1. Per ROC using the kalman filter, different edge cuts
1508 // 2. Per ROC linear in y and z
1509 // 3. Per ROC quadratic in y and z
1510 // 4. Per track offset for each sector, linear for each sector, common quadratic
1511 // store x, y, z information for all models and the cluster to calculate the residuals
1514 // The clusters which do not fulfill given criteria are skipped
1516 // Cluters removed from fit
1517 // 1. Extended shape > kShapeCut
1518 // 2. In saturation Max > kMax
1519 // 3. Distance to edge < cutEdge
1521 // Clusters not used for the calibration:
1523 // 1. Extended shape > kShapeCut
1524 // 2. In saturation Max > kMax
1527 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1528 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1529 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1531 AliTPCclusterMI dummyCl;
1534 Int_t kMaxTracklets=2;
1535 Float_t kShapeCut = 1.3;
1536 Float_t kRatioCut = 2.;
1538 Float_t kMax = 900.;
1541 //=============================================//
1542 // Linear Fitters for the Different Approaches //
1543 //=============================================//
1544 //linear fit model in y and z; inner - outer sector, combined with offset
1545 static TLinearFitter fy1I(2,"hyp1");
1546 static TLinearFitter fy1O(2,"hyp1");
1547 static TLinearFitter fz1I(2,"hyp1");
1548 static TLinearFitter fz1O(2,"hyp1");
1549 static TLinearFitter fy1IO(3,"hyp2");
1550 static TLinearFitter fz1IO(3,"hyp2");
1551 //quadratic fit model in y and z; inner - sector
1552 static TLinearFitter fy2I(3,"hyp2");
1553 static TLinearFitter fy2O(3,"hyp2");
1554 static TLinearFitter fz2I(3,"hyp2");
1555 static TLinearFitter fz2O(3,"hyp2");
1556 //common quadratic fit for IROC and OROC in y and z
1557 static TLinearFitter fy4(5,"hyp4");
1558 static TLinearFitter fz4(5,"hyp4");
1568 //=============================//
1569 // Loop over all Tracklet Cuts //
1570 //=============================//
1571 for (Int_t icut=0; icut<fNcuts; icut++){
1572 Float_t xinMin = 2500, xinMax=-90;
1573 Float_t xoutMin=2500, xoutMax=-90;
1574 Float_t msigmaYIn=0;
1575 Float_t msigmaYOut=0;
1576 Float_t mqratioIn=0;
1577 Float_t mqratioOut=0;
1580 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1582 Double_t edgeCutX = fEdgeXcuts[icut];
1583 Double_t edgeCutY = fEdgeYcuts[icut];
1584 Int_t nclCut = (Int_t)fNClCuts[icut];
1585 //=========================//
1586 // Parameters to calculate //
1587 //=========================//
1588 //fit parameter inner, outer tracklet and combined fit
1589 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1590 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1592 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1593 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1594 TVectorD vecy4res(5),vecz4res(5);
1595 TVectorD vecy1resIO(3),vecz1resIO(3);
1596 // cluster and track positions for each row - used for residuals
1597 TVectorD vecgX(159); // global X
1598 TVectorD vecgY(159); // global Y
1599 TVectorD vecgZ(159); // global Z
1601 TVectorD vecX(159); // x is the same for all (row center)
1602 TVectorD vecYkalman(159); // y from kalman fit
1603 TVectorD vecZkalman(159); // z from kalman fit
1604 TVectorD vecY1(159); // y from pol1 fit per ROC
1605 TVectorD vecZ1(159); // z from pol1 fit per ROC
1606 TVectorD vecY1IO(159); // y from pol1 fit per ROC
1607 TVectorD vecZ1IO(159); // z from pol1 fit per ROC
1608 TVectorD vecY2(159); // y from pol2 fit per ROC
1609 TVectorD vecZ2(159); // z from pol2 fit per ROC
1610 TVectorD vecY4(159); // y from sector fit
1611 TVectorD vecZ4(159); // z from sector fit
1612 TVectorD vecClY(159); // y cluster position
1613 TVectorD vecClZ(159); // z cluster position
1614 TVectorD vecSec(159); // sector for each row
1615 TVectorD isReject(159); // flag - cluster to be rejected
1617 Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
1618 Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
1619 Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
1620 Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
1621 Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
1622 Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
1623 Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
1624 Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
1625 Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
1626 Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
1627 Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
1628 Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
1630 Int_t innerSector = -1; // number of inner sector
1631 Int_t outerSector = -1; // number of outer sector
1632 Int_t nclI=0; // number of clusters (inner)
1633 Int_t nclO=0; // number of clusters (outer)
1634 //=================================================//
1635 // Perform the Kalman Fit using the Tracklet Class //
1636 //=================================================//
1637 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1638 TObjArray tracklets=
1639 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1640 kFALSE,nclCut,kMaxTracklets);
1641 // tracklet pointers
1642 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1643 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1644 AliTPCTracklet *tr=0x0;
1645 AliTPCTracklet dummy;
1646 //continue if we didn't find a tracklet
1647 if ( !trInner && !trOuter ) continue;
1648 //================================================================================//
1649 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1650 //================================================================================//
1651 if ( trInner && trOuter ){
1652 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1653 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1657 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1661 if ( !trInner->GetInner() ) continue;
1663 if ( trInner->GetSector()>35 ){
1668 if ( !trOuter->GetInner() ) continue;
1670 if ( trOuter->GetSector()<36 ){
1676 innerSector = trInner->GetSector();
1677 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1678 outerSector = trOuter->GetSector();
1679 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1681 // array of clusters
1682 TClonesArray arrCl("AliTPCclusterMI",159);
1683 arrCl.ExpandCreateFast(159);
1684 //=======================================//
1685 // fill fitters with cluster information //
1686 //=======================================//
1687 AliDebug(3,"Filing Cluster Information");
1692 for (Int_t irow=158;irow>-1;--irow) {
1693 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1694 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1699 Float_t meanY=0, sumY=0;
1700 for (Int_t drow=-1;drow<=1;drow++) {
1701 if (irow+drow<0) continue;
1702 if (irow+drow>158) continue;
1703 AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1704 if (!ccurrent) continue;
1705 Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1706 if ( roc!=innerSector && roc!=outerSector ) continue;
1707 if (ccurrent->GetX()<10) continue;
1708 if (ccurrent->GetY()==0) continue;
1709 meanY+=ccurrent->GetY();
1712 if (sumY>0) meanY/=sumY;
1717 Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1718 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1721 Int_t roc = static_cast<Int_t>(c->GetDetector());
1722 if ( roc!=innerSector && roc!=outerSector ) continue;
1724 //store clusters in clones array
1727 if (c->GetMax()<4) continue; // noise cluster?
1728 if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
1730 vecX[irow] = c->GetX();
1731 vecClY[irow] = c->GetY();
1732 vecClZ[irow] = c->GetZ();
1735 // c->GetGlobalXYZ(gxyz);
1736 // vecgX[irow] = gxyz[0];
1737 // vecgY[irow] = gxyz[1];
1738 // vecgZ[irow] = gxyz[2];
1740 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1741 Double_t y = vecClY[irow];
1742 Double_t z = vecClZ[irow];
1744 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
1745 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1746 Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
1747 if ( roc == innerSector ) {
1748 x4[0]=1; //offset inner - outer sector
1749 x4[1]=x; //slope parameter inner sector
1752 x4[2]=x; //slope parameter outer sector
1754 x4[3]=x*x; //common parabolic shape
1755 if (pedgeX < fEdgeXcuts[icut]) continue;
1756 if (pedgeY < fEdgeYcuts[icut]) continue;
1757 if (c->GetMax()>900) continue; // cluster in saturation
1759 if ( roc==innerSector ){
1760 fy1I.AddPoint(x2,y);
1761 fz1I.AddPoint(x2,z);
1762 fy2I.AddPoint(x2,y);
1763 fz2I.AddPoint(x2,z);
1765 if (c->GetX()<xinMin) xinMin=c->GetX();
1766 if (c->GetX()>xinMax) xinMax=c->GetX();
1768 msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1769 mqratioIn +=c->GetMax()/c->GetQ();
1772 if ( roc==outerSector ){
1773 fy1O.AddPoint(x2,y);
1774 fz1O.AddPoint(x2,z);
1775 fy2O.AddPoint(x2,y);
1776 fz2O.AddPoint(x2,z);
1778 if (c->GetX()<xoutMin) xoutMin=c->GetX();
1779 if (c->GetX()>xoutMax) xoutMax=c->GetX();
1780 msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1781 mqratioOut +=c->GetMax()/c->GetQ();
1786 fy1IO.AddPoint(xIO,y);
1787 fz1IO.AddPoint(xIO,z);
1797 //======================================//
1798 // Evaluate and retrieve fit parameters //
1799 //======================================//
1800 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1802 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1807 fy1I.GetParameters(vecy1resInner);
1808 fz1I.GetParameters(vecz1resInner);
1809 fy2I.GetParameters(vecy2resInner);
1810 fz2I.GetParameters(vecz2resInner);
1811 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1812 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1813 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1814 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1817 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1822 fy1O.GetParameters(vecy1resOuter);
1823 fz1O.GetParameters(vecz1resOuter);
1824 fy2O.GetParameters(vecy2resOuter);
1825 fz2O.GetParameters(vecz2resOuter);
1826 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1827 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1828 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1829 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1832 if ( innerSector>0 && outerSector>0 ){
1833 if (fy4.GetNpoints()>0) {
1835 fy4.GetParameters(vecy4res);
1836 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1838 if (fz4.GetNpoints()>0) {
1840 fz4.GetParameters(vecz4res);
1841 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1843 if (fy1IO.GetNpoints()>0) {
1845 fy1IO.GetParameters(vecy1resIO);
1846 chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1848 if (fz1IO.GetNpoints()>0) {
1850 fz1IO.GetParameters(vecz1resIO);
1851 chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1855 fy4.ClearPoints(); fz4.ClearPoints();
1856 fy1I.ClearPoints(); fy1O.ClearPoints();
1857 fz1I.ClearPoints(); fz1O.ClearPoints();
1858 fy2I.ClearPoints(); fy2O.ClearPoints();
1859 fz2I.ClearPoints(); fz2O.ClearPoints();
1860 fy1IO.ClearPoints(); fz1IO.ClearPoints();
1861 //==============================//
1862 // calculate tracklet positions //
1863 //==============================//
1864 AliDebug(4,"Calculate tracklet positions");
1865 for (Int_t irow=158;irow>-1;--irow) {
1867 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1868 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1872 if (!c) { //no cluster
1875 if (c->GetMax()>kMax){ //saturation
1878 if ( vecSec[irow] == outerSector ) { //extended shape
1879 if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
1880 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
1882 if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
1883 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
1889 if ( vecSec[irow]==-1 ) continue; //no cluster info
1890 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1892 Double_t x=vecX[irow];
1893 Double_t xref=x-133.4;
1895 Double_t yoffInner=0;
1896 Double_t zoffInner=0;
1897 Double_t yoffInner1=0;
1898 Double_t zoffInner1=0;
1899 Double_t yslopeInner=0;
1900 Double_t yslopeOuter=0;
1901 Double_t zslopeInner=0;
1902 Double_t zslopeOuter=0;
1903 //positions of hyperplane fits
1904 if ( vecSec[irow] == outerSector ) {
1906 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1907 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1908 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1909 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1910 yslopeOuter=vecy4res[3];
1911 zslopeOuter=vecz4res[3];
1914 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1915 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1916 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1917 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1918 yoffInner=vecy4res[1];
1919 zoffInner=vecz4res[1];
1920 yoffInner1=vecy1resIO[1];
1921 zoffInner1=vecz1resIO[1];
1922 yslopeInner=vecy4res[2];
1923 zslopeInner=vecz4res[2];
1925 vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1926 vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1927 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1928 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1929 //positions of kalman fits
1930 Double_t gxyz[3],xyz[3];
1931 AliExternalTrackParam *param = 0x0;
1933 param=tr->GetInner();
1935 param->GetXYZ(gxyz);
1936 Float_t bz = AliTracker::GetBz(gxyz);
1937 param->GetYAt(x, bz, xyz[1]);
1938 param->GetZAt(x, bz, xyz[2]);
1939 vecYkalman[irow]=xyz[1];
1940 vecZkalman[irow]=xyz[2];
1947 //=====================================================================//
1948 // write results from the different tracklet fits with debug streamers //
1949 //=====================================================================//
1950 if (fStreamLevel>4){
1951 TTreeSRedirector *cstream = GetDebugStreamer();
1953 Float_t dedx = track->GetdEdx();
1954 (*cstream)<<"FitModels"<<
1955 "run="<<fRun<< // run number
1956 "event="<<fEvent<< // event number
1957 "time="<<fTime<< // time stamp of event
1958 "trigger="<<fTrigger<< // trigger
1959 "mag="<<fMagF<< // magnetic field
1962 "edgeCutX=" << edgeCutX <<
1963 "edgeCutY=" << edgeCutY <<
1964 "nclCut=" << nclCut <<
1965 "innerSector="<< innerSector <<
1966 "outerSector="<< outerSector <<
1969 "Tr.=" << extparam <<
1970 "yPol1In.=" << &vecy1resInner <<
1971 "zPol1In.=" << &vecz1resInner <<
1972 "yPol1InOut.="<< &vecy1resIO <<
1973 "zPol1InOut.="<< &vecz1resIO <<
1974 "yPol2In.=" << &vecy2resInner <<
1975 "zPol2In.=" << &vecz2resInner <<
1976 "yPol1Out.=" << &vecy1resOuter <<
1977 "zPol1Out.=" << &vecz1resOuter <<
1978 "yPol2Out.=" << &vecy2resOuter <<
1979 "zPol2Out.=" << &vecz2resOuter <<
1980 "yInOut.=" << &vecy4res <<
1981 "zInOut.=" << &vecz4res <<
1982 "chi2y1In=" << chi2I1y <<
1983 "chi2z1In=" << chi2I1z <<
1984 "chi2y1InOut="<< chi2IO1y <<
1985 "chi2z1InOut="<< chi2IO1z <<
1986 "chi2y1Out=" << chi2O1y <<
1987 "chi2z1Out=" << chi2O1z <<
1988 "chi2y2In=" << chi2I2y <<
1989 "chi2z2In=" << chi2I2z <<
1990 "chi2y2Out=" << chi2O2y <<
1991 "chi2z2Out=" << chi2O2z <<
1992 "chi2yInOut=" << chi2IOy <<
1993 "chi2zInOut=" << chi2IOz <<
1994 "trletIn.=" << trInner <<
1995 "trletOut.=" << trOuter <<
1998 "xinMin=" << xinMin<<
1999 "xinMax=" << xinMax<<
2000 "xoutMin=" << xoutMin<<
2001 "xoutMax=" << xoutMax<<
2002 "msigmaYIn=" <<msigmaYIn<<
2003 "msigmaYOut=" <<msigmaYOut<<
2004 "mqratioIn=" <<mqratioIn<<
2005 "mqratioOut=" << mqratioOut <<
2010 // wirte residuals information
2011 if (fStreamLevel>5){
2012 TTreeSRedirector *cstream = GetDebugStreamer();
2014 Float_t dedx = track->GetdEdx();
2015 (*cstream)<<"Residuals"<<
2016 "run="<<fRun<< // run number
2017 "event="<<fEvent<< // event number
2018 "time="<<fTime<< // time stamp of event
2019 "trigger="<<fTrigger<< // trigger
2020 "mag="<<fMagF<< // magnetic field
2023 "edgeCutX=" << edgeCutX <<
2024 "edgeCutY=" << edgeCutY <<
2025 "nclCut=" << nclCut <<
2027 "Tr.=" << extparam<<
2030 "vX.=" << &vecgX<< // global x
2031 "vY.=" << &vecgY<< // global y
2032 "vZ.=" << &vecgZ<< // global z
2034 "TrYpol1.=" << &vecY1 <<
2035 "TrZpol1.=" << &vecZ1 <<
2036 "TrYpol2.=" << &vecY2 <<
2037 "TrZpol2.=" << &vecZ2 <<
2038 "TrYpol1InOut.="<< &vecY1IO <<
2039 "TrZpol1InOut.="<< &vecZ1IO <<
2040 "TrYInOut.=" << &vecY4 <<
2041 "TrZInOut.=" << &vecZ4 <<
2042 "ClY.=" << &vecClY <<
2043 "ClZ.=" << &vecClZ <<
2044 "isReject.=" << &isReject<<
2045 "sec.=" << &vecSec <<
2048 "xinMin=" << xinMin<<
2049 "xinMax=" << xinMax<<
2050 "xoutMin=" << xoutMin<<
2051 "xoutMax=" << xoutMax<<
2052 "msigmaYIn=" <<msigmaYIn<<
2053 "msigmaYOut=" <<msigmaYOut<<
2054 "mqratioIn=" <<mqratioIn<<
2055 "mqratioOut=" << mqratioOut <<
2056 "yInOut.=" << &vecy4res <<
2057 "zInOut.=" << &vecz4res <<
2059 "chi2y1In=" << chi2I1y << //
2060 "chi2z1In=" << chi2I1z <<
2061 "chi2y1Out=" << chi2O1y <<
2062 "chi2z1Out=" << chi2O1z <<
2063 "chi2y1InOut="<< chi2IO1y <<
2064 "chi2z1InOut="<< chi2IO1z <<
2065 "chi2y2In=" << chi2I2y <<
2066 "chi2z2In=" << chi2I2z <<
2067 "chi2y2Out=" << chi2O2y <<
2068 "chi2z2Out=" << chi2O2z <<
2069 "chi2yInOut=" << chi2IOy <<
2070 "chi2zInOut=" << chi2IOz <<
2072 "yPol1In.=" << &vecy1resInner <<
2073 "zPol1In.=" << &vecz1resInner <<
2074 "yPol2In.=" << &vecy2resInner <<
2075 "zPol2In.=" << &vecz2resInner <<
2076 "yPol1Out.=" << &vecy1resOuter <<
2077 "zPol1Out.=" << &vecz1resOuter <<
2078 "yPol1InOut.="<< &vecy1resIO <<
2079 "zPol1InOut.="<< &vecz1resIO <<
2080 "yPol2Out.=" << &vecy2resOuter <<
2081 "zPol2Out.=" << &vecz2resOuter <<
2087 //==========================//
2088 // Fill Residual Histograms //
2089 //==========================//
2090 if (!fHisNclIn) MakeFitHistos();
2092 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2093 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2094 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2095 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2096 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2097 TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2098 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2099 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2101 for (Int_t irow=158;irow>-1;--irow) {
2102 if (vecSec[irow]==-1)continue; // no cluster info
2103 if (isReject[irow]>0.5) continue; //
2104 Double_t ycl = vecClY[irow];
2105 Double_t yfit = vecY1[irow];
2106 Double_t yfit2 = vecY2[irow];
2107 Double_t x = vecX[irow];
2108 Double_t yabsbeam = -1000;
2109 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2110 yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2111 else if(innerSector==fBeamSectorInner[id])
2112 yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2114 // Double_t yfit3 = vecY2[irow];
2115 Double_t zcl = vecClZ[irow];
2116 Double_t zfit = vecZ1[irow];
2117 Double_t zfit2 = vecZ2[irow];
2118 //Double_t zfit3 = vecZ2[irow];
2121 // The expressions for zcorrected has been obtained by
2122 // inverting the fits in the FitDriftV() method (ignoring the
2123 // global y dependence for now):
2125 // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2127 // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2130 // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
2132 // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2134 Double_t dzabs = -1000;
2135 Double_t zcorrected = -1000;
2136 if (ltrp->GetSide()==0){
2137 if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2138 // ignore global y dependence for now
2140 if(!fUseFixedDriftV)
2141 zcorrected = (zcl + (*fFitAside)[0] -
2142 (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2144 zcorrected = (zcl + fFixedFitAside0 -
2145 (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2146 // zcorrected = zcl;
2147 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2148 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2149 else if(innerSector==fBeamSectorInner[id])
2150 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2153 if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2155 if(!fUseFixedDriftV)
2156 zcorrected = (zcl - (*fFitCside)[0] +
2157 (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2159 zcorrected = (zcl - fFixedFitCside0 +
2160 (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2162 // zcorrected = zcl;
2163 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2164 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2165 else if(innerSector==fBeamSectorInner[id])
2166 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2170 if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2172 profy->Fill(irow,ycl-yfit);
2173 profy2->Fill(irow,ycl-yfit2);
2175 fHisYAbsErrors->Fill(id);
2176 // profyabs->Fill(irow,-0.99);
2178 profyabs->Fill(irow,ycl-yabsbeam);
2180 // profy3->Fill(irow,ycl-yfit3);
2183 profz->Fill(irow,zcl-zfit);
2184 profz2->Fill(irow,zcl-zfit2);
2185 //profz3->Fill(irow,zcl-zfit3);
2188 fHisZAbsErrors->Fill(id);
2190 profzabs->Fill(irow,dzabs);
2196 // Fill laser fit histograms
2198 Float_t dedx = track->GetdEdx();
2199 if (nclI>20&&nclO>20){
2200 fHisNclIn->Fill(id,nclI); //->Number of clusters inner
2201 fHisNclOut->Fill(id,nclO); //->Number of clusters outer
2202 fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
2204 fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
2205 fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
2206 fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
2208 fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2209 fHisdZfit->Fill(id,fFitZ[id]);
2212 fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
2213 fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
2214 fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
2215 fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
2216 fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
2219 fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
2220 fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
2221 fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
2222 fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
2223 fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
2226 fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
2227 fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
2228 fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
2230 fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
2231 fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
2232 fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
2234 fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
2235 fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
2236 fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
2237 fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
2239 fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
2240 fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
2241 fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
2242 fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
2245 fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
2246 fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
2250 fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
2251 fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
2259 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2261 // Dump information about laser beams
2262 // isOK variable indicates usability of the beam
2263 // Beam is not usable if:
2264 // a. No entries in range (krmsCut0)
2265 // b. Big sperad (krmscut1)
2266 // c. RMSto fit sigma bigger then (kmultiCut)
2267 // d. Too big angular spread
2270 const Float_t krmsCut0=0.001;
2271 const Float_t krmsCut1=0.16;
2272 const Float_t kmultiCut=2;
2273 const Float_t kcutP0=0.002;
2274 AliMagF* magF= dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
2275 Double_t xyz[3]={90,0,10}; // tmp. global position
2276 Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical
2277 Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local
2279 AliTPCcalibLaser *laser = this;
2280 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2281 TF1 fg("fg","gaus");
2282 AliTPCParam * tpcparam = 0;
2283 // start set up for absolute residuals analysis
2285 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2286 tpcparam = calib->GetParameters();
2287 if (!tpcparam) tpcparam = new AliTPCParamSR;
2289 AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2291 Float_t bfield = 0, bz=0;
2294 Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
2295 current = grp->GetL3Current((AliGRPObject::Stats)0);
2296 bfield = polarity*5*current/30000.;
2297 bz = polarity*5*current/30000.;
2298 printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2301 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2302 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2303 TLinearFitter lfabsyInner(2);
2304 lfabsyInner.SetFormula("1 ++ x");
2305 TLinearFitter lfabszInner(2);
2306 lfabszInner.SetFormula("1 ++ x");
2308 TLinearFitter lfabsyOuter(2);
2309 lfabsyOuter.SetFormula("1 ++ x");
2310 TLinearFitter lfabszOuter(2);
2311 lfabszOuter.SetFormula("1 ++ x");
2312 // end set up for absolute residuals analysis
2316 for (Int_t id=0; id<336; id++){
2318 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
2319 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2320 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
2321 TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
2322 TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
2323 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2324 //if (!hisphi) continue;
2325 Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2326 //if (entries<minEntries) continue;
2328 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2330 AliTPCLaserTrack::LoadTracks();
2331 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2333 ltrp->UpdatePoints();
2334 pcstream->GetFile()->cd();
2335 if (hisphi) hisphi->Write();
2336 if (hisphiP) hisphiP->Write();
2337 if (hisZ) hisZ->Write();
2338 if (hisP3) hisP3->Write();
2339 if (hisP4) hisP4->Write();
2341 Float_t meanphi = hisphi->GetMean();
2342 Float_t rmsphi = hisphi->GetRMS();
2344 Float_t meanphiP = hisphiP->GetMean();
2345 Float_t rmsphiP = hisphiP->GetRMS();
2346 Float_t meanZ = hisZ->GetMean();
2347 Float_t rmsZ = hisZ->GetRMS();
2348 if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2349 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2350 Double_t gphi1 = fg.GetParameter(1);
2351 Double_t gphi2 = fg.GetParameter(2);
2352 if (hisphiP->GetRMS()>0)
2353 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2354 Double_t gphiP1 = fg.GetParameter(1);
2355 Double_t gphiP2 = fg.GetParameter(2);
2357 if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2358 hisZ->Fit(&fg,"","");
2359 Double_t gz1 = fg.GetParameter(1);
2360 Double_t gz2 = fg.GetParameter(2);
2362 if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2363 hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2364 Double_t gp31 = fg.GetParameter(1);
2365 Double_t gp32 = fg.GetParameter(2);
2366 Double_t meanp3 = hisP3->GetMean();
2367 Double_t rmsp3 = hisP3->GetRMS();
2369 if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2370 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2371 Double_t gp41 = fg.GetParameter(1);
2372 Double_t gp42 = fg.GetParameter(2);
2373 Double_t meanp4 = hisP4->GetMean();
2374 Double_t rmsp4 = hisP4->GetRMS();
2376 Float_t meanS=hisS->GetMean();
2381 ltrp->GetPxPyPz(lpxyz);
2383 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2384 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2385 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2386 if (gphiP2>kcutP0) isOK=kFALSE;
2392 his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2393 Float_t mnclIn = his->GetMean();
2395 his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2396 Float_t mnclOut = his->GetMean();
2398 his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2399 Float_t mnclIO = his->GetMean();
2401 his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2402 Float_t mLclIn = his->GetMean();
2404 his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2405 Float_t mLclOut = his->GetMean();
2407 his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2408 Float_t mLclIO = his->GetMean();
2411 his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2412 Float_t mdEdx = his->GetMean();
2418 his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2419 Float_t mChi2YIn1= his->GetMean();
2421 his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2422 Float_t mChi2YOut1 = his->GetMean();
2424 his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2425 Float_t mChi2YIn2 = his->GetMean();
2427 his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2428 Float_t mChi2YOut2 = his->GetMean();
2430 his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2431 Float_t mChi2YIO1 = his->GetMean();
2433 his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2434 Float_t mChi2ZIn1 = his->GetMean();
2436 his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2437 Float_t mChi2ZOut1 = his->GetMean();
2439 his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2440 Float_t mChi2ZIn2 = his->GetMean();
2442 his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2443 Float_t mChi2ZOut2 = his->GetMean();
2445 his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2446 Float_t mChi2ZIO1 = his->GetMean();
2451 his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2452 Float_t edZfit = his->GetEntries();
2453 Float_t mdZfit = his->GetMean();
2454 Float_t rdZfit = his->GetRMS();
2457 his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2458 Float_t ePy1vP0 = his->GetEntries();
2459 Float_t mPy1vP0 = his->GetMean();
2460 Float_t rPy1vP0 = his->GetRMS();
2463 his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2464 Float_t ePy2vP0 = his->GetEntries();
2465 Float_t mPy2vP0 = his->GetMean();
2466 Float_t rPy2vP0 = his->GetRMS();
2469 his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2470 Float_t ePy3vP0 = his->GetEntries();
2471 Float_t mPy3vP0 = his->GetMean();
2472 Float_t rPy3vP0 = his->GetRMS();
2475 his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2476 Float_t ePy1vP1 = his->GetEntries();
2477 Float_t mPy1vP1 = his->GetMean();
2478 Float_t rPy1vP1 = his->GetRMS();
2481 his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2482 Float_t ePy2vP1 = his->GetEntries();
2483 Float_t mPy2vP1 = his->GetMean();
2484 Float_t rPy2vP1 = his->GetRMS();
2487 his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2488 Float_t ePy3vP1 = his->GetEntries();
2489 Float_t mPy3vP1 = his->GetMean();
2490 Float_t rPy3vP1 = his->GetRMS();
2493 his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2494 Float_t ePy2vP2In = his->GetEntries();
2495 Float_t mPy2vP2In = his->GetMean();
2496 Float_t rPy2vP2In = his->GetRMS();
2499 his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2500 Float_t ePy2vP2Out = his->GetEntries();
2501 Float_t mPy2vP2Out = his->GetMean();
2502 Float_t rPy2vP2Out = his->GetRMS();
2505 his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2506 Float_t ePy3vP2IO = his->GetEntries();
2507 Float_t mPy3vP2IO = his->GetMean();
2508 Float_t rPy3vP2IO = his->GetRMS();
2513 his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2514 Float_t ePz1vP0 = his->GetEntries();
2515 Float_t mPz1vP0 = his->GetMean();
2516 Float_t rPz1vP0 = his->GetRMS();
2519 his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2520 Float_t ePz2vP0 = his->GetEntries();
2521 Float_t mPz2vP0 = his->GetMean();
2522 Float_t rPz2vP0 = his->GetRMS();
2525 his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2526 Float_t ePz3vP0 = his->GetEntries();
2527 Float_t mPz3vP0 = his->GetMean();
2528 Float_t rPz3vP0 = his->GetRMS();
2531 his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2532 Float_t ePz1vP1 = his->GetEntries();
2533 Float_t mPz1vP1 = his->GetMean();
2534 Float_t rPz1vP1 = his->GetRMS();
2537 his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2538 Float_t ePz2vP1 = his->GetEntries();
2539 Float_t mPz2vP1 = his->GetMean();
2540 Float_t rPz2vP1 = his->GetRMS();
2543 his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2544 Float_t ePz3vP1 = his->GetEntries();
2545 Float_t mPz3vP1 = his->GetMean();
2546 Float_t rPz3vP1 = his->GetRMS();
2549 his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2550 Float_t ePz2vP2In = his->GetEntries();
2551 Float_t mPz2vP2In = his->GetMean();
2552 Float_t rPz2vP2In = his->GetRMS();
2555 his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2556 Float_t ePz2vP2Out = his->GetEntries();
2557 Float_t mPz2vP2Out = his->GetMean();
2558 Float_t rPz2vP2Out = his->GetRMS();
2561 his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2562 Float_t ePz3vP2IO = his->GetEntries();
2563 Float_t mPz3vP2IO = his->GetMean();
2564 Float_t rPz3vP2IO = his->GetRMS();
2567 // Fit absolute laser residuals
2568 TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2569 TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2571 Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2572 Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2574 TVectorD vecX(159); // X
2575 TVectorD vecY(159); // Y
2576 TVectorD vecR(159); // R
2577 TVectorD vecDY(159); // absolute residuals in Y
2578 TVectorD vecDZ(159); // absolute residuals in Z
2579 TVectorD vecN(159); // number of clusters
2580 TVectorD vecEy(159); //error y
2581 TVectorD vecEz(159); //error z
2582 TVectorD vecPhi(159); // local tangent
2583 TVectorD vecPhiR(159); // local tangent
2584 // magnetic field integrals
2585 TVectorD vecIBR(159); // radial
2586 TVectorD vecIBRPhi(159); // r-phi
2587 TVectorD vecIBLX(159); // local x
2588 TVectorD vecIBLY(159); // local y
2589 TVectorD vecIBGX(159); // local x
2590 TVectorD vecIBGY(159); // local y
2591 TVectorD vecIBZ(159); // z
2593 for (Int_t irow=0;irow<159;irow++){
2601 Double_t gx =(*(ltrp->fVecGX))[irow];
2602 Double_t gy =(*(ltrp->fVecGY))[irow];
2603 Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
2604 Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2605 Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2606 xyz[2]=(*(ltrp->fVecGZ))[irow];
2607 xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2608 xyz[1]=TMath::ATan2(gy,gx);
2609 Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2611 magF->GetTPCIntCyl(xyz,bxyz);
2612 magF->GetTPCInt(gxyz,bgxyz);
2613 vecIBR[irow]=bxyz[0];
2614 vecIBRPhi[irow]=bxyz[1];
2616 vecIBGX[irow]=bgxyz[0];
2617 vecIBGY[irow]=bgxyz[1];
2619 vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
2620 vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2623 vecIBZ[irow]=bxyz[2];
2628 lfabsyInner.ClearPoints();
2629 lfabszInner.ClearPoints();
2630 lfabsyOuter.ClearPoints();
2631 lfabszOuter.ClearPoints();
2634 Float_t yAbsInOffset = -100;
2635 Float_t yAbsInSlope = -100;
2636 Float_t yAbsInDeltay = -100;
2638 Float_t zAbsInOffset = -100;
2639 Float_t zAbsInSlope = -100;
2640 Float_t zAbsInDeltaz = -100;
2642 Float_t yAbsOutOffset = -100;
2643 Float_t yAbsOutSlope = -100;
2644 Float_t yAbsOutDeltay = -100;
2646 Float_t zAbsOutOffset = -100;
2647 Float_t zAbsOutSlope = -100;
2648 Float_t zAbsOutDeltaz = -100;
2650 Float_t lasTanPhiLocIn = -100;
2651 Float_t lasTanPhiLocOut = -100;
2653 if(histAbsY && histAbsY->GetEntries()>0) {
2655 Double_t rotAngOut = 10;
2656 Double_t rotAngIn = 10;
2657 if((secInner%36)!=(secOuter%36))
2658 rotAngIn += 20; // 30 degrees
2660 // Calculate laser mirror X position in local frame
2661 Double_t laserposOut =
2662 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2663 Double_t laserposIn =
2664 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2666 // Calculate laser tan phi in local frame
2667 lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2668 if(lasTanPhiLocOut<0) {
2669 lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2670 lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2673 lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2674 lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2677 lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2678 lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2680 TProfile* yprof = histAbsY->ProfileX("yprof");
2681 TProfile* zprof = histAbsZ->ProfileX("zprof");
2683 for(Int_t bin = 1; bin<=159; bin++) {
2685 if(yprof->GetBinEntries(bin)<10&&
2686 zprof->GetBinEntries(bin)<10) {
2690 // There is a problem in defining inner and outer sectors for
2691 // the outer beams (0 and 6) where both sectors are OROCs. To
2692 // make sure there is no overlap row 94 to 99 are cutted.
2693 if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2696 Int_t row = (bin-1);
2700 Bool_t isOuter = kTRUE;
2701 Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2704 (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2707 sector = TMath::Nint(fBeamSectorInner[id]);
2711 Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2712 vecN[bin-1] =yprof->GetBinEntries(bin);
2713 vecEy[bin-1]=yprof->GetBinError(bin);
2714 vecEz[bin-1]=zprof->GetBinError(bin);
2716 vecDY[bin-1] = yprof->GetBinContent(bin);
2717 vecDZ[bin-1] = zprof->GetBinContent(bin);
2718 if(!isOuter) { // inner
2719 vecPhi[bin-1]=lasTanPhiLocIn;
2720 // Calculate local y from residual and database
2721 Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2724 Double_t r = TMath::Sqrt(x*x + y*y);
2726 // Find angle between laser vector and R vector
2727 // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2728 Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2729 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2730 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2731 if(lasTanPhiLocIn<0)
2732 vecPhiR[bin-1]*=-1; // to have the same sign
2734 if(yprof->GetBinEntries(bin)>=10) {
2735 lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2736 TMath::Max(yprof->GetBinError(bin), 0.001));
2738 if(zprof->GetBinEntries(bin)>=10) {
2739 lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2740 TMath::Max(zprof->GetBinError(bin), 0.001));
2743 vecPhi[bin-1]=lasTanPhiLocOut;
2744 // Calculate local y from residual and database
2745 Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2748 Double_t r = TMath::Sqrt(x*x + y*y);
2751 Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2752 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2753 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2754 if(lasTanPhiLocOut<0)
2755 vecPhiR[bin-1]*=-1; // to have the same sign
2757 if(yprof->GetBinEntries(bin)>=10) {
2758 lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2759 TMath::Max(yprof->GetBinError(bin), 0.001));
2761 if(zprof->GetBinEntries(bin)>=10) {
2762 lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2763 TMath::Max(zprof->GetBinError(bin), 0.001));
2770 delete yprof; delete zprof;
2773 // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2774 nClInY = lfabsyInner.GetNpoints();
2775 if(lfabsyInner.GetNpoints()>10) {
2776 lfabsyInner.EvalRobust(0.95);
2779 lfabsyInner.GetParameters(result);
2780 yAbsInOffset = result[0];
2781 yAbsInSlope = result[1];
2782 yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2784 nClInZ = lfabszInner.GetNpoints();
2785 if(lfabszInner.GetNpoints()>10) {
2786 lfabszInner.EvalRobust(0.95);
2789 lfabszInner.GetParameters(result);
2790 zAbsInOffset = result[0];
2791 zAbsInSlope = result[1];
2792 zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2794 nClOutY = lfabsyOuter.GetNpoints();
2795 if(lfabsyOuter.GetNpoints()>10) {
2796 lfabsyOuter.EvalRobust(0.95);
2799 lfabsyOuter.GetParameters(result);
2800 yAbsOutOffset = result[0];
2801 yAbsOutSlope = result[1];
2802 yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2804 nClOutZ = lfabszOuter.GetNpoints();
2805 if(lfabszOuter.GetNpoints()>10) {
2806 lfabszOuter.EvalRobust(0.95);
2809 lfabszOuter.GetParameters(result);
2810 zAbsOutOffset = result[0];
2811 zAbsOutSlope = result[1];
2812 zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2818 Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2819 Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2820 Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2821 Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2823 Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2824 Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2826 Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2827 Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2830 (*pcstream)<<"Mean"<<
2843 "id="<<id<< // track id
2844 "entries="<<entries<< // number of entries
2845 "bz="<<bfield<< // bfield
2846 "LTr.="<<ltrp<< // refernece track
2848 "mnclIn="<<mnclIn<< // mean number of clusters in inner
2849 "mnclOut="<<mnclOut<< // mean number of clusters in outer
2850 "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2851 "mLclIn="<<mLclIn<< // mean number of clusters in inner
2852 "mLclOut="<<mLclOut<< // mean number of clusters in outer
2853 "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2854 "mdEdx="<<mdEdx<< // mean dEdx
2855 "edZfit="<<edZfit<< // entries z fit
2856 "mdZfit="<<mdZfit<< // mean z fit
2857 "rdZfit="<<rdZfit<< // RMS z fit
2860 "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2861 "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2862 "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2863 "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2864 "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2865 "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2866 "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2867 "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2868 "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2869 "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2872 "ePy1vP0="<<ePy1vP0<<
2873 "mPy1vP0="<<mPy1vP0<<
2874 "rPy1vP0="<<rPy1vP0<<
2875 "ePy2vP0="<<ePy2vP0<<
2876 "mPy2vP0="<<mPy2vP0<<
2877 "rPy2vP0="<<rPy2vP0<<
2878 "ePy3vP0="<<ePy3vP0<<
2879 "mPy3vP0="<<mPy3vP0<<
2880 "rPy3vP0="<<rPy3vP0<<
2881 "ePy1vP1="<<ePy1vP1<<
2882 "mPy1vP1="<<mPy1vP1<<
2883 "rPy1vP1="<<rPy1vP1<<
2884 "ePy2vP1="<<ePy2vP1<<
2885 "mPy2vP1="<<mPy2vP1<<
2886 "rPy2vP1="<<rPy2vP1<<
2887 "ePy3vP1="<<ePy3vP1<<
2888 "mPy3vP1="<<mPy3vP1<<
2889 "rPy3vP1="<<rPy3vP1<<
2890 "ePy2vP2In="<<ePy2vP2In<<
2891 "mPy2vP2In="<<mPy2vP2In<<
2892 "rPy2vP2In="<<rPy2vP2In<<
2893 "ePy2vP2Out="<<ePy2vP2Out<<
2894 "mPy2vP2Out="<<mPy2vP2Out<<
2895 "rPy2vP2Out="<<rPy2vP2Out<<
2896 "ePy3vP2IO="<<ePy3vP2IO<<
2897 "mPy3vP2IO="<<mPy3vP2IO<<
2898 "rPy3vP2IO="<<rPy3vP2IO<<
2901 "ePz1vP0="<<ePz1vP0<<
2902 "mPz1vP0="<<mPz1vP0<<
2903 "rPz1vP0="<<rPz1vP0<<
2904 "ePz2vP0="<<ePz2vP0<<
2905 "mPz2vP0="<<mPz2vP0<<
2906 "rPz2vP0="<<rPz2vP0<<
2907 "ePz3vP0="<<ePz3vP0<<
2908 "mPz3vP0="<<mPz3vP0<<
2909 "rPz3vP0="<<rPz3vP0<<
2910 "ePz1vP1="<<ePz1vP1<<
2911 "mPz1vP1="<<mPz1vP1<<
2912 "rPz1vP1="<<rPz1vP1<<
2913 "ePz2vP1="<<ePz2vP1<<
2914 "mPz2vP1="<<mPz2vP1<<
2915 "rPz2vP1="<<rPz2vP1<<
2916 "ePz3vP1="<<ePz3vP1<<
2917 "mPz3vP1="<<mPz3vP1<<
2918 "rPz3vP1="<<rPz3vP1<<
2919 "ePz2vP2In="<<ePz2vP2In<<
2920 "mPz2vP2In="<<mPz2vP2In<<
2921 "rPz2vP2In="<<rPz2vP2In<<
2922 "ePz2vP2Out="<<ePz2vP2Out<<
2923 "mPz2vP2Out="<<mPz2vP2Out<<
2924 "rPz2vP2Out="<<rPz2vP2Out<<
2925 "ePz3vP2IO="<<ePz3vP2IO<<
2926 "mPz3vP2IO="<<mPz3vP2IO<<
2927 "rPz3vP2IO="<<rPz3vP2IO<<
2931 "lx0="<<lxyz[0]<< // reference x
2932 "lx1="<<lxyz[1]<< // reference y
2933 "lx2="<<lxyz[2]<< // refernece z
2934 "lpx0="<<lpxyz[0]<< // reference x
2935 "lpx1="<<lpxyz[1]<< // reference y
2936 "lpx2="<<lpxyz[2]<< // refernece z
2940 "mphi="<<meanphi<< //
2941 "rmsphi="<<rmsphi<< //
2945 "mphiP="<<meanphiP<< //
2946 "rmsphiP="<<rmsphiP<< //
2955 "gp31="<<gp31<< //gaus mean - tgl
2956 "gp32="<<gp32<< //gaus rms -tgl
2957 "meanp3="<<meanp3<< //mean - tgl
2958 "rmsp3="<<rmsp3<< //rms -tgl
2959 "gp41="<<gp41<< //gaus mean - P4
2960 "gp42="<<gp42<< //gaus rms - P4
2961 "meanp4="<<meanp4<< //mean - P4
2962 "rmsp4="<<rmsp4<< //rms - P4
2963 // Parameters from abs res analysis
2964 "SecIn="<<secInner<< // inner sector
2965 "SecOut="<<secOuter<< // outer sector
2966 "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2967 "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
2968 "ibr.="<<&vecIBR<< // radial filed integral
2969 "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
2970 "ibr.="<<&vecIBR<< // radial filed integral
2971 "ibz.="<<&vecIBZ<< // z filed integral
2973 "iblx.="<<&vecIBLX<< // local bx integral
2974 "ibly.="<<&vecIBLY<< // local by integral
2975 "ibgx.="<<&vecIBGX<< // global bx integral
2976 "ibgy.="<<&vecIBGY<< // global by integral
2978 "X.="<<&vecX<< // local x
2979 "Y.="<<&vecY<< // local y
2980 "R.="<<&vecR<< // radius
2981 "dY.="<<&vecDY<< // abs y residuals
2982 "dZ.="<<&vecDZ<< // abs z residuals
2983 "eY.="<<&vecEy<< // error of y residuals
2984 "eZ.="<<&vecEz<< // error of z residuals
2985 "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
2986 "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
2987 "nCl.="<<&vecN<< // number of clusters
2989 "nClInY="<<nClInY<< // Number of clusters for inner y
2990 "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2991 "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
2992 "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2993 "nClInZ="<<nClInZ<< // Number of clusters for inner z
2994 "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2995 "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
2996 "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2998 "nClOutY="<<nClOutY<< // Number of clusters for outer y
2999 "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
3000 "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
3001 "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
3002 "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
3003 "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
3004 "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
3005 "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
3012 TFile fmean("laserMean.root");
3022 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3026 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3027 TFile * f = pcstream->GetFile();
3029 f->mkdir("dirphiP");
3031 TF1 fp("p1","pol1");
3036 char grnamefull[1000];
3039 Double_t mphiP[100];
3040 Double_t smphi[100];
3041 Double_t smphiP[100];
3052 for (Int_t id=0; id<336; id++){
3054 snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
3055 Int_t entries = chain->Draw("bz",cut,"goff");
3056 if (entries<3) continue;
3057 AliTPCLaserTrack *ltrp = 0;
3058 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3059 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3063 ltrp->GetPxPyPz(lpxyz);
3065 chain->Draw("bz",cut,"goff");
3066 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3067 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3068 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3070 chain->Draw("gphi1",cut,"goff");
3071 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3072 chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3073 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3075 chain->Draw("gphiP1",cut,"goff");
3076 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3077 chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3078 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3080 chain->Draw("gz1",cut,"goff");
3081 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3082 chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3083 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3086 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3087 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3094 for (Int_t ientry=0; ientry<entries; ientry++){
3095 if (TMath::Abs(bz[ientry])<0.05){
3096 phiB0 = mphi[ientry];
3097 phiPB0 = mphiP[ientry];
3101 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3104 pphi[0] = fp.GetParameter(0); // offset
3105 pphi[1] = fp.GetParameter(1); // slope
3106 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3107 sprintf(grname,"phi_id%d",id);
3108 grphi->SetName(grname); grphi->SetTitle(grnamefull);
3109 grphi->GetXaxis()->SetTitle("b_{z} (T)");
3110 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3111 grphi->SetMaximum(1.2);
3112 grphi->SetMinimum(-1.2);
3116 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3119 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3122 pphiP[0] = fp.GetParameter(0); // offset
3123 pphiP[1] = fp.GetParameter(1); // slope
3124 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3125 sprintf(grname,"phiP_id%d",id);
3126 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3127 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3128 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3129 grphiP->SetMaximum(pphiP[0]+0.005);
3130 grphiP->SetMinimum(pphiP[0]-0.005);
3132 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3137 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3140 pmZ[0] = fp.GetParameter(0); // offset
3141 pmZ[1] = fp.GetParameter(1); // slope
3142 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3143 sprintf(grname,"mZ_id%d",id);
3144 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3145 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3146 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3148 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3154 for (Int_t ientry=0; ientry<entries; ientry++){
3155 (*pcstream)<<"Mean"<<
3158 "entries="<<entries<<
3160 "lx0="<<lxyz[0]<< // reference x
3161 "lx1="<<lxyz[1]<< // reference y
3162 "lx2="<<lxyz[2]<< // refernece z
3163 "lpx0="<<lpxyz[0]<< // reference x
3164 "lpx1="<<lpxyz[1]<< // reference y
3165 "lpx2="<<lpxyz[2]<< // refernece z
3167 "phiB0="<<phiB0<< // position shift at 0 field
3168 "phiPB0="<<phiPB0<< // angular shift at 0 field
3169 "zB0="<<zB0<< // z shift for 0 field
3171 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3172 "pphi0="<<pphi[0]<< // offset
3173 "pphi1="<<pphi[1]<< // slope
3174 "pphi2="<<pphi[2]<< // norm chi2
3176 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3177 "pphiP0="<<pphiP[0]<< // offset
3178 "pphiP1="<<pphiP[1]<< // slope
3179 "pphiP2="<<pphiP[2]<< // norm chi2
3181 "gz1="<<mZ[ientry]<<
3182 "pmZ0="<<pmZ[0]<< // offset
3183 "pmZ1="<<pmZ[1]<< // slope
3184 "pmZ2="<<pmZ[2]<< // norm chi2
3194 void AliTPCcalibLaser::Analyze(){
3201 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3203 TIterator* iter = li->MakeIterator();
3204 AliTPCcalibLaser* cal = 0;
3205 static Int_t counter0=0;
3206 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3207 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3208 // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3211 printf("Marging number %d\n", counter0);
3214 MergeFitHistos(cal);
3219 // TProfile *hp=0x0;
3220 //TProfile *hpm=0x0;
3222 for (Int_t id=0; id<336; id++){
3223 // merge fDeltaZ histograms
3224 hm = (TH1F*)cal->fDeltaZ.At(id);
3225 h = (TH1F*)fDeltaZ.At(id);
3227 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3229 fDeltaZ.AddAt(h,id);
3232 // merge fP3 histograms
3233 hm = (TH1F*)cal->fDeltaP3.At(id);
3234 h = (TH1F*)fDeltaP3.At(id);
3236 h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3238 fDeltaP3.AddAt(h,id);
3241 // merge fP4 histograms
3242 hm = (TH1F*)cal->fDeltaP4.At(id);
3243 h = (TH1F*)fDeltaP4.At(id);
3245 h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3247 fDeltaP4.AddAt(h,id);
3252 // merge fDeltaPhi histograms
3253 hm = (TH1F*)cal->fDeltaPhi.At(id);
3254 h = (TH1F*)fDeltaPhi.At(id);
3256 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3258 fDeltaPhi.AddAt(h,id);
3261 // merge fDeltaPhiP histograms
3262 hm = (TH1F*)cal->fDeltaPhiP.At(id);
3263 h = (TH1F*)fDeltaPhiP.At(id);
3265 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3267 fDeltaPhiP.AddAt(h,id);
3270 // merge fSignals histograms
3271 hm = (TH1F*)cal->fSignals.At(id);
3272 h = (TH1F*)fSignals.At(id);
3274 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3276 fSignals.AddAt(h,id);
3281 // merge ProfileY histograms -0
3282 h2m = (TH2F*)cal->fDeltaYres.At(id);
3283 h2 = (TH2F*)fDeltaYres.At(id);
3284 if (h2m&&h2) h2->Add(h2m);
3286 h2m = (TH2F*)cal->fDeltaZres.At(id);
3287 h2 = (TH2F*)fDeltaZres.At(id);
3288 if (h2m&&h2) h2->Add(h2m);
3289 // merge ProfileY histograms - 2
3290 h2m = (TH2F*)cal->fDeltaYres2.At(id);
3291 h2 = (TH2F*)fDeltaYres2.At(id);
3292 if (h2m&&h2) h2->Add(h2m);
3294 h2m = (TH2F*)cal->fDeltaZres2.At(id);
3295 h2 = (TH2F*)fDeltaZres2.At(id);
3296 if (h2m&&h2) h2->Add(h2m);
3298 // merge ProfileY histograms - abs
3299 h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3300 h2 = (TH2F*)fDeltaYresAbs.At(id);
3301 if (h2m&&h2) h2->Add(h2m);
3303 h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3304 h2 = (TH2F*)fDeltaZresAbs.At(id);
3305 if (h2m&&h2) h2->Add(h2m);
3306 // merge ProfileY histograms - 3
3307 //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3308 //h2 = (TH2F*)fDeltaYres3.At(id);
3309 //if (h2m) h2->Add(h2m);
3311 //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3312 //h2 = (TH2F*)fDeltaZres3.At(id);
3313 //if (h2m) h->Add(h2m);
3321 void AliTPCcalibLaser::MakeFitHistos(){
3323 // Make a fit histograms
3325 // Number of clusters
3327 //TH2F *fHisNclIn; //->Number of clusters inner
3328 //TH2F *fHisNclOut; //->Number of clusters outer
3329 //TH2F *fHisNclIO; //->Number of cluster inner outer
3330 // TH2F *fHisdEdx; //->dEdx histo
3331 fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3332 fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3333 fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3335 fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3336 fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3337 fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3339 fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3340 fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3345 // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3346 // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3347 // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3348 // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3349 // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3350 fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3351 fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3352 fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3353 fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3354 fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3355 // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3356 // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3357 // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3358 // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3359 // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3360 fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3361 fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3362 fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3363 fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3364 fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3369 // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3370 // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3371 // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3372 // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3373 // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3374 // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3375 // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3376 // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3377 // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3378 fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3379 fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3380 fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3381 fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3382 fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3383 fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3384 fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3385 fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3386 fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3389 // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3390 // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3391 // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3392 // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3393 // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3394 // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3395 // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3396 // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3397 // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3398 fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3399 fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3400 fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3401 fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3402 fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3403 fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3404 fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3405 fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3406 fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3408 fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3409 fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3411 fHisNclIn->SetDirectory(0); //->Number of clusters inner
3412 fHisNclOut->SetDirectory(0); //->Number of clusters outer
3413 fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3414 fHisLclIn->SetDirectory(0); //->Level arm inner
3415 fHisLclOut->SetDirectory(0); //->Level arm outer
3416 fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3417 fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3418 fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3421 fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3422 fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3423 fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3424 fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3425 fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3426 fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3427 fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3428 fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3429 fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3430 fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3433 fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3434 fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3435 fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3436 fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3437 fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3438 fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3439 fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3440 fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3441 fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3444 fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3445 fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3446 fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3447 fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3448 fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3449 fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3450 fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3451 fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3452 fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3454 fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3455 fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3462 for (Int_t id=0; id<336;id++){
3463 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3464 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3465 //TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3467 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
3468 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3469 TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3470 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3471 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3473 profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3474 profy->SetDirectory(0);
3475 fDeltaYres.AddAt(profy,id);
3476 profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3477 profy2->SetDirectory(0);
3478 fDeltaYres2.AddAt(profy2,id);
3479 if(!fUseFixedDriftV)
3480 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
3482 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
3483 profyabs->SetDirectory(0);
3484 fDeltaYresAbs.AddAt(profyabs,id);
3485 //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3486 //profy3->SetDirectory(0);
3487 //fDeltaYres3.AddAt(profy3,id);
3490 profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3491 profz->SetDirectory(0);
3492 fDeltaZres.AddAt(profz,id);
3493 profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3494 profz2->SetDirectory(0);
3495 fDeltaZres2.AddAt(profz2,id);
3496 if(!fUseFixedDriftV)
3497 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
3499 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
3500 profzabs->SetDirectory(0);
3501 fDeltaZresAbs.AddAt(profzabs,id);
3502 //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3503 //profz3->SetDirectory(0);
3504 //fDeltaZres3.AddAt(profz3,id);
3509 for (Int_t id=0; id<336;id++){
3510 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3511 //TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3513 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3515 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3516 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3517 TH1F * hisSignal = (TH1F*)fSignals.At(id);
3520 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3521 hisdz->SetDirectory(0);
3522 fDeltaZ.AddAt(hisdz,id);
3524 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3525 hisP3->SetDirectory(0);
3526 fDeltaP3.AddAt(hisP3,id);
3528 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3529 hisP4->SetDirectory(0);
3530 fDeltaP4.AddAt(hisP4,id);
3533 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3534 hisdphi->SetDirectory(0);
3535 fDeltaPhi.AddAt(hisdphi,id);
3537 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3538 hisdphiP->SetDirectory(0);
3539 fDeltaPhiP.AddAt(hisdphiP,id);
3540 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3541 hisSignal->SetDirectory(0);
3542 fSignals.AddAt(hisSignal,id);
3546 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3547 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3548 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3549 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3553 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3554 Int_t binsLaser[12]= {336, //id
3566 Double_t xminLaser[12]= {0, //id
3578 Double_t xmaxLaser[12]= {336, //id
3591 TString nameLaser[12]= {"id",
3603 TString titleLaser[12]= {"id",
3615 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3616 for (Int_t iaxis=1; iaxis<12; iaxis++){
3617 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3618 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3622 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3624 // Merge content of histograms
3626 // Only first histogram is checked - all other should be the same
3627 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3629 if (!laser->fHisNclIn) return; // empty histograms
3630 if (!fHisNclIn) MakeFitHistos();
3632 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3633 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3634 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3635 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3636 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3637 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3638 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3639 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3642 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3643 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3644 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3645 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3646 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3647 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3648 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3649 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3650 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3651 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3654 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3655 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3656 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3657 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3658 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3659 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3660 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3661 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3662 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3665 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3666 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3667 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3668 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3669 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3670 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3671 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3672 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3673 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3674 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3675 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3688 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3690 // Dump fit information - collect information from the streamers
3693 TChain * chainFit=0;
3694 TChain * chainTrack=0;
3697 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3698 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3699 AliXRDPROOFtoolkit tool;
3700 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3701 chainTrack->Lookup();
3702 chainTrack->SetProof(kTRUE);
3703 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3704 chainDrift->Lookup();
3705 chainDrift->SetProof(kTRUE);
3707 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3709 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3711 chainFit->SetProof(kTRUE);
3712 chain->SetProof(kTRUE);
3713 AliTPCLaserTrack::LoadTracks();
3714 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3720 TCut cutP3("abs(Tr.fP[3])<0.1");
3721 TCut cutP4("abs(Tr.fP[4])<0.5");
3722 TCut cutPx = cutP3+cutP4;
3723 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3724 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3725 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3726 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3728 TCut cutdEdx("sqrt(dEdx)>3");
3729 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3730 TCut cutN("nclO>20&&nclI>20");
3731 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3735 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3736 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3737 TCut cutClX("abs(Cl[].fX)>10");
3738 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3739 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3740 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3741 TCut cutQ("sqrt(Cl[].fMax)>4");
3742 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3746 TH1F * phisAccept = 0;
3748 TProfile * pdEdx = 0;
3756 TProfile * pNclI = 0;
3757 TProfile * pNclO = 0;
3759 TProfile * pchi2YIn =0;
3760 TProfile * pchi2ZIn =0;
3761 TProfile * pchi2YOut =0;
3762 TProfile * pchi2ZOut =0;
3763 TProfile * pchi2YInOut =0;
3764 TProfile * pchi2ZInOut =0;;
3766 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3767 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3768 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3769 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3770 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3771 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3773 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3774 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3777 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3778 pP0 = (TProfile*)gROOT->FindObject("hP0");
3779 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3780 pP1 = (TProfile*)gROOT->FindObject("hP1");
3781 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3782 pP2 = (TProfile*)gROOT->FindObject("hP2");
3783 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3784 pP3 = (TProfile*)gROOT->FindObject("hP3");
3785 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3786 pP4 = (TProfile*)gROOT->FindObject("hP4");
3788 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3789 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3790 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3791 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3794 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3795 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3796 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3797 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3798 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3799 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3800 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3801 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3802 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3803 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3804 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3805 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3807 // second derivatives
3809 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3810 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3811 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3812 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3813 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3814 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3816 phisPy2In->FitSlicesY();
3817 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3818 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3819 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3821 phisPy2Out->FitSlicesY();
3822 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3823 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3824 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3826 phisPy2InOut->FitSlicesY();
3827 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3828 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3829 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3831 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3832 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3833 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3834 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3835 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3836 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3838 phisPz2In->FitSlicesY();
3839 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3840 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3841 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3843 phisPz2Out->FitSlicesY();
3844 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3845 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3846 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3848 phisPz2InOut->FitSlicesY();
3849 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3850 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3851 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3858 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3859 for (Int_t ilaser=0; ilaser<336; ilaser++){
3860 Float_t all=phisAl->GetBinContent(ilaser+1);
3861 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3862 Float_t out=phisOut->GetBinContent(ilaser+1);
3863 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3864 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3865 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3866 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3867 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3868 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3871 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3872 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3874 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3875 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3876 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3877 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3878 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3879 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3881 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3882 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3883 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3885 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3886 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3887 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3889 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3890 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3891 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3893 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3894 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3895 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3897 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3898 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3899 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3901 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3902 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3903 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3905 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3906 (*pcstream)<<"Scan"<<
3920 "chi2YIn="<<chi2YIn<<
3921 "chi2YOut="<<chi2YOut<<
3922 "chi2YInOut="<<chi2YInOut<<
3923 "chi2ZIn="<<chi2ZIn<<
3924 "chi2ZOut="<<chi2ZOut<<
3925 "chi2ZInOut="<<chi2ZInOut<<
3927 "nPy2In="<<entriesPy2In<<
3928 "mPy2In="<<meanPy2In<<
3929 "sPy2In="<<sigmaPy2In<<
3931 "nPy2Out="<<entriesPy2Out<<
3932 "mPy2Out="<<meanPy2Out<<
3933 "sPy2Out="<<sigmaPy2Out<<
3935 "nPy2InOut="<<entriesPy2InOut<<
3936 "mPy2InOut="<<meanPy2InOut<<
3937 "sPy2InOut="<<sigmaPy2InOut<<
3939 "nPz2In="<<entriesPz2In<<
3940 "mPz2In="<<meanPz2In<<
3941 "sPz2In="<<sigmaPz2In<<
3943 "nPz2Out="<<entriesPz2Out<<
3944 "mPz2Out="<<meanPz2Out<<
3945 "sPz2Out="<<sigmaPz2Out<<
3947 "nPz2InOut="<<entriesPz2InOut<<
3948 "mPz2InOut="<<meanPz2InOut<<
3949 "sPz2InOut="<<sigmaPz2InOut<<
3957 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3958 TVectorD& meanSlope,
3959 TVectorD& sectorArray,
3962 // This method should ideally go in AliTPCLaser
3963 // option == 0 (pads outer - closest to beam)
3964 // option == 1 (pads inner)
3965 // option == 2 (time outer)
3966 // option == 3 (time inner)
3967 Int_t nFailures = 0;
3969 for(Int_t id = 0; id < 336; id++) {
3971 if (!AliTPCLaserTrack::GetTracks())
3972 AliTPCLaserTrack::LoadTracks();
3973 AliTPCLaserTrack *ltrp =
3974 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3976 AliExternalTrackParam trackParam(*ltrp);
3978 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3979 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3980 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3982 Double_t angle = trackParam.GetAlpha();
3984 angle += 2*TMath::Pi();
3985 if(trackParam.GetSnp()>0) // track points to sector "before"
3986 angle -= deltaangle*TMath::DegToRad();
3987 else // track points to sector "after"
3988 angle += deltaangle*TMath::DegToRad();
3990 Bool_t success = trackParam.Rotate(angle);
3993 // cout << "WARNING: Rotate failed for ID: " << id << endl;
3997 angle *= TMath::RadToDeg();
3999 Int_t sector = TMath::Nint((angle-10.0)/20.0);
4004 if(ltrp->GetSide()==1) // C side
4006 if(option==0 || option==2)
4008 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
4011 sectorArray[id] = sector;
4013 const Double_t x0 = 0;
4015 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
4016 Double_t slopez = trackParam.GetTgl();
4017 // One needs a factor sqrt(1+slopey**2) to take into account the
4018 // longer path length
4019 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
4020 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
4022 // Double_t offsetz = trackParam.GetZ();
4023 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4024 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4025 if(option==2 || option==3) {
4026 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4028 meanOffset[id] = offsety; meanSlope[id] = slopey;
4033 AliWarning(Form("Rotate method failed %d times", nFailures));
4038 TFile f("vscan.root");
4043 chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
4045 chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
4048 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);
4051 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)
4060 // check edge effects
4061 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4063 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4065 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)
4069 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4070 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")
4081 dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
4084 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)
4086 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)
4092 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)
4094 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)
4098 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)
4100 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)
4109 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutDY,"prof")
4111 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&<r.fP[1]<0"+cutA+cutDY,"prof")
4115 chainFit->Draw("LTr.fId","nclI>10",100000)
4117 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4119 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4121 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4131 gSystem->Load("libSTAT.so")
4132 TStatToolkit toolkit;
4138 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
4143 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
4144 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
4145 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
4146 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
4148 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
4149 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
4150 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
4151 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
4153 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
4154 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
4155 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
4156 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
4161 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4163 treeT->SetAlias("fit",strq0->Data());
4166 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4168 treeT->SetAlias("fitP",strqP->Data());
4171 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4172 treeT->SetAlias("fitD",strqDrift->Data());
4175 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4177 for (Int_t i=0; i<6;i++){
4178 treeT->SetLineColor(i+2);
4179 treeT->SetMarkerSize(1);
4180 treeT->SetMarkerStyle(22+i);
4181 treeT->SetMarkerColor(i+2);
4183 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4191 TTree * tree = (TTree*)f.Get("FitModels");
4193 TEventList listLFit0("listLFit0","listLFit0");
4194 TEventList listLFit1("listLFit1","listLFit1");
4195 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4196 tree->SetEventList(&listLFit0);
4201 gSystem->Load("libSTAT.so")
4202 TStatToolkit toolkit;
4208 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4209 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4213 fstring+="cos(dp)++";
4214 fstring+="sin(dp)++";
4215 fstring+="cos(dt)++";
4216 fstring+="sin(dt)++";
4218 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4231 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4232 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4233 AliXRDPROOFtoolkit tool;
4234 TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4235 chainTrack->Lookup();
4236 chainTrack->SetProof(kTRUE);
4238 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4240 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4242 chainFit->SetProof(kTRUE);
4243 chain->SetProof(kTRUE);
4247 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4248 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4249 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4250 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4252 TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4253 TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4254 TCut cutN("nclO>20&&nclI>20");
4255 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4259 TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4260 TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4261 TCut cutClX("abs(Cl.fX)>10");
4262 TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4263 TCut cutCl=cutClY+cutClZ+cutClX;
4266 // check edge effects
4267 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4269 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4271 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)
4275 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4276 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")