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(0)
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("laser.txt","* driftvN",1)
73 TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
75 TChain * chainDriftN = tool.MakeChainRandom("laser.txt.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"
119 #include "TTreeStream.h"
122 #include "AliTPCLaserTrack.h"
123 #include "AliTPCcalibDB.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCParamSR.h"
126 #include "TTimeStamp.h"
127 #include "AliDCSSensorArray.h"
128 #include "AliDCSSensor.h"
129 #include "AliGRPObject.h"
133 ClassImp(AliTPCcalibLaser)
135 AliTPCcalibLaser::AliTPCcalibLaser():
142 fTracksEsdParam(336),
152 fHisLaser(0), // N dim histogram of laser
153 fHisNclIn(0), //->Number of clusters inner
154 fHisNclOut(0), //->Number of clusters outer
155 fHisNclIO(0), //->Number of cluster inner outer
156 fHisLclIn(0), //->Level arm inner
157 fHisLclOut(0), //->Level arm outer
158 fHisLclIO(0), //->Number of cluster inner outer
159 fHisdEdx(0), //->dEdx histo
160 fHisdZfit(0), //->distance to the mirror after linear fit
163 fHisChi2YIn1(0), //->chi2 y inner - line
164 fHisChi2YOut1(0), //->chi2 y inner - line
165 fHisChi2YIn2(0), //->chi2 y inner - parabola
166 fHisChi2YOut2(0), //->chi2 y inner - parabola
167 fHisChi2YIO1(0), //->chi2 y IO - common
168 fHisChi2ZIn1(0), //->chi2 z inner - line
169 fHisChi2ZOut1(0), //->chi2 z inner - line
170 fHisChi2ZIn2(0), //->chi2 z inner - parabola
171 fHisChi2ZOut2(0), //->chi2 z inner - parabola
172 fHisChi2ZIO1(0), //->chi2 z IO - common
175 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
176 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
177 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
178 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
179 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
180 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
181 fHisPy2vP2In(0), //-> Curv P2inner - parabola
182 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
183 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
186 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
187 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
188 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
189 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
190 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
191 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
192 fHisPz2vP2In(0), //-> Curv P2inner - parabola
193 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
194 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
196 fDeltaYres(336), //->2D histo of residuals
197 fDeltaZres(336), //->2D histo fo residuals
198 fDeltaYres2(336), //->2D histo of residuals
199 fDeltaZres2(336), //->2D histo fo residuals
200 fDeltaYresAbs(336), //->2D histo of residuals
201 fHisYAbsErrors(0), //-> total errors per beam in the abs res y analysis
202 fDeltaZresAbs(336), //->2D histo of residuals
203 fHisZAbsErrors(0), //-> total errors per beam in the abs res z analysis
204 //fDeltaYres3(336), //->2D histo of residuals
205 //fDeltaZres3(336), //->2D histo fo residuals
206 fFitAside(new TVectorD(5)),
207 fFitCside(new TVectorD(5)),
208 fFitACside(new TVectorD(6)),
213 fBeamSectorOuter(336),
214 fBeamSectorInner(336),
215 fBeamOffsetYOuter(336),
216 fBeamSlopeYOuter(336),
217 fBeamOffsetYInner(336),
218 fBeamSlopeYInner(336),
219 fBeamOffsetZOuter(336),
220 fBeamSlopeZOuter(336),
221 fBeamOffsetZInner(336),
222 fBeamSlopeZInner(336),
223 fInverseSlopeZ(kTRUE),
225 fFixedFitAside0(0.0),
226 fFixedFitAside1(1.0),
227 fFixedFitCside0(0.0),
233 fTracksEsdParam.SetOwner(kTRUE);
236 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
243 fTracksEsdParam(336),
247 fDeltaZ(336), // array of histograms of delta z for each track
248 fDeltaP3(336), // array of histograms of delta z for each track
249 fDeltaP4(336), // array of histograms of P3 for each track
250 fDeltaPhi(336), // array of histograms of P4 for each track
251 fDeltaPhiP(336), // array of histograms of delta z for each track
252 fSignals(336), // array of dedx signals
255 fHisLaser(0), // N dim histogram of laser
256 fHisNclIn(0), //->Number of clusters inner
257 fHisNclOut(0), //->Number of clusters outer
258 fHisNclIO(0), //->Number of cluster inner outer
259 fHisLclIn(0), //->Level arm inner
260 fHisLclOut(0), //->Level arm outer
261 fHisLclIO(0), //->Number of cluster inner outer
262 fHisdEdx(0), //->dEdx histo
263 fHisdZfit(0), //->distance to the mirror after linear fit
266 fHisChi2YIn1(0), //->chi2 y inner - line
267 fHisChi2YOut1(0), //->chi2 y inner - line
268 fHisChi2YIn2(0), //->chi2 y inner - parabola
269 fHisChi2YOut2(0), //->chi2 y inner - parabola
270 fHisChi2YIO1(0), //->chi2 y IO - common
271 fHisChi2ZIn1(0), //->chi2 z inner - line
272 fHisChi2ZOut1(0), //->chi2 z inner - line
273 fHisChi2ZIn2(0), //->chi2 z inner - parabola
274 fHisChi2ZOut2(0), //->chi2 z inner - parabola
275 fHisChi2ZIO1(0), //->chi2 z IO - common
278 fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
279 fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
280 fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
281 fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
282 fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
283 fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
284 fHisPy2vP2In(0), //-> Curv P2inner - parabola
285 fHisPy2vP2Out(0), //-> Curv P2outer - parabola
286 fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
289 fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
290 fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
291 fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
292 fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
293 fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
294 fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
295 fHisPz2vP2In(0), //-> Curv P2inner - parabola
296 fHisPz2vP2Out(0), //-> Curv P2outer - parabola
297 fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
311 fFitAside(new TVectorD(5)), // drift fit - A side
312 fFitCside(new TVectorD(5)), // drift fit - C- side
313 fFitACside(new TVectorD(6)), // drift fit - AC- side
314 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
315 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
316 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
317 fNcuts(0), // number of cuts
318 fBeamSectorOuter(336),
319 fBeamSectorInner(336),
320 fBeamOffsetYOuter(336),
321 fBeamSlopeYOuter(336),
322 fBeamOffsetYInner(336),
323 fBeamSlopeYInner(336),
324 fBeamOffsetZOuter(336),
325 fBeamSlopeZOuter(336),
326 fBeamOffsetZInner(336),
327 fBeamSlopeZInner(336),
328 fInverseSlopeZ(kTRUE),
330 fFixedFitAside0(0.0),
331 fFixedFitAside1(1.0),
332 fFixedFitCside0(0.0),
340 fTracksEsdParam.SetOwner(kTRUE);
343 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
344 AliTPCcalibBase(calibLaser),
350 fTracksEsdParam(336),
352 fFullCalib(calibLaser.fFullCalib),
354 fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
355 fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
356 fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
357 fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
358 fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
359 fSignals(((calibLaser.fSignals))), // array of dedx signals
362 fHisLaser(0), // N dim histogram of laser
363 fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
364 fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
365 fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
366 fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
367 fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
368 fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
369 fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
370 fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
373 fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
374 fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
375 fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
376 fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
377 fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
378 fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
379 fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
380 fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
381 fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
382 fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
385 fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
386 fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
387 fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
388 fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
389 fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
390 fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
391 fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
392 fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
393 fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
396 fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
397 fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
398 fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
399 fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
400 fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
401 fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
402 fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
403 fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
404 fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
407 fDeltaYres(((calibLaser.fDeltaYres))),
408 fDeltaZres(((calibLaser.fDeltaZres))),
409 fDeltaYres2(((calibLaser.fDeltaYres))),
410 fDeltaZres2(((calibLaser.fDeltaZres))),
411 fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
412 fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
413 fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
414 fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
415 // fDeltaYres3(((calibLaser.fDeltaYres))),
416 //fDeltaZres3(((calibLaser.fDeltaZres))),
417 fFitAside(new TVectorD(5)), // drift fit - A side
418 fFitCside(new TVectorD(5)), // drift fit - C- side
419 fFitACside(new TVectorD(6)), // drift fit - C- side
420 fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
421 fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
422 fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
423 fNcuts(0), // number of cuts
424 fBeamSectorOuter(336),
425 fBeamSectorInner(336),
426 fBeamOffsetYOuter(336),
427 fBeamSlopeYOuter(336),
428 fBeamOffsetYInner(336),
429 fBeamSlopeYInner(336),
430 fBeamOffsetZOuter(336),
431 fBeamSlopeZOuter(336),
432 fBeamOffsetZInner(336),
433 fBeamSlopeZInner(336),
434 fInverseSlopeZ(calibLaser.fInverseSlopeZ),
435 fUseFixedDriftV(calibLaser.fUseFixedDriftV),
436 fFixedFitAside0(calibLaser.fFixedFitAside0),
437 fFixedFitAside1(calibLaser.fFixedFitAside1),
438 fFixedFitCside0(calibLaser.fFixedFitCside0),
439 fFixedFitCside1(calibLaser.fFixedFitCside1)
448 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
450 // assgnment operator
452 if (this != &calibLaser) {
453 new (this) AliTPCcalibLaser(calibLaser);
462 AliTPCcalibLaser::~AliTPCcalibLaser() {
467 delete fHisLaser; //->
468 delete fHisNclIn; //->Number of clusters inner
469 delete fHisNclOut; //->Number of clusters outer
470 delete fHisNclIO; //->Number of cluster inner outer
471 delete fHisLclIn; //->Level arm inner
472 delete fHisLclOut; //->Level arm outer
473 delete fHisLclIO; //->Number of cluster inner outer
478 delete fHisChi2YIn1; //->chi2 y inner - line
479 delete fHisChi2YOut1; //->chi2 y inner - line
480 delete fHisChi2YIn2; //->chi2 y inner - parabola
481 delete fHisChi2YOut2; //->chi2 y inner - parabola
482 delete fHisChi2YIO1; //->chi2 y IO - common
483 delete fHisChi2ZIn1; //->chi2 z inner - line
484 delete fHisChi2ZOut1; //->chi2 z inner - line
485 delete fHisChi2ZIn2; //->chi2 z inner - parabola
486 delete fHisChi2ZOut2; //->chi2 z inner - parabola
487 delete fHisChi2ZIO1; //->chi2 z IO - common
490 delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
491 delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
492 delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
493 delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
494 delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
495 delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
496 delete fHisPy2vP2In; //-> Curv P2inner - parabola
497 delete fHisPy2vP2Out; //-> Curv P2outer - parabola
498 delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
501 delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
502 delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
503 delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
504 delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
505 delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
506 delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
507 delete fHisPz2vP2In; //-> Curv P2inner - parabola
508 delete fHisPz2vP2Out; //-> Curv P2outer - parabola
509 delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
514 fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
515 fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
516 fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
517 fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
518 fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
519 fSignals.SetOwner(); //->Array of dedx signals
521 fDeltaZ.Delete(); //-> array of histograms of delta z for each track
522 fDeltaP3.Delete(); //-> array of histograms of P3 for each track
523 fDeltaP4.Delete(); //-> array of histograms of P4 for each track
524 fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
525 fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
526 fSignals.Delete(); //->Array of dedx signals
528 fDeltaYres.SetOwner();
530 delete fHisYAbsErrors;
531 fDeltaZres.SetOwner();
533 delete fHisZAbsErrors;
534 fDeltaYres2.SetOwner();
535 fDeltaYres2.Delete();
536 fDeltaZres2.SetOwner();
537 fDeltaZres2.Delete();
539 fDeltaYresAbs.SetOwner();
540 fDeltaYresAbs.Delete();
541 fDeltaZresAbs.SetOwner();
542 fDeltaZresAbs.Delete();
547 void AliTPCcalibLaser::Process(AliESDEvent * event) {
550 // Loop over tracks and call Process function
552 const Int_t kMinTracks=20;
553 const Int_t kMinClusters=40;
559 fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
563 if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
564 AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
566 // find CE background if present
568 if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
569 TH1D hisCE("hhisCE","hhisCE",100,-100,100);
570 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
571 AliESDtrack *track=fESD->GetTrack(i);
572 if (!track) continue;
573 hisCE.Fill(track->GetZ());
574 hisCE.Fill(track->GetZ()+2);
575 hisCE.Fill(track->GetZ()-2);
583 fTracksEsdParam.Delete();
584 for (Int_t id=0; id<336;id++) {
586 fClusterCounter[id]=0;
590 Int_t n=fESD->GetNumberOfTracks();
592 for (Int_t i=0;i<n;++i) {
593 AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
594 if (!friendTrack) continue;
595 AliESDtrack *track=fESD->GetTrack(i);
596 Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
597 if (binC>336) continue; //remove CE background
598 TObject *calibObject=0;
600 for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
601 if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
603 if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
605 Int_t id = FindMirror(track,seed);
606 if (id>=0) counter++;
611 if (counter<kMinTracks) return;
615 if (!fFullCalib) return;
616 static Bool_t init=kFALSE;
618 init = kTRUE; // way around for PROOF - to be investigated
622 for (Int_t id=0; id<336; id++){
624 if (!fTracksEsdParam.At(id)) continue;
625 if (fClusterSatur[id]>0.3) continue; // tracks in saturation
627 if ( AcceptLaser(id) && fFitZ[id]<0.5){
635 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
639 // for (Int_t id=0; id<336; id++){
642 if (!fTracksEsdParam.At(id)) return;
643 if (!AcceptLaser(id)) return;
644 Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
647 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
648 if (!hisdz) MakeFitHistos();
649 hisdz = (TH1F*)fDeltaZ.At(id);
650 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
651 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
652 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
653 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
654 TH1F * hisSignal = (TH1F*)fSignals.At(id);
657 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
658 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
659 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
668 param->GetPxPyPz(pxyz);
670 ltrp->GetPxPyPz(lpxyz);
672 Float_t dz = param->GetZ()-ltrp->GetZ();
673 Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
674 Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
675 if (hisdz) hisdz->Fill(dz);
676 if (hisP3) hisP3->Fill(param->GetParameter()[3]);
677 if (hisP4) hisP4->Fill(param->GetParameter()[4]);
678 if (hisdphi) hisdphi->Fill(dphi);
679 if (hisdphiP) hisdphiP->Fill(dphiP);
680 if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
682 xhis[0] = ltrp->GetId();
683 xhis[1] = ltrp->GetSide();
684 xhis[2] = ltrp->GetRod();
685 xhis[3] = ltrp->GetBundle();
686 xhis[4] = ltrp->GetBeam();
689 xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
690 xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
691 xhis[9] = param->GetParameter()[4];
692 xhis[10]= track->GetTPCNcls();
693 xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
695 fHisLaser->Fill(xhis);
698 void AliTPCcalibLaser::FitDriftV(){
700 // Fit corrections to the drift velocity - linear approximation in the z and global y
701 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
706 z = s* (z0 - vd*(t-t0))
710 vd - nominal drift velocity
711 zs - miscalibrated position
713 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
714 vr - relative change of the drift velocity
719 zs ~ z - s*vr*(z0-s*z)+s*dzt
720 --------------------------------
721 1. Correction function vr constant:
724 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
725 dzs/dl = dz/dl +s*s*vr*dz/dl
729 const Float_t kZCut = 240; // remove the closest laser beam
730 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
731 const Float_t kDistCut = 3; // distance sigma cut
732 const Float_t kDistCutAbs = 0.25;
733 const Float_t kMinClusters = 60; // minimal amount of the clusters
734 const Float_t kMinSignal = 16; // minimal mean height of the signal
735 const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
736 static TLinearFitter fdriftA(3,"hyp2");
737 static TLinearFitter fdriftC(3,"hyp2");
738 static TLinearFitter fdriftAC(4,"hyp3");
739 TVectorD fitA(3),fitC(3),fitAC(4);
741 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
742 AliTPCParam * tpcparam = calib->GetParameters();
744 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
754 for (Int_t iter=0; iter<3; iter++){
755 fdriftA.ClearPoints();
756 fdriftC.ClearPoints();
757 fdriftAC.ClearPoints();
759 for (Int_t id=0; id<336; id++){
760 if (!fTracksEsdParam.At(id)) continue;
761 if (!AcceptLaser(id)) continue;
762 if ( fClusterSatur[id]>kSaturCut) continue;
763 if ( fClusterCounter[id]<kMinClusters) continue;
764 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
765 if (track->GetTPCsignal()<kMinSignal) continue;
766 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
767 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
774 param->GetPxPyPz(pxyz);
776 ltrp->GetPxPyPz(lpxyz);
777 if (TMath::Abs(lxyz[2])>kZCut) continue;
778 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
779 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
780 if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
781 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
784 Double_t zlength = tpcparam->GetZLength(0);
785 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
786 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
787 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
788 if (ltrp->GetSide()==0){
789 fdriftA.AddPoint(xxx,mdrift,1);
791 fdriftC.AddPoint(xxx,mdrift,1);
793 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
794 fdriftAC.AddPoint(xxx2,mdrift,1);
797 if (fdriftA.GetNpoints()>10){
800 if (iter==0) fdriftA.EvalRobust(0.7);
801 else fdriftA.EvalRobust(0.8);
802 fdriftA.GetParameters(fitA);
803 npointsA= fdriftA.GetNpoints();
804 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
805 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
806 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
807 (*fFitAside)[0] = fitA[0];
808 (*fFitAside)[1] = fitA[1];
809 (*fFitAside)[2] = fitA[2];
810 (*fFitAside)[3] = fdriftA.GetNpoints();
811 (*fFitAside)[4] = chi2A;
814 if (fdriftC.GetNpoints()>10){
816 if (iter==0) fdriftC.EvalRobust(0.7);
817 else fdriftC.EvalRobust(0.8);
819 fdriftC.GetParameters(fitC);
820 npointsC= fdriftC.GetNpoints();
821 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
822 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
823 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
824 (*fFitCside)[0] = fitC[0];
825 (*fFitCside)[1] = fitC[1];
826 (*fFitCside)[2] = fitC[2];
827 (*fFitCside)[3] = fdriftC.GetNpoints();
828 (*fFitCside)[4] = chi2C;
832 if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
834 if (iter==0) fdriftAC.EvalRobust(0.7);
835 else fdriftAC.EvalRobust(0.8);
837 fdriftAC.GetParameters(fitAC);
838 npointsAC= fdriftAC.GetNpoints();
839 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
840 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
843 for (Int_t id=0; id<336; id++){
844 if (!fTracksEsdParam.At(id)) continue;
846 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
847 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
853 param->GetPxPyPz(pxyz);
855 ltrp->GetPxPyPz(lpxyz);
856 Double_t zlength = tpcparam->GetZLength(0);
857 Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
858 Double_t mdrift = zlength-TMath::Abs(xyz[2]);
861 if (ltrp->GetSide()==0){
862 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
864 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
867 fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
872 TTreeSRedirector *cstream = GetDebugStreamer();
873 TTimeStamp tstamp(fTime);
874 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
875 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
876 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
877 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
878 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
879 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
880 TVectorD vecGoofie(20);
881 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
883 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
884 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
885 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
889 (*cstream)<<"driftv"<<
890 "run="<<fRun<< // run number
891 "event="<<fEvent<< // event number
892 "time="<<fTime<< // time stamp of event
893 "trigger="<<fTrigger<< // trigger
894 "mag="<<fMagF<< // magnetic field
895 // Environment values
896 "press0="<<valuePressure0<<
897 "press1="<<valuePressure1<<
898 "pt0="<<ptrelative0<<
899 "pt1="<<ptrelative1<<
902 "vecGoofie.="<<&vecGoofie<<
906 "driftA.="<<fFitAside<<
907 "driftC.="<<fFitCside<<
908 "driftAC.="<<fFitACside<<
920 Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
922 // Fit corrections to the drift velocity - linear approximation in the z and global y
923 //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
925 // Source of outlyers :
926 // 0. Track in the saturation - postpeak
927 // 1. gating grid close the part of the signal for first bundle
929 // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
930 // 1. Robust fit is used in the itteration number 0
931 // only fraction of laser uted
932 // 2. Only the tracks close to the fit used in the second itteration
936 z = s* (z0 - vd*(t-t0))
940 vd - nominal drift velocity
941 zs - miscalibrated position
943 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
944 vr - relative change of the drift velocity
949 zs ~ z - s*vr*(z0-s*z)+s*dzt
950 --------------------------------
951 1. Correction function vr constant:
954 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
955 dzs/dl = dz/dl +s*s*vr*dz/dl
958 const Int_t knLaser = 336; //n laser tracks
959 const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
961 const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
962 const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
963 const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
964 const Float_t kMinClusters = 60.; // minimal amount of the clusters
965 const Float_t kMinSignal = 10.; // minimal mean height of the signal
966 const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
968 static TLinearFitter fdriftA(3,"hyp2");
969 static TLinearFitter fdriftC(3,"hyp2");
970 static TLinearFitter fdriftAC(4,"hyp3");
971 TVectorD fitA(3),fitC(3),fitAC(4);
973 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
974 AliTPCParam * tpcparam = calib->GetParameters();
978 for (Int_t id=0; id<336; id++) fFitZ[id]=0;
979 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
980 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
981 for (Int_t i=0;i<5; i++){
993 Int_t nbA[4]={0,0,0,0};
994 Int_t nbC[4]={0,0,0,0};
995 TVectorD vecZM(336); // measured z potion of laser
996 TVectorD vecA(336); // accepted laser
997 TVectorD vecZF(336); // fitted position
998 TVectorD vecDz(336); // deltaZ
999 TVectorD vecZS(336); // surveyed position of laser
1000 // additional variable to cut
1001 TVectorD vecdEdx(336); // dEdx
1002 TVectorD vecSy(336); // shape y
1003 TVectorD vecSz(336); // shape z
1006 for (Int_t id=0; id<336; id++){
1008 AliTPCLaserTrack *ltrp =
1009 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1010 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1011 AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1012 vecZM(id)= (param==0) ? 0:param->GetZ();
1013 vecZS(id)= ltrp->GetZ();
1016 vecdEdx(id)=(seed)?seed->GetdEdx():0;
1017 vecSy(id) =(seed)?seed->CookShape(1):0;
1018 vecSz(id) =(seed)?seed->CookShape(2):0;
1021 if (!fTracksEsdParam.At(id)) reject|=1;
1022 if (!param) continue;
1023 if (!AcceptLaser(id)) reject|=2;
1024 if ( fClusterSatur[id]>kSaturCut) reject|=4;
1025 if ( fClusterCounter[id]<kMinClusters) reject|=8;
1027 if (reject>0) continue;
1028 if (ltrp->GetSide()==0){
1030 nbA[ltrp->GetBundle()]++;
1032 if (ltrp->GetSide()==1){
1034 nbC[ltrp->GetBundle()]++;
1038 // reject "bad events"
1043 if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1045 for (Int_t i=0;i<4;i++){
1046 //count accepted for all layers
1047 if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1048 if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1050 if (naA<3 &&naC<3) isOK=kFALSE;
1051 if (isOK==kFALSE) return kFALSE;
1056 for (Int_t iter=0; iter<2; iter++){
1057 fdriftA.ClearPoints();
1058 fdriftC.ClearPoints();
1059 fdriftAC.ClearPoints();
1060 npointsA=0; npointsC=0; npointsAC=0;
1062 for (Int_t id=0; id<336; id++){
1063 if (!fTracksEsdParam.At(id)) continue;
1064 if (!AcceptLaser(id)) continue;
1065 if ( fClusterSatur[id]>kSaturCut) continue;
1066 if ( fClusterCounter[id]<kMinClusters) continue;
1067 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1068 if (track->GetTPCsignal()<kMinSignal) continue;
1069 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1070 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1076 param->GetPxPyPz(pxyz);
1078 ltrp->GetPxPyPz(lpxyz);
1079 Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1080 if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1081 if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1082 if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1084 // // drift distance
1085 // Double_t zlength = tpcparam->GetZLength(0);
1086 // Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1087 // Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1089 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1090 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1091 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1094 Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1095 if (iter==0 &<rp->GetBundle()==0) continue;
1096 // skip bundle 0 - can be wrong in the 0 iteration
1097 if (ltrp->GetSide()==0){
1098 fdriftA.AddPoint(xxx,mdrift,1);
1100 fdriftC.AddPoint(xxx,mdrift,1);
1102 Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1103 fdriftAC.AddPoint(xxx2,mdrift,1);
1106 if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1109 npointsA= fdriftA.GetNpoints();
1110 chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1111 fdriftA.EvalRobust(kFraction[iter]);
1112 fdriftA.GetParameters(fitA);
1113 if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1114 if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1115 (*fFitAside)[0] = fitA[0];
1116 (*fFitAside)[1] = fitA[1];
1117 (*fFitAside)[2] = fitA[2];
1118 (*fFitAside)[3] = fdriftA.GetNpoints();
1119 (*fFitAside)[4] = chi2A;
1122 if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1124 npointsC= fdriftC.GetNpoints();
1125 chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1126 fdriftC.EvalRobust(kFraction[iter]);
1127 fdriftC.GetParameters(fitC);
1128 if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1129 if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1130 (*fFitCside)[0] = fitC[0];
1131 (*fFitCside)[1] = fitC[1];
1132 (*fFitCside)[2] = fitC[2];
1133 (*fFitCside)[3] = fdriftC.GetNpoints();
1134 (*fFitCside)[4] = chi2C;
1138 if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1140 npointsAC= fdriftAC.GetNpoints();
1141 chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1142 fdriftAC.EvalRobust(kFraction[iter]);
1143 fdriftAC.GetParameters(fitAC);
1144 if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1145 (*fFitACside)[0] = fitAC[0];
1146 (*fFitACside)[1] = fitAC[1];
1147 (*fFitACside)[2] = fitAC[2];
1148 (*fFitACside)[3] = fdriftAC.GetNpoints();
1149 (*fFitACside)[4] = chi2AC;
1152 for (Int_t id=0; id<336; id++){
1153 if (!fTracksEsdParam.At(id)) continue;
1155 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1156 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1162 param->GetPxPyPz(pxyz);
1164 ltrp->GetPxPyPz(lpxyz);
1165 //Double_t zlength = tpcparam->GetZLength(0);
1166 //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1167 //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1168 Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1169 Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1170 Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1174 if (ltrp->GetSide()==0){
1175 fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1177 fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1180 //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1182 fFitZ[id]=mdrift-fz;
1184 vecDz[id]=mdrift-fz;
1186 if (fStreamLevel>0){
1187 TTreeSRedirector *cstream = GetDebugStreamer();
1188 TTimeStamp tstamp(fTime);
1189 Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1190 Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1191 Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1192 Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1193 Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1194 Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1195 TVectorD vecGoofie(20);
1196 AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1198 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1199 AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1200 if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1204 (*cstream)<<"driftvN"<<
1205 "run="<<fRun<< // run number
1206 "event="<<fEvent<< // event number
1207 "time="<<fTime<< // time stamp of event
1208 "trigger="<<fTrigger<< // trigger
1209 "mag="<<fMagF<< // magnetic field
1210 // Environment values
1211 "press0="<<valuePressure0<<
1212 "press1="<<valuePressure1<<
1213 "pt0="<<ptrelative0<<
1214 "pt1="<<ptrelative1<<
1217 "vecGoofie.="<<&vecGoofie<<
1220 "vecZM.="<<&vecZM<< // measured z position
1221 "vecZS.="<<&vecZS<< // surveyed z position
1222 "vecZF.="<<&vecZF<< // fitted z position
1223 "vecDz.="<<&vecDz<< // fitted z position
1224 "vecA.="<<&vecA<< // accept laser flag
1225 "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
1226 "vecSy.="<<&vecSy<< // shape y - to cut on
1227 "vecSz.="<<&vecSz<< // shape z - to cut on
1230 "driftA.="<<fFitAside<<
1231 "driftC.="<<fFitCside<<
1232 "driftAC.="<<fFitACside<<
1242 variables to check in debug mode:
1244 chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1245 chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1246 chainDriftN->SetAlias("driftF","vecZF.fElements");
1247 chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
1248 TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1249 TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1258 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1260 // get distance between mirror and laser track
1271 dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1274 // apply drift correction if already exist
1277 if (ltrp->GetSide()==0){
1278 if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1280 if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1282 if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1283 dist+=TMath::Abs(dz);
1285 // phi dist - divergence on 50 cm
1287 dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1292 Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
1297 TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1298 TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1299 TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1300 TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1302 TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1304 AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1305 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1306 //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1311 if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1312 if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
1313 if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
1314 if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
1321 Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1323 // Find corresponding mirror
1324 // add the corresponding tracks
1326 if (!track->GetOuterParam()) return -1;
1328 Float_t kRadius0 = 252;
1329 Float_t kRadius = 254.2;
1331 Float_t counterSatur=0;
1334 for (Int_t irow=158;irow>-1;--irow) {
1335 AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1337 Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1338 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1339 if (pedgeY<3) continue;
1340 if (pedgeX<3) continue;
1342 if (c->GetDetector()%36<18) csideA++;
1343 if (c->GetDetector()%36>=18) csideC++;
1344 if (c->GetMax()>900) counterSatur++;
1346 counterSatur/=(countercl+1);
1350 if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
1353 if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1356 AliExternalTrackParam param(*(track->GetOuterParam()));
1357 AliTracker::PropagateTrackTo(¶m,kRadius0,0.10566,3,kTRUE);
1358 AliTracker::PropagateTrackTo(¶m,kRadius,0.10566,0.1,kTRUE);
1359 AliTPCLaserTrack ltr;
1360 AliTPCLaserTrack *ltrp=0x0;
1361 // AliTPCLaserTrack *ltrpjw=0x0;
1363 Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
1364 // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
1365 //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1367 if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1368 ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1372 if (id<0) return -1;
1373 if (ltrp->GetSide()!=side) return -1;
1379 if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1382 Float_t radius=TMath::Abs(ltrp->GetX());
1383 param.Rotate(ltrp->GetAlpha());
1384 AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
1386 if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1387 Bool_t accept=kTRUE;
1389 // choose closer track
1391 AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1393 Float_t dist0=GetDistance(param0,ltrp);
1394 Float_t dist1=GetDistance(¶m,ltrp);
1395 if (dist0<dist1) accept=kFALSE;
1399 fClusterCounter[id]=countercl;
1400 fTracksEsdParam.AddAt(param.Clone(),id);
1401 fTracksEsd.AddAt(track,id);
1402 fTracksTPC.AddAt(seed,id);
1409 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1411 // Dump Laser info to the tree
1413 AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1414 AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1415 AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1424 param->GetPxPyPz(pxyz);
1426 ltrp->GetPxPyPz(lpxyz);
1427 Float_t dist3D = GetDistance(param,ltrp);
1428 Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1429 Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1432 if (fStreamLevel>0){
1433 TTreeSRedirector *cstream = GetDebugStreamer();
1434 Int_t time = fESD->GetTimeStamp();
1435 Bool_t accept = AcceptLaser(id);
1437 (*cstream)<<"Track"<<
1439 "run="<<fRun<< // run number
1440 "event="<<fEvent<< // event number
1441 "time="<<fTime<< // time stamp of event
1442 "trigger="<<fTrigger<< // trigger
1443 "mag="<<fMagF<< // magnetic field
1447 "driftA.="<<fFitAside<<
1448 "driftC.="<<fFitCside<<
1452 "distphi="<<distphi<<
1455 "counter="<<fCounter[id]<<
1456 "clcounter="<<fClusterCounter[id]<<
1457 "satur="<<fClusterSatur[id]<<
1481 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1483 // Refit the track with different tracklet models:
1484 // 1. Per ROC using the kalman filter, different edge cuts
1485 // 2. Per ROC linear in y and z
1486 // 3. Per ROC quadratic in y and z
1487 // 4. Per track offset for each sector, linear for each sector, common quadratic
1488 // store x, y, z information for all models and the cluster to calculate the residuals
1491 // The clusters which do not fulfill given criteria are skipped
1493 // Cluters removed from fit
1494 // 1. Extended shape > kShapeCut
1495 // 2. In saturation Max > kMax
1496 // 3. Distance to edge < cutEdge
1498 // Clusters not used for the calibration:
1500 // 1. Extended shape > kShapeCut
1501 // 2. In saturation Max > kMax
1504 AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1505 AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1506 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1508 AliTPCclusterMI dummyCl;
1511 Int_t kMaxTracklets=2;
1512 Float_t kShapeCut = 1.3;
1513 Float_t kRatioCut = 2.;
1515 Float_t kMax = 900.;
1518 //=============================================//
1519 // Linear Fitters for the Different Approaches //
1520 //=============================================//
1521 //linear fit model in y and z; inner - outer sector, combined with offset
1522 static TLinearFitter fy1I(2,"hyp1");
1523 static TLinearFitter fy1O(2,"hyp1");
1524 static TLinearFitter fz1I(2,"hyp1");
1525 static TLinearFitter fz1O(2,"hyp1");
1526 static TLinearFitter fy1IO(3,"hyp2");
1527 static TLinearFitter fz1IO(3,"hyp2");
1528 //quadratic fit model in y and z; inner - sector
1529 static TLinearFitter fy2I(3,"hyp2");
1530 static TLinearFitter fy2O(3,"hyp2");
1531 static TLinearFitter fz2I(3,"hyp2");
1532 static TLinearFitter fz2O(3,"hyp2");
1533 //common quadratic fit for IROC and OROC in y and z
1534 static TLinearFitter fy4(5,"hyp4");
1535 static TLinearFitter fz4(5,"hyp4");
1545 //=============================//
1546 // Loop over all Tracklet Cuts //
1547 //=============================//
1548 for (Int_t icut=0; icut<fNcuts; icut++){
1549 Float_t xinMin = 2500, xinMax=-90;
1550 Float_t xoutMin=2500, xoutMax=-90;
1551 Float_t msigmaYIn=0;
1552 Float_t msigmaYOut=0;
1553 Float_t mqratioIn=0;
1554 Float_t mqratioOut=0;
1557 AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1559 Double_t edgeCutX = fEdgeXcuts[icut];
1560 Double_t edgeCutY = fEdgeYcuts[icut];
1561 Int_t nclCut = (Int_t)fNClCuts[icut];
1562 //=========================//
1563 // Parameters to calculate //
1564 //=========================//
1565 //fit parameter inner, outer tracklet and combined fit
1566 TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1567 TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1569 TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1570 TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1571 TVectorD vecy4res(5),vecz4res(5);
1572 TVectorD vecy1resIO(3),vecz1resIO(3);
1573 // cluster and track positions for each row - used for residuals
1574 TVectorD vecgX(159); // global X
1575 TVectorD vecgY(159); // global Y
1576 TVectorD vecgZ(159); // global Z
1578 TVectorD vecX(159); // x is the same for all (row center)
1579 TVectorD vecYkalman(159); // y from kalman fit
1580 TVectorD vecZkalman(159); // z from kalman fit
1581 TVectorD vecY1(159); // y from pol1 fit per ROC
1582 TVectorD vecZ1(159); // z from pol1 fit per ROC
1583 TVectorD vecY1IO(159); // y from pol1 fit per ROC
1584 TVectorD vecZ1IO(159); // z from pol1 fit per ROC
1585 TVectorD vecY2(159); // y from pol2 fit per ROC
1586 TVectorD vecZ2(159); // z from pol2 fit per ROC
1587 TVectorD vecY4(159); // y from sector fit
1588 TVectorD vecZ4(159); // z from sector fit
1589 TVectorD vecClY(159); // y cluster position
1590 TVectorD vecClZ(159); // z cluster position
1591 TVectorD vecSec(159); // sector for each row
1592 TVectorD isReject(159); // flag - cluster to be rejected
1594 Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
1595 Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
1596 Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
1597 Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
1598 Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
1599 Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
1600 Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
1601 Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
1602 Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
1603 Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
1604 Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
1605 Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
1607 Int_t innerSector = -1; // number of inner sector
1608 Int_t outerSector = -1; // number of outer sector
1609 Int_t nclI=0; // number of clusters (inner)
1610 Int_t nclO=0; // number of clusters (outer)
1611 //=================================================//
1612 // Perform the Kalman Fit using the Tracklet Class //
1613 //=================================================//
1614 AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1615 TObjArray tracklets=
1616 AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1617 kFALSE,nclCut,kMaxTracklets);
1618 // tracklet pointers
1619 AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1620 AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1621 AliTPCTracklet *tr=0x0;
1622 AliTPCTracklet dummy;
1623 //continue if we didn't find a tracklet
1624 if ( !trInner && !trOuter ) continue;
1625 //================================================================================//
1626 // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1627 //================================================================================//
1628 if ( trInner && trOuter ){
1629 if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1630 if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1634 AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1638 if ( !trInner->GetInner() ) continue;
1640 if ( trInner->GetSector()>35 ){
1645 if ( !trOuter->GetInner() ) continue;
1647 if ( trOuter->GetSector()<36 ){
1653 innerSector = trInner->GetSector();
1654 if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1655 outerSector = trOuter->GetSector();
1656 if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1658 // array of clusters
1659 TClonesArray arrCl("AliTPCclusterMI",159);
1660 arrCl.ExpandCreateFast(159);
1661 //=======================================//
1662 // fill fitters with cluster information //
1663 //=======================================//
1664 AliDebug(3,"Filing Cluster Information");
1669 for (Int_t irow=158;irow>-1;--irow) {
1670 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1671 AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1676 Float_t meanY=0, sumY=0;
1677 for (Int_t drow=-1;drow<=1;drow++) {
1678 if (irow+drow<0) continue;
1679 if (irow+drow>158) continue;
1680 AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1681 if (!ccurrent) continue;
1682 Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1683 if ( roc!=innerSector && roc!=outerSector ) continue;
1684 if (ccurrent->GetX()<10) continue;
1685 if (ccurrent->GetY()==0) continue;
1686 meanY+=ccurrent->GetY();
1689 if (sumY>0) meanY/=sumY;
1694 Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1695 Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1698 Int_t roc = static_cast<Int_t>(c->GetDetector());
1699 if ( roc!=innerSector && roc!=outerSector ) continue;
1701 //store clusters in clones array
1704 if (c->GetMax()<4) continue; // noise cluster?
1705 if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
1707 vecX[irow] = c->GetX();
1708 vecClY[irow] = c->GetY();
1709 vecClZ[irow] = c->GetZ();
1712 // c->GetGlobalXYZ(gxyz);
1713 // vecgX[irow] = gxyz[0];
1714 // vecgY[irow] = gxyz[1];
1715 // vecgZ[irow] = gxyz[2];
1717 Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1718 Double_t y = vecClY[irow];
1719 Double_t z = vecClZ[irow];
1721 Double_t x2[2]={x,x*x}; //linear and parabolic parameters
1722 Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1723 Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
1724 if ( roc == innerSector ) {
1725 x4[0]=1; //offset inner - outer sector
1726 x4[1]=x; //slope parameter inner sector
1729 x4[2]=x; //slope parameter outer sector
1731 x4[3]=x*x; //common parabolic shape
1732 if (pedgeX < fEdgeXcuts[icut]) continue;
1733 if (pedgeY < fEdgeYcuts[icut]) continue;
1734 if (c->GetMax()>900) continue; // cluster in saturation
1736 if ( roc==innerSector ){
1737 fy1I.AddPoint(x2,y);
1738 fz1I.AddPoint(x2,z);
1739 fy2I.AddPoint(x2,y);
1740 fz2I.AddPoint(x2,z);
1742 if (c->GetX()<xinMin) xinMin=c->GetX();
1743 if (c->GetX()>xinMax) xinMax=c->GetX();
1745 msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1746 mqratioIn +=c->GetMax()/c->GetQ();
1749 if ( roc==outerSector ){
1750 fy1O.AddPoint(x2,y);
1751 fz1O.AddPoint(x2,z);
1752 fy2O.AddPoint(x2,y);
1753 fz2O.AddPoint(x2,z);
1755 if (c->GetX()<xoutMin) xoutMin=c->GetX();
1756 if (c->GetX()>xoutMax) xoutMax=c->GetX();
1757 msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1758 mqratioOut +=c->GetMax()/c->GetQ();
1763 fy1IO.AddPoint(xIO,y);
1764 fz1IO.AddPoint(xIO,z);
1774 //======================================//
1775 // Evaluate and retrieve fit parameters //
1776 //======================================//
1777 AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1779 if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1784 fy1I.GetParameters(vecy1resInner);
1785 fz1I.GetParameters(vecz1resInner);
1786 fy2I.GetParameters(vecy2resInner);
1787 fz2I.GetParameters(vecz2resInner);
1788 chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1789 chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1790 chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1791 chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1794 if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1799 fy1O.GetParameters(vecy1resOuter);
1800 fz1O.GetParameters(vecz1resOuter);
1801 fy2O.GetParameters(vecy2resOuter);
1802 fz2O.GetParameters(vecz2resOuter);
1803 chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1804 chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1805 chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1806 chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1809 if ( innerSector>0 && outerSector>0 ){
1810 if (fy4.GetNpoints()>0) {
1812 fy4.GetParameters(vecy4res);
1813 chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1815 if (fz4.GetNpoints()>0) {
1817 fz4.GetParameters(vecz4res);
1818 chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1820 if (fy1IO.GetNpoints()>0) {
1822 fy1IO.GetParameters(vecy1resIO);
1823 chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1825 if (fz1IO.GetNpoints()>0) {
1827 fz1IO.GetParameters(vecz1resIO);
1828 chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1832 fy4.ClearPoints(); fz4.ClearPoints();
1833 fy1I.ClearPoints(); fy1O.ClearPoints();
1834 fz1I.ClearPoints(); fz1O.ClearPoints();
1835 fy2I.ClearPoints(); fy2O.ClearPoints();
1836 fz2I.ClearPoints(); fz2O.ClearPoints();
1837 fy1IO.ClearPoints(); fz1IO.ClearPoints();
1838 //==============================//
1839 // calculate tracklet positions //
1840 //==============================//
1841 AliDebug(4,"Calculate tracklet positions");
1842 for (Int_t irow=158;irow>-1;--irow) {
1844 AliTPCclusterMI *c=track->GetClusterPointer(irow);
1845 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1849 if (!c) { //no cluster
1852 if (c->GetMax()>kMax){ //saturation
1855 if ( vecSec[irow] == outerSector ) { //extended shape
1856 if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
1857 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
1859 if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
1860 if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
1866 if ( vecSec[irow]==-1 ) continue; //no cluster info
1867 if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1869 Double_t x=vecX[irow];
1870 Double_t xref=x-133.4;
1872 Double_t yoffInner=0;
1873 Double_t zoffInner=0;
1874 Double_t yoffInner1=0;
1875 Double_t zoffInner1=0;
1876 Double_t yslopeInner=0;
1877 Double_t yslopeOuter=0;
1878 Double_t zslopeInner=0;
1879 Double_t zslopeOuter=0;
1880 //positions of hyperplane fits
1881 if ( vecSec[irow] == outerSector ) {
1883 vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1884 vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1885 vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1886 vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1887 yslopeOuter=vecy4res[3];
1888 zslopeOuter=vecz4res[3];
1891 vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1892 vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1893 vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1894 vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1895 yoffInner=vecy4res[1];
1896 zoffInner=vecz4res[1];
1897 yoffInner1=vecy1resIO[1];
1898 zoffInner1=vecz1resIO[1];
1899 yslopeInner=vecy4res[2];
1900 zslopeInner=vecz4res[2];
1902 vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1903 vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1904 vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1905 vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1906 //positions of kalman fits
1907 Double_t gxyz[3],xyz[3];
1908 AliExternalTrackParam *param = 0x0;
1910 param=tr->GetInner();
1912 param->GetXYZ(gxyz);
1913 Float_t bz = AliTracker::GetBz(gxyz);
1914 param->GetYAt(x, bz, xyz[1]);
1915 param->GetZAt(x, bz, xyz[2]);
1916 vecYkalman[irow]=xyz[1];
1917 vecZkalman[irow]=xyz[2];
1924 //=====================================================================//
1925 // write results from the different tracklet fits with debug streamers //
1926 //=====================================================================//
1927 if (fStreamLevel>4){
1928 TTreeSRedirector *cstream = GetDebugStreamer();
1930 Float_t dedx = track->GetdEdx();
1931 (*cstream)<<"FitModels"<<
1932 "run="<<fRun<< // run number
1933 "event="<<fEvent<< // event number
1934 "time="<<fTime<< // time stamp of event
1935 "trigger="<<fTrigger<< // trigger
1936 "mag="<<fMagF<< // magnetic field
1939 "edgeCutX=" << edgeCutX <<
1940 "edgeCutY=" << edgeCutY <<
1941 "nclCut=" << nclCut <<
1942 "innerSector="<< innerSector <<
1943 "outerSector="<< outerSector <<
1946 "Tr.=" << extparam <<
1947 "yPol1In.=" << &vecy1resInner <<
1948 "zPol1In.=" << &vecz1resInner <<
1949 "yPol1InOut.="<< &vecy1resIO <<
1950 "zPol1InOut.="<< &vecz1resIO <<
1951 "yPol2In.=" << &vecy2resInner <<
1952 "zPol2In.=" << &vecz2resInner <<
1953 "yPol1Out.=" << &vecy1resOuter <<
1954 "zPol1Out.=" << &vecz1resOuter <<
1955 "yPol2Out.=" << &vecy2resOuter <<
1956 "zPol2Out.=" << &vecz2resOuter <<
1957 "yInOut.=" << &vecy4res <<
1958 "zInOut.=" << &vecz4res <<
1959 "chi2y1In=" << chi2I1y <<
1960 "chi2z1In=" << chi2I1z <<
1961 "chi2y1InOut="<< chi2IO1y <<
1962 "chi2z1InOut="<< chi2IO1z <<
1963 "chi2y1Out=" << chi2O1y <<
1964 "chi2z1Out=" << chi2O1z <<
1965 "chi2y2In=" << chi2I2y <<
1966 "chi2z2In=" << chi2I2z <<
1967 "chi2y2Out=" << chi2O2y <<
1968 "chi2z2Out=" << chi2O2z <<
1969 "chi2yInOut=" << chi2IOy <<
1970 "chi2zInOut=" << chi2IOz <<
1971 "trletIn.=" << trInner <<
1972 "trletOut.=" << trOuter <<
1975 "xinMin=" << xinMin<<
1976 "xinMax=" << xinMax<<
1977 "xoutMin=" << xoutMin<<
1978 "xoutMax=" << xoutMax<<
1979 "msigmaYIn=" <<msigmaYIn<<
1980 "msigmaYOut=" <<msigmaYOut<<
1981 "mqratioIn=" <<mqratioIn<<
1982 "mqratioOut=" << mqratioOut <<
1987 // wirte residuals information
1988 if (fStreamLevel>5){
1989 TTreeSRedirector *cstream = GetDebugStreamer();
1991 Float_t dedx = track->GetdEdx();
1992 (*cstream)<<"Residuals"<<
1993 "run="<<fRun<< // run number
1994 "event="<<fEvent<< // event number
1995 "time="<<fTime<< // time stamp of event
1996 "trigger="<<fTrigger<< // trigger
1997 "mag="<<fMagF<< // magnetic field
2000 "edgeCutX=" << edgeCutX <<
2001 "edgeCutY=" << edgeCutY <<
2002 "nclCut=" << nclCut <<
2004 "Tr.=" << extparam<<
2007 "vX.=" << &vecgX<< // global x
2008 "vY.=" << &vecgY<< // global y
2009 "vZ.=" << &vecgZ<< // global z
2011 "TrYpol1.=" << &vecY1 <<
2012 "TrZpol1.=" << &vecZ1 <<
2013 "TrYpol2.=" << &vecY2 <<
2014 "TrZpol2.=" << &vecZ2 <<
2015 "TrYpol1InOut.="<< &vecY1IO <<
2016 "TrZpol1InOut.="<< &vecZ1IO <<
2017 "TrYInOut.=" << &vecY4 <<
2018 "TrZInOut.=" << &vecZ4 <<
2019 "ClY.=" << &vecClY <<
2020 "ClZ.=" << &vecClZ <<
2021 "isReject.=" << &isReject<<
2022 "sec.=" << &vecSec <<
2025 "xinMin=" << xinMin<<
2026 "xinMax=" << xinMax<<
2027 "xoutMin=" << xoutMin<<
2028 "xoutMax=" << xoutMax<<
2029 "msigmaYIn=" <<msigmaYIn<<
2030 "msigmaYOut=" <<msigmaYOut<<
2031 "mqratioIn=" <<mqratioIn<<
2032 "mqratioOut=" << mqratioOut <<
2033 "yInOut.=" << &vecy4res <<
2034 "zInOut.=" << &vecz4res <<
2036 "chi2y1In=" << chi2I1y << //
2037 "chi2z1In=" << chi2I1z <<
2038 "chi2y1Out=" << chi2O1y <<
2039 "chi2z1Out=" << chi2O1z <<
2040 "chi2y1InOut="<< chi2IO1y <<
2041 "chi2z1InOut="<< chi2IO1z <<
2042 "chi2y2In=" << chi2I2y <<
2043 "chi2z2In=" << chi2I2z <<
2044 "chi2y2Out=" << chi2O2y <<
2045 "chi2z2Out=" << chi2O2z <<
2046 "chi2yInOut=" << chi2IOy <<
2047 "chi2zInOut=" << chi2IOz <<
2049 "yPol1In.=" << &vecy1resInner <<
2050 "zPol1In.=" << &vecz1resInner <<
2051 "yPol2In.=" << &vecy2resInner <<
2052 "zPol2In.=" << &vecz2resInner <<
2053 "yPol1Out.=" << &vecy1resOuter <<
2054 "zPol1Out.=" << &vecz1resOuter <<
2055 "yPol1InOut.="<< &vecy1resIO <<
2056 "zPol1InOut.="<< &vecz1resIO <<
2057 "yPol2Out.=" << &vecy2resOuter <<
2058 "zPol2Out.=" << &vecz2resOuter <<
2064 //==========================//
2065 // Fill Residual Histograms //
2066 //==========================//
2067 if (!fHisNclIn) MakeFitHistos();
2069 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2070 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2071 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2072 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2073 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2074 TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2075 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2076 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2078 for (Int_t irow=158;irow>-1;--irow) {
2079 if (vecSec[irow]==-1)continue; // no cluster info
2080 if (isReject[irow]>0.5) continue; //
2081 Double_t ycl = vecClY[irow];
2082 Double_t yfit = vecY1[irow];
2083 Double_t yfit2 = vecY2[irow];
2084 Double_t x = vecX[irow];
2085 Double_t yabsbeam = -1000;
2086 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2087 yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2088 else if(innerSector==fBeamSectorInner[id])
2089 yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2091 // Double_t yfit3 = vecY2[irow];
2092 Double_t zcl = vecClZ[irow];
2093 Double_t zfit = vecZ1[irow];
2094 Double_t zfit2 = vecZ2[irow];
2095 //Double_t zfit3 = vecZ2[irow];
2098 // The expressions for zcorrected has been obtained by
2099 // inverting the fits in the FitDriftV() method (ignoring the
2100 // global y dependence for now):
2102 // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2104 // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2107 // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
2109 // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2111 Double_t dzabs = -1000;
2112 Double_t zcorrected = -1000;
2113 if (ltrp->GetSide()==0){
2114 if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2115 // ignore global y dependence for now
2117 if(!fUseFixedDriftV)
2118 zcorrected = (zcl + (*fFitAside)[0] -
2119 (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2121 zcorrected = (zcl + fFixedFitAside0 -
2122 (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2123 // zcorrected = zcl;
2124 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2125 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2126 else if(innerSector==fBeamSectorInner[id])
2127 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2130 if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2132 if(!fUseFixedDriftV)
2133 zcorrected = (zcl - (*fFitCside)[0] +
2134 (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2136 zcorrected = (zcl - fFixedFitCside0 +
2137 (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2139 // zcorrected = zcl;
2140 if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2141 dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2142 else if(innerSector==fBeamSectorInner[id])
2143 dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2147 if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2149 profy->Fill(irow,ycl-yfit);
2150 profy2->Fill(irow,ycl-yfit2);
2152 fHisYAbsErrors->Fill(id);
2153 // profyabs->Fill(irow,-0.99);
2155 profyabs->Fill(irow,ycl-yabsbeam);
2157 // profy3->Fill(irow,ycl-yfit3);
2160 profz->Fill(irow,zcl-zfit);
2161 profz2->Fill(irow,zcl-zfit2);
2162 //profz3->Fill(irow,zcl-zfit3);
2165 fHisZAbsErrors->Fill(id);
2167 profzabs->Fill(irow,dzabs);
2173 // Fill laser fit histograms
2175 Float_t dedx = track->GetdEdx();
2176 if (nclI>20&&nclO>20){
2177 fHisNclIn->Fill(id,nclI); //->Number of clusters inner
2178 fHisNclOut->Fill(id,nclO); //->Number of clusters outer
2179 fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
2181 fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
2182 fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
2183 fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
2185 fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2186 fHisdZfit->Fill(id,fFitZ[id]);
2189 fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
2190 fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
2191 fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
2192 fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
2193 fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
2196 fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
2197 fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
2198 fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
2199 fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
2200 fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
2203 fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
2204 fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
2205 fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
2207 fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
2208 fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
2209 fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
2211 fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
2212 fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
2213 fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
2214 fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
2216 fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
2217 fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
2218 fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
2219 fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
2222 fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
2223 fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
2227 fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
2228 fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
2236 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2238 // Dump information about laser beams
2239 // isOK variable indicates usability of the beam
2240 // Beam is not usable if:
2241 // a. No entries in range (krmsCut0)
2242 // b. Big sperad (krmscut1)
2243 // c. RMSto fit sigma bigger then (kmultiCut)
2244 // d. Too big angular spread
2247 const Float_t krmsCut0=0.001;
2248 const Float_t krmsCut1=0.16;
2249 const Float_t kmultiCut=2;
2250 const Float_t kcutP0=0.002;
2252 AliTPCcalibLaser *laser = this;
2253 TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2254 TF1 fg("fg","gaus");
2255 AliTPCParam * tpcparam = 0;
2256 // start set up for absolute residuals analysis
2258 AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2259 tpcparam = calib->GetParameters();
2260 if (!tpcparam) tpcparam = new AliTPCParamSR;
2262 AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2264 Float_t bfield = 0, bz=0;
2266 current = grp->GetL3Current((AliGRPObject::Stats)0);
2267 bfield = 5*current/30000.;
2268 bz = 5*current/30000.;
2269 printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2272 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2273 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2274 TLinearFitter lfabsyInner(2);
2275 lfabsyInner.SetFormula("1 ++ x");
2276 TLinearFitter lfabszInner(2);
2277 lfabszInner.SetFormula("1 ++ x");
2279 TLinearFitter lfabsyOuter(2);
2280 lfabsyOuter.SetFormula("1 ++ x");
2281 TLinearFitter lfabszOuter(2);
2282 lfabszOuter.SetFormula("1 ++ x");
2283 // end set up for absolute residuals analysis
2287 for (Int_t id=0; id<336; id++){
2289 TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
2290 TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2291 TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
2292 TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
2293 TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
2294 TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2295 //if (!hisphi) continue;
2296 Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2297 //if (entries<minEntries) continue;
2299 AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2301 AliTPCLaserTrack::LoadTracks();
2302 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2304 pcstream->GetFile()->cd();
2305 if (hisphi) hisphi->Write();
2306 if (hisphiP) hisphiP->Write();
2307 if (hisZ) hisZ->Write();
2308 if (hisP3) hisP3->Write();
2309 if (hisP4) hisP4->Write();
2311 Float_t meanphi = hisphi->GetMean();
2312 Float_t rmsphi = hisphi->GetRMS();
2314 Float_t meanphiP = hisphiP->GetMean();
2315 Float_t rmsphiP = hisphiP->GetRMS();
2316 Float_t meanZ = hisZ->GetMean();
2317 Float_t rmsZ = hisZ->GetRMS();
2318 if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2319 hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2320 Double_t gphi1 = fg.GetParameter(1);
2321 Double_t gphi2 = fg.GetParameter(2);
2322 if (hisphiP->GetRMS()>0)
2323 hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2324 Double_t gphiP1 = fg.GetParameter(1);
2325 Double_t gphiP2 = fg.GetParameter(2);
2327 if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2328 hisZ->Fit(&fg,"","");
2329 Double_t gz1 = fg.GetParameter(1);
2330 Double_t gz2 = fg.GetParameter(2);
2332 if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2333 hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2334 Double_t gp31 = fg.GetParameter(1);
2335 Double_t gp32 = fg.GetParameter(2);
2337 if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2338 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2339 Double_t gp41 = fg.GetParameter(1);
2340 Double_t gp42 = fg.GetParameter(2);
2342 Float_t meanS=hisS->GetMean();
2347 ltrp->GetPxPyPz(lpxyz);
2349 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2350 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2351 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2352 if (gphiP2>kcutP0) isOK=kFALSE;
2358 his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2359 Float_t mnclIn = his->GetMean();
2361 his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2362 Float_t mnclOut = his->GetMean();
2364 his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2365 Float_t mnclIO = his->GetMean();
2367 his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2368 Float_t mLclIn = his->GetMean();
2370 his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2371 Float_t mLclOut = his->GetMean();
2373 his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2374 Float_t mLclIO = his->GetMean();
2377 his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2378 Float_t mdEdx = his->GetMean();
2384 his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2385 Float_t mChi2YIn1= his->GetMean();
2387 his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2388 Float_t mChi2YOut1 = his->GetMean();
2390 his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2391 Float_t mChi2YIn2 = his->GetMean();
2393 his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2394 Float_t mChi2YOut2 = his->GetMean();
2396 his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2397 Float_t mChi2YIO1 = his->GetMean();
2399 his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2400 Float_t mChi2ZIn1 = his->GetMean();
2402 his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2403 Float_t mChi2ZOut1 = his->GetMean();
2405 his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2406 Float_t mChi2ZIn2 = his->GetMean();
2408 his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2409 Float_t mChi2ZOut2 = his->GetMean();
2411 his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2412 Float_t mChi2ZIO1 = his->GetMean();
2417 his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2418 Float_t edZfit = his->GetEntries();
2419 Float_t mdZfit = his->GetMean();
2420 Float_t rdZfit = his->GetRMS();
2423 his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2424 Float_t ePy1vP0 = his->GetEntries();
2425 Float_t mPy1vP0 = his->GetMean();
2426 Float_t rPy1vP0 = his->GetRMS();
2429 his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2430 Float_t ePy2vP0 = his->GetEntries();
2431 Float_t mPy2vP0 = his->GetMean();
2432 Float_t rPy2vP0 = his->GetRMS();
2435 his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2436 Float_t ePy3vP0 = his->GetEntries();
2437 Float_t mPy3vP0 = his->GetMean();
2438 Float_t rPy3vP0 = his->GetRMS();
2441 his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2442 Float_t ePy1vP1 = his->GetEntries();
2443 Float_t mPy1vP1 = his->GetMean();
2444 Float_t rPy1vP1 = his->GetRMS();
2447 his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2448 Float_t ePy2vP1 = his->GetEntries();
2449 Float_t mPy2vP1 = his->GetMean();
2450 Float_t rPy2vP1 = his->GetRMS();
2453 his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2454 Float_t ePy3vP1 = his->GetEntries();
2455 Float_t mPy3vP1 = his->GetMean();
2456 Float_t rPy3vP1 = his->GetRMS();
2459 his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2460 Float_t ePy2vP2In = his->GetEntries();
2461 Float_t mPy2vP2In = his->GetMean();
2462 Float_t rPy2vP2In = his->GetRMS();
2465 his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2466 Float_t ePy2vP2Out = his->GetEntries();
2467 Float_t mPy2vP2Out = his->GetMean();
2468 Float_t rPy2vP2Out = his->GetRMS();
2471 his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2472 Float_t ePy3vP2IO = his->GetEntries();
2473 Float_t mPy3vP2IO = his->GetMean();
2474 Float_t rPy3vP2IO = his->GetRMS();
2479 his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2480 Float_t ePz1vP0 = his->GetEntries();
2481 Float_t mPz1vP0 = his->GetMean();
2482 Float_t rPz1vP0 = his->GetRMS();
2485 his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2486 Float_t ePz2vP0 = his->GetEntries();
2487 Float_t mPz2vP0 = his->GetMean();
2488 Float_t rPz2vP0 = his->GetRMS();
2491 his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2492 Float_t ePz3vP0 = his->GetEntries();
2493 Float_t mPz3vP0 = his->GetMean();
2494 Float_t rPz3vP0 = his->GetRMS();
2497 his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2498 Float_t ePz1vP1 = his->GetEntries();
2499 Float_t mPz1vP1 = his->GetMean();
2500 Float_t rPz1vP1 = his->GetRMS();
2503 his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2504 Float_t ePz2vP1 = his->GetEntries();
2505 Float_t mPz2vP1 = his->GetMean();
2506 Float_t rPz2vP1 = his->GetRMS();
2509 his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2510 Float_t ePz3vP1 = his->GetEntries();
2511 Float_t mPz3vP1 = his->GetMean();
2512 Float_t rPz3vP1 = his->GetRMS();
2515 his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2516 Float_t ePz2vP2In = his->GetEntries();
2517 Float_t mPz2vP2In = his->GetMean();
2518 Float_t rPz2vP2In = his->GetRMS();
2521 his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2522 Float_t ePz2vP2Out = his->GetEntries();
2523 Float_t mPz2vP2Out = his->GetMean();
2524 Float_t rPz2vP2Out = his->GetRMS();
2527 his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2528 Float_t ePz3vP2IO = his->GetEntries();
2529 Float_t mPz3vP2IO = his->GetMean();
2530 Float_t rPz3vP2IO = his->GetRMS();
2533 // Fit absolute laser residuals
2534 TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2535 TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2537 Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2538 Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2540 TVectorD vecX(159); // X
2541 TVectorD vecY(159); // Y
2542 TVectorD vecR(159); // R
2543 TVectorD vecDY(159); // absolute residuals in Y
2544 TVectorD vecDZ(159); // absolute residuals in Z
2545 TVectorD vecN(159); // number of clusters
2546 TVectorD vecEy(159); //error y
2547 TVectorD vecEz(159); //error z
2548 TVectorD vecPhi(159); // local tangent
2549 TVectorD vecPhiR(159); // local tangent
2551 lfabsyInner.ClearPoints();
2552 lfabszInner.ClearPoints();
2553 lfabsyOuter.ClearPoints();
2554 lfabszOuter.ClearPoints();
2557 Float_t yAbsInOffset = -100;
2558 Float_t yAbsInSlope = -100;
2559 Float_t yAbsInDeltay = -100;
2561 Float_t zAbsInOffset = -100;
2562 Float_t zAbsInSlope = -100;
2563 Float_t zAbsInDeltaz = -100;
2565 Float_t yAbsOutOffset = -100;
2566 Float_t yAbsOutSlope = -100;
2567 Float_t yAbsOutDeltay = -100;
2569 Float_t zAbsOutOffset = -100;
2570 Float_t zAbsOutSlope = -100;
2571 Float_t zAbsOutDeltaz = -100;
2573 Float_t lasTanPhiLocIn = -100;
2574 Float_t lasTanPhiLocOut = -100;
2576 if(histAbsY->GetEntries()>0) {
2578 Double_t rotAngOut = 10;
2579 Double_t rotAngIn = 10;
2580 if((secInner%36)!=(secOuter%36))
2581 rotAngIn += 20; // 30 degrees
2583 // Calculate laser mirror X position in local frame
2584 Double_t laserposOut =
2585 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2586 Double_t laserposIn =
2587 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2589 // Calculate laser tan phi in local frame
2590 lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2591 if(lasTanPhiLocOut<0) {
2592 lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2593 lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2596 lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2597 lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2600 lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2601 lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2603 TProfile* yprof = histAbsY->ProfileX("yprof");
2604 TProfile* zprof = histAbsZ->ProfileX("zprof");
2606 for(Int_t bin = 1; bin<=159; bin++) {
2608 if(yprof->GetBinEntries(bin)<10&&
2609 zprof->GetBinEntries(bin)<10) {
2613 // There is a problem in defining inner and outer sectors for
2614 // the outer beams (0 and 6) where both sectors are OROCs. To
2615 // make sure there is no overlap row 94 to 99 are cutted.
2616 if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2619 Int_t row = (bin-1);
2623 Bool_t isOuter = kTRUE;
2624 Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2627 (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2630 sector = TMath::Nint(fBeamSectorInner[id]);
2634 Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2635 vecN[bin-1] =yprof->GetBinEntries(bin);
2636 vecEy[bin-1]=yprof->GetBinError(bin);
2637 vecEz[bin-1]=zprof->GetBinError(bin);
2639 vecDY[bin-1] = yprof->GetBinContent(bin);
2640 vecDZ[bin-1] = zprof->GetBinContent(bin);
2641 if(!isOuter) { // inner
2642 vecPhi[bin-1]=lasTanPhiLocIn;
2643 // Calculate local y from residual and database
2644 Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2647 Double_t r = TMath::Sqrt(x*x + y*y);
2649 // Find angle between laser vector and R vector
2650 // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2651 Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2652 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2653 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2654 if(lasTanPhiLocIn<0)
2655 vecPhiR[bin-1]*=-1; // to have the same sign
2657 if(yprof->GetBinEntries(bin)>=10) {
2658 lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2659 TMath::Max(yprof->GetBinError(bin), 0.001));
2661 if(zprof->GetBinEntries(bin)>=10) {
2662 lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2663 TMath::Max(zprof->GetBinError(bin), 0.001));
2666 vecPhi[bin-1]=lasTanPhiLocOut;
2667 // Calculate local y from residual and database
2668 Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2671 Double_t r = TMath::Sqrt(x*x + y*y);
2674 Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2675 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2676 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2677 if(lasTanPhiLocOut<0)
2678 vecPhiR[bin-1]*=-1; // to have the same sign
2680 if(yprof->GetBinEntries(bin)>=10) {
2681 lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2682 TMath::Max(yprof->GetBinError(bin), 0.001));
2684 if(zprof->GetBinEntries(bin)>=10) {
2685 lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2686 TMath::Max(zprof->GetBinError(bin), 0.001));
2692 delete yprof; delete zprof;
2695 // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2696 nClInY = lfabsyInner.GetNpoints();
2697 if(lfabsyInner.GetNpoints()>10) {
2698 lfabsyInner.EvalRobust(0.95);
2701 lfabsyInner.GetParameters(result);
2702 yAbsInOffset = result[0];
2703 yAbsInSlope = result[1];
2704 yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2706 nClInZ = lfabszInner.GetNpoints();
2707 if(lfabszInner.GetNpoints()>10) {
2708 lfabszInner.EvalRobust(0.95);
2711 lfabszInner.GetParameters(result);
2712 zAbsInOffset = result[0];
2713 zAbsInSlope = result[1];
2714 zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2716 nClOutY = lfabsyOuter.GetNpoints();
2717 if(lfabsyOuter.GetNpoints()>10) {
2718 lfabsyOuter.EvalRobust(0.95);
2721 lfabsyOuter.GetParameters(result);
2722 yAbsOutOffset = result[0];
2723 yAbsOutSlope = result[1];
2724 yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2726 nClOutZ = lfabszOuter.GetNpoints();
2727 if(lfabszOuter.GetNpoints()>10) {
2728 lfabszOuter.EvalRobust(0.95);
2731 lfabszOuter.GetParameters(result);
2732 zAbsOutOffset = result[0];
2733 zAbsOutSlope = result[1];
2734 zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2740 Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2741 Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2742 Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2743 Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2745 Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2746 Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2748 Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2749 Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2752 (*pcstream)<<"Mean"<<
2765 "id="<<id<< // track id
2766 "entries="<<entries<< // number of entries
2767 "bz="<<bfield<< // bfield
2768 "LTr.="<<ltrp<< // refernece track
2770 "mnclIn="<<mnclIn<< // mean number of clusters in inner
2771 "mnclOut="<<mnclOut<< // mean number of clusters in outer
2772 "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2773 "mLclIn="<<mLclIn<< // mean number of clusters in inner
2774 "mLclOut="<<mLclOut<< // mean number of clusters in outer
2775 "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2776 "mdEdx="<<mdEdx<< // mean dEdx
2777 "edZfit="<<edZfit<< // entries z fit
2778 "mdZfit="<<mdZfit<< // mean z fit
2779 "rdZfit="<<rdZfit<< // RMS z fit
2782 "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2783 "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2784 "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2785 "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2786 "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2787 "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2788 "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2789 "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2790 "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2791 "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2794 "ePy1vP0="<<ePy1vP0<<
2795 "mPy1vP0="<<mPy1vP0<<
2796 "rPy1vP0="<<rPy1vP0<<
2797 "ePy2vP0="<<ePy2vP0<<
2798 "mPy2vP0="<<mPy2vP0<<
2799 "rPy2vP0="<<rPy2vP0<<
2800 "ePy3vP0="<<ePy3vP0<<
2801 "mPy3vP0="<<mPy3vP0<<
2802 "rPy3vP0="<<rPy3vP0<<
2803 "ePy1vP1="<<ePy1vP1<<
2804 "mPy1vP1="<<mPy1vP1<<
2805 "rPy1vP1="<<rPy1vP1<<
2806 "ePy2vP1="<<ePy2vP1<<
2807 "mPy2vP1="<<mPy2vP1<<
2808 "rPy2vP1="<<rPy2vP1<<
2809 "ePy3vP1="<<ePy3vP1<<
2810 "mPy3vP1="<<mPy3vP1<<
2811 "rPy3vP1="<<rPy3vP1<<
2812 "ePy2vP2In="<<ePy2vP2In<<
2813 "mPy2vP2In="<<mPy2vP2In<<
2814 "rPy2vP2In="<<rPy2vP2In<<
2815 "ePy2vP2Out="<<ePy2vP2Out<<
2816 "mPy2vP2Out="<<mPy2vP2Out<<
2817 "rPy2vP2Out="<<rPy2vP2Out<<
2818 "ePy3vP2IO="<<ePy3vP2IO<<
2819 "mPy3vP2IO="<<mPy3vP2IO<<
2820 "rPy3vP2IO="<<rPy3vP2IO<<
2823 "ePz1vP0="<<ePz1vP0<<
2824 "mPz1vP0="<<mPz1vP0<<
2825 "rPz1vP0="<<rPz1vP0<<
2826 "ePz2vP0="<<ePz2vP0<<
2827 "mPz2vP0="<<mPz2vP0<<
2828 "rPz2vP0="<<rPz2vP0<<
2829 "ePz3vP0="<<ePz3vP0<<
2830 "mPz3vP0="<<mPz3vP0<<
2831 "rPz3vP0="<<rPz3vP0<<
2832 "ePz1vP1="<<ePz1vP1<<
2833 "mPz1vP1="<<mPz1vP1<<
2834 "rPz1vP1="<<rPz1vP1<<
2835 "ePz2vP1="<<ePz2vP1<<
2836 "mPz2vP1="<<mPz2vP1<<
2837 "rPz2vP1="<<rPz2vP1<<
2838 "ePz3vP1="<<ePz3vP1<<
2839 "mPz3vP1="<<mPz3vP1<<
2840 "rPz3vP1="<<rPz3vP1<<
2841 "ePz2vP2In="<<ePz2vP2In<<
2842 "mPz2vP2In="<<mPz2vP2In<<
2843 "rPz2vP2In="<<rPz2vP2In<<
2844 "ePz2vP2Out="<<ePz2vP2Out<<
2845 "mPz2vP2Out="<<mPz2vP2Out<<
2846 "rPz2vP2Out="<<rPz2vP2Out<<
2847 "ePz3vP2IO="<<ePz3vP2IO<<
2848 "mPz3vP2IO="<<mPz3vP2IO<<
2849 "rPz3vP2IO="<<rPz3vP2IO<<
2853 "lx0="<<lxyz[0]<< // reference x
2854 "lx1="<<lxyz[1]<< // reference y
2855 "lx2="<<lxyz[2]<< // refernece z
2856 "lpx0="<<lpxyz[0]<< // reference x
2857 "lpx1="<<lpxyz[1]<< // reference y
2858 "lpx2="<<lpxyz[2]<< // refernece z
2862 "mphi="<<meanphi<< //
2863 "rmsphi="<<rmsphi<< //
2867 "mphiP="<<meanphiP<< //
2868 "rmsphiP="<<rmsphiP<< //
2877 "gp31="<<gp31<< //gaus mean - tgl
2878 "gp32="<<gp32<< //gaus rms -tgl
2879 "gp41="<<gp41<< //gaus mean - P4
2880 "gp42="<<gp42<< //gaus rms - P4
2881 // Parameters from abs res analysis
2882 "SecIn="<<secInner<< // inner sector
2883 "SecOut="<<secOuter<< // outer sector
2884 "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2885 "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
2886 "X.="<<&vecX<< // local x
2887 "Y.="<<&vecY<< // local y
2888 "R.="<<&vecR<< // radius
2889 "dY.="<<&vecDY<< // abs y residuals
2890 "dZ.="<<&vecDZ<< // abs z residuals
2891 "eY.="<<&vecEy<< // error of y residuals
2892 "eZ.="<<&vecEz<< // error of z residuals
2893 "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
2894 "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
2895 "nCl.="<<&vecN<< // number of clusters
2897 "nClInY="<<nClInY<< // Number of clusters for inner y
2898 "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2899 "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
2900 "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2901 "nClInZ="<<nClInZ<< // Number of clusters for inner z
2902 "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2903 "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
2904 "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2906 "nClOutY="<<nClOutY<< // Number of clusters for outer y
2907 "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
2908 "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
2909 "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
2910 "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
2911 "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
2912 "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
2913 "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
2920 TFile fmean("laserMean.root");
2930 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
2934 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
2935 TFile * f = pcstream->GetFile();
2937 f->mkdir("dirphiP");
2939 TF1 fp("p1","pol1");
2944 char grnamefull[1000];
2947 Double_t mphiP[100];
2948 Double_t smphi[100];
2949 Double_t smphiP[100];
2960 for (Int_t id=0; id<336; id++){
2962 sprintf(cut,"fId==%d&&%s",id,cutUser);
2963 Int_t entries = chain->Draw("bz",cut,"goff");
2964 if (entries<3) continue;
2965 AliTPCLaserTrack *ltrp = 0;
2966 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
2967 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2971 ltrp->GetPxPyPz(lpxyz);
2973 chain->Draw("bz",cut,"goff");
2974 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
2975 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
2976 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
2978 chain->Draw("gphi1",cut,"goff");
2979 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
2980 chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
2981 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
2983 chain->Draw("gphiP1",cut,"goff");
2984 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
2985 chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
2986 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
2988 chain->Draw("gz1",cut,"goff");
2989 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
2990 chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
2991 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
2994 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
2995 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3002 for (Int_t ientry=0; ientry<entries; ientry++){
3003 if (TMath::Abs(bz[ientry])<0.05){
3004 phiB0 = mphi[ientry];
3005 phiPB0 = mphiP[ientry];
3009 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3012 pphi[0] = fp.GetParameter(0); // offset
3013 pphi[1] = fp.GetParameter(1); // slope
3014 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3015 sprintf(grname,"phi_id%d",id);
3016 grphi->SetName(grname); grphi->SetTitle(grnamefull);
3017 grphi->GetXaxis()->SetTitle("b_{z} (T)");
3018 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3019 grphi->SetMaximum(1.2);
3020 grphi->SetMinimum(-1.2);
3024 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3027 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3030 pphiP[0] = fp.GetParameter(0); // offset
3031 pphiP[1] = fp.GetParameter(1); // slope
3032 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3033 sprintf(grname,"phiP_id%d",id);
3034 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3035 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3036 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3037 grphiP->SetMaximum(pphiP[0]+0.005);
3038 grphiP->SetMinimum(pphiP[0]-0.005);
3040 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3045 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3048 pmZ[0] = fp.GetParameter(0); // offset
3049 pmZ[1] = fp.GetParameter(1); // slope
3050 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3051 sprintf(grname,"mZ_id%d",id);
3052 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3053 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3054 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3056 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3062 for (Int_t ientry=0; ientry<entries; ientry++){
3063 (*pcstream)<<"Mean"<<
3066 "entries="<<entries<<
3068 "lx0="<<lxyz[0]<< // reference x
3069 "lx1="<<lxyz[1]<< // reference y
3070 "lx2="<<lxyz[2]<< // refernece z
3071 "lpx0="<<lpxyz[0]<< // reference x
3072 "lpx1="<<lpxyz[1]<< // reference y
3073 "lpx2="<<lpxyz[2]<< // refernece z
3075 "phiB0="<<phiB0<< // position shift at 0 field
3076 "phiPB0="<<phiPB0<< // angular shift at 0 field
3077 "zB0="<<zB0<< // z shift for 0 field
3079 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3080 "pphi0="<<pphi[0]<< // offset
3081 "pphi1="<<pphi[1]<< // slope
3082 "pphi2="<<pphi[2]<< // norm chi2
3084 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3085 "pphiP0="<<pphiP[0]<< // offset
3086 "pphiP1="<<pphiP[1]<< // slope
3087 "pphiP2="<<pphiP[2]<< // norm chi2
3089 "gz1="<<mZ[ientry]<<
3090 "pmZ0="<<pmZ[0]<< // offset
3091 "pmZ1="<<pmZ[1]<< // slope
3092 "pmZ2="<<pmZ[2]<< // norm chi2
3102 void AliTPCcalibLaser::Analyze(){
3109 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3111 TIterator* iter = li->MakeIterator();
3112 AliTPCcalibLaser* cal = 0;
3113 static Int_t counter0=0;
3114 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3115 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3116 // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3119 printf("Marging number %d\n", counter0);
3122 MergeFitHistos(cal);
3127 // TProfile *hp=0x0;
3128 //TProfile *hpm=0x0;
3130 for (Int_t id=0; id<336; id++){
3131 // merge fDeltaZ histograms
3132 hm = (TH1F*)cal->fDeltaZ.At(id);
3133 h = (TH1F*)fDeltaZ.At(id);
3135 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3137 fDeltaZ.AddAt(h,id);
3140 // merge fP3 histograms
3141 hm = (TH1F*)cal->fDeltaP3.At(id);
3142 h = (TH1F*)fDeltaP3.At(id);
3144 h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3146 fDeltaP3.AddAt(h,id);
3149 // merge fP4 histograms
3150 hm = (TH1F*)cal->fDeltaP4.At(id);
3151 h = (TH1F*)fDeltaP4.At(id);
3153 h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3155 fDeltaP4.AddAt(h,id);
3160 // merge fDeltaPhi histograms
3161 hm = (TH1F*)cal->fDeltaPhi.At(id);
3162 h = (TH1F*)fDeltaPhi.At(id);
3164 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3166 fDeltaPhi.AddAt(h,id);
3169 // merge fDeltaPhiP histograms
3170 hm = (TH1F*)cal->fDeltaPhiP.At(id);
3171 h = (TH1F*)fDeltaPhiP.At(id);
3173 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3175 fDeltaPhiP.AddAt(h,id);
3178 // merge fSignals histograms
3179 hm = (TH1F*)cal->fSignals.At(id);
3180 h = (TH1F*)fSignals.At(id);
3182 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3184 fSignals.AddAt(h,id);
3189 // merge ProfileY histograms -0
3190 h2m = (TH2F*)cal->fDeltaYres.At(id);
3191 h2 = (TH2F*)fDeltaYres.At(id);
3192 if (h2m&&h2) h2->Add(h2m);
3194 h2m = (TH2F*)cal->fDeltaZres.At(id);
3195 h2 = (TH2F*)fDeltaZres.At(id);
3196 if (h2m&&h2) h2->Add(h2m);
3197 // merge ProfileY histograms - 2
3198 h2m = (TH2F*)cal->fDeltaYres2.At(id);
3199 h2 = (TH2F*)fDeltaYres2.At(id);
3200 if (h2m&&h2) h2->Add(h2m);
3202 h2m = (TH2F*)cal->fDeltaZres2.At(id);
3203 h2 = (TH2F*)fDeltaZres2.At(id);
3204 if (h2m&&h2) h2->Add(h2m);
3206 // merge ProfileY histograms - abs
3207 h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3208 h2 = (TH2F*)fDeltaYresAbs.At(id);
3209 if (h2m&&h2) h2->Add(h2m);
3211 h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3212 h2 = (TH2F*)fDeltaZresAbs.At(id);
3213 if (h2m&&h2) h2->Add(h2m);
3214 // merge ProfileY histograms - 3
3215 //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3216 //h2 = (TH2F*)fDeltaYres3.At(id);
3217 //if (h2m) h2->Add(h2m);
3219 //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3220 //h2 = (TH2F*)fDeltaZres3.At(id);
3221 //if (h2m) h->Add(h2m);
3229 void AliTPCcalibLaser::MakeFitHistos(){
3231 // Make a fit histograms
3233 // Number of clusters
3235 //TH2F *fHisNclIn; //->Number of clusters inner
3236 //TH2F *fHisNclOut; //->Number of clusters outer
3237 //TH2F *fHisNclIO; //->Number of cluster inner outer
3238 // TH2F *fHisdEdx; //->dEdx histo
3239 fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3240 fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3241 fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3243 fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3244 fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3245 fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3247 fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3248 fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3253 // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3254 // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3255 // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3256 // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3257 // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3258 fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3259 fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3260 fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3261 fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3262 fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3263 // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3264 // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3265 // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3266 // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3267 // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3268 fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3269 fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3270 fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3271 fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3272 fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3277 // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3278 // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3279 // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3280 // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3281 // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3282 // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3283 // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3284 // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3285 // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3286 fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3287 fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3288 fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3289 fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3290 fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3291 fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3292 fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3293 fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3294 fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3297 // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3298 // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3299 // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3300 // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3301 // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3302 // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3303 // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3304 // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3305 // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3306 fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3307 fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3308 fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3309 fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3310 fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3311 fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3312 fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3313 fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3314 fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3316 fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3317 fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3319 fHisNclIn->SetDirectory(0); //->Number of clusters inner
3320 fHisNclOut->SetDirectory(0); //->Number of clusters outer
3321 fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3322 fHisLclIn->SetDirectory(0); //->Level arm inner
3323 fHisLclOut->SetDirectory(0); //->Level arm outer
3324 fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3325 fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3326 fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3329 fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3330 fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3331 fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3332 fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3333 fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3334 fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3335 fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3336 fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3337 fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3338 fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3341 fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3342 fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3343 fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3344 fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3345 fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3346 fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3347 fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3348 fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3349 fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3352 fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3353 fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3354 fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3355 fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3356 fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3357 fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3358 fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3359 fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3360 fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3362 fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3363 fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3370 for (Int_t id=0; id<336;id++){
3371 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3372 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3373 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3374 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
3375 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3376 TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3377 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3378 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3380 profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3381 profy->SetDirectory(0);
3382 fDeltaYres.AddAt(profy,id);
3383 profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3384 profy2->SetDirectory(0);
3385 fDeltaYres2.AddAt(profy2,id);
3386 if(!fUseFixedDriftV)
3387 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
3389 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
3390 profyabs->SetDirectory(0);
3391 fDeltaYresAbs.AddAt(profyabs,id);
3392 //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3393 //profy3->SetDirectory(0);
3394 //fDeltaYres3.AddAt(profy3,id);
3397 profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3398 profz->SetDirectory(0);
3399 fDeltaZres.AddAt(profz,id);
3400 profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3401 profz2->SetDirectory(0);
3402 fDeltaZres2.AddAt(profz2,id);
3403 if(!fUseFixedDriftV)
3404 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
3406 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
3407 profzabs->SetDirectory(0);
3408 fDeltaZresAbs.AddAt(profzabs,id);
3409 //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3410 //profz3->SetDirectory(0);
3411 //fDeltaZres3.AddAt(profz3,id);
3416 for (Int_t id=0; id<336;id++){
3417 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3418 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3419 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3421 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3422 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3423 TH1F * hisSignal = (TH1F*)fSignals.At(id);
3426 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3427 hisdz->SetDirectory(0);
3428 fDeltaZ.AddAt(hisdz,id);
3430 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3431 hisP3->SetDirectory(0);
3432 fDeltaP3.AddAt(hisP3,id);
3434 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3435 hisP4->SetDirectory(0);
3436 fDeltaP4.AddAt(hisP4,id);
3439 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3440 hisdphi->SetDirectory(0);
3441 fDeltaPhi.AddAt(hisdphi,id);
3443 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3444 hisdphiP->SetDirectory(0);
3445 fDeltaPhiP.AddAt(hisdphiP,id);
3446 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3447 hisSignal->SetDirectory(0);
3448 fSignals.AddAt(hisSignal,id);
3452 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3453 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3454 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3455 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3459 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3460 Int_t binsLaser[12]= {336, //id
3472 Double_t xminLaser[12]= {0, //id
3484 Double_t xmaxLaser[12]= {336, //id
3497 TString nameLaser[12]= {"id",
3509 TString titleLaser[12]= {"id",
3521 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3522 for (Int_t iaxis=1; iaxis<12; iaxis++){
3523 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3524 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3528 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3530 // Merge content of histograms
3532 // Only first histogram is checked - all other should be the same
3533 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3535 if (!laser->fHisNclIn) return; // empty histograms
3536 if (!fHisNclIn) MakeFitHistos();
3538 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3539 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3540 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3541 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3542 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3543 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3544 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3545 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3548 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3549 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3550 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3551 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3552 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3553 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3554 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3555 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3556 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3557 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3560 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3561 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3562 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3563 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3564 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3565 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3566 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3567 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3568 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3571 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3572 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3573 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3574 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3575 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3576 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3577 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3578 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3579 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3580 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3581 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3594 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3596 // Dump fit information - collect information from the streamers
3599 TChain * chainFit=0;
3600 TChain * chainTrack=0;
3603 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3604 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3605 AliXRDPROOFtoolkit tool;
3606 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3607 chainTrack->Lookup();
3608 chainTrack->SetProof(kTRUE);
3609 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3610 chainDrift->Lookup();
3611 chainDrift->SetProof(kTRUE);
3613 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3615 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3617 chainFit->SetProof(kTRUE);
3618 chain->SetProof(kTRUE);
3619 AliTPCLaserTrack::LoadTracks();
3620 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3626 TCut cutP3("abs(Tr.fP[3])<0.1");
3627 TCut cutP4("abs(Tr.fP[4])<0.5");
3628 TCut cutPx = cutP3+cutP4;
3629 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3630 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3631 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3632 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3634 TCut cutdEdx("sqrt(dEdx)>3");
3635 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3636 TCut cutN("nclO>20&&nclI>20");
3637 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3641 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3642 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3643 TCut cutClX("abs(Cl[].fX)>10");
3644 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3645 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3646 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3647 TCut cutQ("sqrt(Cl[].fMax)>4");
3648 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3652 TH1F * phisAccept = 0;
3654 TProfile * pdEdx = 0;
3662 TProfile * pNclI = 0;
3663 TProfile * pNclO = 0;
3665 TProfile * pchi2YIn =0;
3666 TProfile * pchi2ZIn =0;
3667 TProfile * pchi2YOut =0;
3668 TProfile * pchi2ZOut =0;
3669 TProfile * pchi2YInOut =0;
3670 TProfile * pchi2ZInOut =0;;
3672 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3673 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3674 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3675 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3676 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3677 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3679 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3680 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3683 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3684 pP0 = (TProfile*)gROOT->FindObject("hP0");
3685 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3686 pP1 = (TProfile*)gROOT->FindObject("hP1");
3687 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3688 pP2 = (TProfile*)gROOT->FindObject("hP2");
3689 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3690 pP3 = (TProfile*)gROOT->FindObject("hP3");
3691 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3692 pP4 = (TProfile*)gROOT->FindObject("hP4");
3694 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3695 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3696 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3697 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3700 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3701 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3702 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3703 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3704 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3705 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3706 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3707 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3708 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3709 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3710 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3711 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3713 // second derivatives
3715 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3716 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3717 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3718 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3719 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3720 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3722 phisPy2In->FitSlicesY();
3723 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3724 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3725 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3727 phisPy2Out->FitSlicesY();
3728 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3729 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3730 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3732 phisPy2InOut->FitSlicesY();
3733 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3734 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3735 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3737 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3738 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3739 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3740 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3741 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3742 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3744 phisPz2In->FitSlicesY();
3745 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3746 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3747 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3749 phisPz2Out->FitSlicesY();
3750 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3751 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3752 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3754 phisPz2InOut->FitSlicesY();
3755 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3756 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3757 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3764 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3765 for (Int_t ilaser=0; ilaser<336; ilaser++){
3766 Float_t all=phisAl->GetBinContent(ilaser+1);
3767 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3768 Float_t out=phisOut->GetBinContent(ilaser+1);
3769 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3770 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3771 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3772 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3773 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3774 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3777 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3778 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3780 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3781 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3782 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3783 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3784 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3785 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3787 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3788 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3789 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3791 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3792 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3793 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3795 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3796 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3797 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3799 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3800 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3801 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3803 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3804 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3805 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3807 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3808 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3809 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3811 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3812 (*pcstream)<<"Scan"<<
3826 "chi2YIn="<<chi2YIn<<
3827 "chi2YOut="<<chi2YOut<<
3828 "chi2YInOut="<<chi2YInOut<<
3829 "chi2ZIn="<<chi2ZIn<<
3830 "chi2ZOut="<<chi2ZOut<<
3831 "chi2ZInOut="<<chi2ZInOut<<
3833 "nPy2In="<<entriesPy2In<<
3834 "mPy2In="<<meanPy2In<<
3835 "sPy2In="<<sigmaPy2In<<
3837 "nPy2Out="<<entriesPy2Out<<
3838 "mPy2Out="<<meanPy2Out<<
3839 "sPy2Out="<<sigmaPy2Out<<
3841 "nPy2InOut="<<entriesPy2InOut<<
3842 "mPy2InOut="<<meanPy2InOut<<
3843 "sPy2InOut="<<sigmaPy2InOut<<
3845 "nPz2In="<<entriesPz2In<<
3846 "mPz2In="<<meanPz2In<<
3847 "sPz2In="<<sigmaPz2In<<
3849 "nPz2Out="<<entriesPz2Out<<
3850 "mPz2Out="<<meanPz2Out<<
3851 "sPz2Out="<<sigmaPz2Out<<
3853 "nPz2InOut="<<entriesPz2InOut<<
3854 "mPz2InOut="<<meanPz2InOut<<
3855 "sPz2InOut="<<sigmaPz2InOut<<
3863 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3864 TVectorD& meanSlope,
3865 TVectorD& sectorArray,
3868 // This method should ideally go in AliTPCLaser
3869 // option == 0 (pads outer - closest to beam)
3870 // option == 1 (pads inner)
3871 // option == 2 (time outer)
3872 // option == 3 (time inner)
3873 Int_t nFailures = 0;
3875 for(Int_t id = 0; id < 336; id++) {
3877 if (!AliTPCLaserTrack::GetTracks())
3878 AliTPCLaserTrack::LoadTracks();
3879 AliTPCLaserTrack *ltrp =
3880 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3882 AliExternalTrackParam trackParam(*ltrp);
3884 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3885 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3886 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3888 Double_t angle = trackParam.GetAlpha();
3890 angle += 2*TMath::Pi();
3891 if(trackParam.GetSnp()>0) // track points to sector "before"
3892 angle -= deltaangle*TMath::DegToRad();
3893 else // track points to sector "after"
3894 angle += deltaangle*TMath::DegToRad();
3896 Bool_t success = trackParam.Rotate(angle);
3899 // cout << "WARNING: Rotate failed for ID: " << id << endl;
3903 angle *= TMath::RadToDeg();
3905 Int_t sector = TMath::Nint((angle-10.0)/20.0);
3910 if(ltrp->GetSide()==1) // C side
3912 if(option==0 || option==2)
3914 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
3917 sectorArray[id] = sector;
3919 const Double_t x0 = 0;
3921 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3922 Double_t slopez = trackParam.GetTgl();
3923 // One needs a factor sqrt(1+slopey**2) to take into account the
3924 // longer path length
3925 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3926 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3928 // Double_t offsetz = trackParam.GetZ();
3929 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
3930 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
3931 if(option==2 || option==3) {
3932 meanOffset[id] = offsetz; meanSlope[id] = slopez;
3934 meanOffset[id] = offsety; meanSlope[id] = slopey;
3939 AliWarning(Form("Rotate method failed %d times", nFailures));
3944 TFile f("vscan.root");
3949 chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
3951 chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
3954 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);
3957 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)
3966 // check edge effects
3967 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
3969 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
3971 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)
3975 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
3976 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")
3987 dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
3990 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)
3992 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)
3998 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)
4000 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)
4004 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)
4006 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)
4015 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutDY,"prof")
4017 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&<r.fP[1]<0"+cutA+cutDY,"prof")
4021 chainFit->Draw("LTr.fId","nclI>10",100000)
4023 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4025 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4027 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4037 gSystem->Load("libSTAT.so")
4038 TStatToolkit toolkit;
4044 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
4049 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
4050 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
4051 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
4052 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
4054 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
4055 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
4056 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
4057 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
4059 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
4060 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
4061 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
4062 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
4067 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4069 treeT->SetAlias("fit",strq0->Data());
4072 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4074 treeT->SetAlias("fitP",strqP->Data());
4077 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4078 treeT->SetAlias("fitD",strqDrift->Data());
4081 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4083 for (Int_t i=0; i<6;i++){
4084 treeT->SetLineColor(i+2);
4085 treeT->SetMarkerSize(1);
4086 treeT->SetMarkerStyle(22+i);
4087 treeT->SetMarkerColor(i+2);
4089 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4097 TTree * tree = (TTree*)f.Get("FitModels");
4099 TEventList listLFit0("listLFit0","listLFit0");
4100 TEventList listLFit1("listLFit1","listLFit1");
4101 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4102 tree->SetEventList(&listLFit0);
4107 gSystem->Load("libSTAT.so")
4108 TStatToolkit toolkit;
4114 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4115 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4119 fstring+="cos(dp)++";
4120 fstring+="sin(dp)++";
4121 fstring+="cos(dt)++";
4122 fstring+="sin(dt)++";
4124 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4137 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4138 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4139 AliXRDPROOFtoolkit tool;
4140 TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4141 chainTrack->Lookup();
4142 chainTrack->SetProof(kTRUE);
4144 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4146 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4148 chainFit->SetProof(kTRUE);
4149 chain->SetProof(kTRUE);
4153 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4154 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4155 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4156 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4158 TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4159 TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4160 TCut cutN("nclO>20&&nclI>20");
4161 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4165 TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4166 TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4167 TCut cutClX("abs(Cl.fX)>10");
4168 TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4169 TCut cutCl=cutClY+cutClZ+cutClX;
4172 // check edge effects
4173 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4175 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4177 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)
4181 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4182 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&<r.fP[1]<0"+cutA+cutD,"prof")