Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
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 : **************************************************************************/
15 :
16 : /*
17 : //
18 : // FUNCTIONALITY:
19 : //
20 : // 1. The laser track is associated with the mirror
21 : // see function FindMirror
22 : //
23 : // 2. The laser track is accepted for the analysis under certain condition
24 : // (see function Accpet laser)
25 : //
26 : // 3. The drift velocity and jitter is calculated event by event
27 : // (see function drift velocity)
28 : //
29 : // 4. The tracks are refitted at different sectors
30 : // Fit model
31 : // 4.a) line
32 : // 4.b) parabola
33 : // 4.c) parabola with common P2 for inner and outer
34 : //
35 : // To make laser scan the user interaction neccessary
36 : //
37 : .x ~/NimStyle.C
38 : gSystem->Load("libANALYSIS");
39 : gSystem->Load("libTPCcalib");
40 : TFile fcalib("CalibObjectsTrain2.root");
41 : AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
42 : laser->DumpMeanInfo(run)
43 : TFile fmean("laserMean.root")
44 : //
45 : // laser track clasification;
46 : //
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);
52 :
53 : //
54 : //
55 : // Analyze LASER scan
56 : //
57 :
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);
62 : chain->Lookup();
63 : AliTPCcalibLaser::DumpScanInfo(chain)
64 : TFile fscan("laserScan.root")
65 : TTree * treeT = (TTree*)fscan.Get("Mean")
66 : //
67 : // Analyze laser
68 : //
69 : gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
70 : gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
71 : AliXRDPROOFtoolkit tool;
72 : AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1)
73 : TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
74 : chainDrift->Lookup();
75 : TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
76 : chainDriftN->Lookup();
77 :
78 :
79 : TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
80 : chain->Lookup();
81 : TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
82 : chainFit->Lookup();
83 : TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
84 : chainTrack->Lookup();
85 :
86 : */
87 :
88 :
89 :
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"
97 : #include "TH1D.h"
98 : #include "TH1F.h"
99 : #include "TProfile.h"
100 : #include "TVectorD.h"
101 : #include "TTreeStream.h"
102 : #include "TFile.h"
103 : #include "TF1.h"
104 : #include "TGraphErrors.h"
105 : #include "AliTPCclusterMI.h"
106 : #include "AliTPCseed.h"
107 : #include "AliTracker.h"
108 : #include "AliLog.h"
109 : #include "TClonesArray.h"
110 : #include "TPad.h"
111 : #include "TSystem.h"
112 : #include "TCut.h"
113 : #include "TChain.h"
114 : #include "TH2F.h"
115 : #include "TStatToolkit.h"
116 : #include "TROOT.h"
117 : #include "TDatabasePDG.h"
118 :
119 :
120 : #include "TTreeStream.h"
121 : #include <iostream>
122 : #include <sstream>
123 : #include "AliTPCLaserTrack.h"
124 : #include "AliTPCcalibDB.h"
125 : #include "AliTPCParam.h"
126 : #include "AliTPCParamSR.h"
127 : #include "TTimeStamp.h"
128 : #include "AliDCSSensorArray.h"
129 : #include "AliDCSSensor.h"
130 : #include "AliGRPObject.h"
131 : #include "AliTPCROC.h"
132 : #include "AliTPCreco.h"
133 :
134 : using namespace std;
135 :
136 6 : ClassImp(AliTPCcalibLaser)
137 :
138 : AliTPCcalibLaser::AliTPCcalibLaser():
139 0 : AliTPCcalibBase(),
140 0 : fESD(0),
141 0 : fESDfriend(0),
142 0 : fNtracks(0),
143 0 : fTracksMirror(336),
144 0 : fTracksEsd(336),
145 0 : fTracksEsdParam(336),
146 0 : fTracksTPC(336),
147 0 : fFullCalib(kTRUE),
148 0 : fDeltaZ(336),
149 0 : fDeltaP3(336),
150 0 : fDeltaP4(336),
151 0 : fDeltaPhi(336),
152 0 : fDeltaPhiP(336),
153 0 : fSignals(336),
154 : //
155 0 : fHisLaser(0), // N dim histogram of laser
156 0 : fHisLaserPad(0), // N dim histogram of laser
157 0 : fHisLaserTime(0), // N dim histogram of laser
158 0 : fHisNclIn(0), //->Number of clusters inner
159 0 : fHisNclOut(0), //->Number of clusters outer
160 0 : fHisNclIO(0), //->Number of cluster inner outer
161 0 : fHisLclIn(0), //->Level arm inner
162 0 : fHisLclOut(0), //->Level arm outer
163 0 : fHisLclIO(0), //->Number of cluster inner outer
164 0 : fHisdEdx(0), //->dEdx histo
165 0 : fHisdZfit(0), //->distance to the mirror after linear fit
166 : //
167 : //
168 0 : fHisChi2YIn1(0), //->chi2 y inner - line
169 0 : fHisChi2YOut1(0), //->chi2 y inner - line
170 0 : fHisChi2YIn2(0), //->chi2 y inner - parabola
171 0 : fHisChi2YOut2(0), //->chi2 y inner - parabola
172 0 : fHisChi2YIO1(0), //->chi2 y IO - common
173 0 : fHisChi2ZIn1(0), //->chi2 z inner - line
174 0 : fHisChi2ZOut1(0), //->chi2 z inner - line
175 0 : fHisChi2ZIn2(0), //->chi2 z inner - parabola
176 0 : fHisChi2ZOut2(0), //->chi2 z inner - parabola
177 0 : fHisChi2ZIO1(0), //->chi2 z IO - common
178 : //
179 : //
180 0 : fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
181 0 : fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
182 0 : fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
183 0 : fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
184 0 : fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
185 0 : fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
186 0 : fHisPy2vP2In(0), //-> Curv P2inner - parabola
187 0 : fHisPy2vP2Out(0), //-> Curv P2outer - parabola
188 0 : fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
189 : //
190 : //
191 0 : fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
192 0 : fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
193 0 : fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
194 0 : fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
195 0 : fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
196 0 : fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
197 0 : fHisPz2vP2In(0), //-> Curv P2inner - parabola
198 0 : fHisPz2vP2Out(0), //-> Curv P2outer - parabola
199 0 : fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
200 : //
201 0 : fDeltaYres(336), //->2D histo of residuals
202 0 : fDeltaZres(336), //->2D histo fo residuals
203 0 : fDeltaYres2(336), //->2D histo of residuals
204 0 : fDeltaZres2(336), //->2D histo fo residuals
205 0 : fDeltaYresAbs(336), //->2D histo of residuals
206 0 : fHisYAbsErrors(0), //-> total errors per beam in the abs res y analysis
207 0 : fDeltaZresAbs(336), //->2D histo of residuals
208 0 : fHisZAbsErrors(0), //-> total errors per beam in the abs res z analysis
209 : //fDeltaYres3(336), //->2D histo of residuals
210 : //fDeltaZres3(336), //->2D histo fo residuals
211 0 : fFitAside(new TVectorD(5)),
212 0 : fFitCside(new TVectorD(5)),
213 0 : fFitACside(new TVectorD(6)),
214 0 : fEdgeXcuts(3),
215 0 : fEdgeYcuts(3),
216 0 : fNClCuts(5),
217 0 : fNcuts(0),
218 0 : fBeamSectorOuter(336),
219 0 : fBeamSectorInner(336),
220 0 : fBeamOffsetYOuter(336),
221 0 : fBeamSlopeYOuter(336),
222 0 : fBeamOffsetYInner(336),
223 0 : fBeamSlopeYInner(336),
224 0 : fBeamOffsetZOuter(336),
225 0 : fBeamSlopeZOuter(336),
226 0 : fBeamOffsetZInner(336),
227 0 : fBeamSlopeZInner(336),
228 0 : fInverseSlopeZ(kTRUE),
229 0 : fUseFixedDriftV(0),
230 0 : fFixedFitAside0(0.0),
231 0 : fFixedFitAside1(1.0),
232 0 : fFixedFitCside0(0.0),
233 0 : fFixedFitCside1(1.0)
234 0 : {
235 : //
236 : // Constructor
237 : //
238 0 : fTracksEsdParam.SetOwner(kTRUE);
239 0 : for (Int_t i=0; i<336; i++) {
240 0 : fFitZ[i]=0;
241 0 : fCounter[i]=0; //! counter of usage
242 0 : fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
243 0 : fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
244 : }
245 0 : }
246 :
247 : AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
248 0 : AliTPCcalibBase(),
249 0 : fESD(0),
250 0 : fESDfriend(0),
251 0 : fNtracks(0),
252 0 : fTracksMirror(336),
253 0 : fTracksEsd(336),
254 0 : fTracksEsdParam(336),
255 0 : fTracksTPC(336),
256 0 : fFullCalib(full),
257 : //
258 0 : fDeltaZ(336), // array of histograms of delta z for each track
259 0 : fDeltaP3(336), // array of histograms of delta z for each track
260 0 : fDeltaP4(336), // array of histograms of P3 for each track
261 0 : fDeltaPhi(336), // array of histograms of P4 for each track
262 0 : fDeltaPhiP(336), // array of histograms of delta z for each track
263 0 : fSignals(336), // array of dedx signals
264 : //
265 : //
266 0 : fHisLaser(0), // N dim histogram of laser
267 0 : fHisLaserPad(0), // N dim histogram of laser
268 0 : fHisLaserTime(0), // N dim histogram of laser
269 :
270 0 : fHisNclIn(0), //->Number of clusters inner
271 0 : fHisNclOut(0), //->Number of clusters outer
272 0 : fHisNclIO(0), //->Number of cluster inner outer
273 0 : fHisLclIn(0), //->Level arm inner
274 0 : fHisLclOut(0), //->Level arm outer
275 0 : fHisLclIO(0), //->Number of cluster inner outer
276 0 : fHisdEdx(0), //->dEdx histo
277 0 : fHisdZfit(0), //->distance to the mirror after linear fit
278 : //
279 : //
280 0 : fHisChi2YIn1(0), //->chi2 y inner - line
281 0 : fHisChi2YOut1(0), //->chi2 y inner - line
282 0 : fHisChi2YIn2(0), //->chi2 y inner - parabola
283 0 : fHisChi2YOut2(0), //->chi2 y inner - parabola
284 0 : fHisChi2YIO1(0), //->chi2 y IO - common
285 0 : fHisChi2ZIn1(0), //->chi2 z inner - line
286 0 : fHisChi2ZOut1(0), //->chi2 z inner - line
287 0 : fHisChi2ZIn2(0), //->chi2 z inner - parabola
288 0 : fHisChi2ZOut2(0), //->chi2 z inner - parabola
289 0 : fHisChi2ZIO1(0), //->chi2 z IO - common
290 : //
291 : //
292 0 : fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
293 0 : fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
294 0 : fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
295 0 : fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
296 0 : fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
297 0 : fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
298 0 : fHisPy2vP2In(0), //-> Curv P2inner - parabola
299 0 : fHisPy2vP2Out(0), //-> Curv P2outer - parabola
300 0 : fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
301 : //
302 : //
303 0 : fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
304 0 : fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
305 0 : fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
306 0 : fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
307 0 : fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
308 0 : fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
309 0 : fHisPz2vP2In(0), //-> Curv P2inner - parabola
310 0 : fHisPz2vP2Out(0), //-> Curv P2outer - parabola
311 0 : fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
312 : //
313 : //
314 : //
315 0 : fDeltaYres(336),
316 0 : fDeltaZres(336),
317 0 : fDeltaYres2(336),
318 0 : fDeltaZres2(336),
319 0 : fDeltaYresAbs(336),
320 0 : fHisYAbsErrors(0),
321 0 : fDeltaZresAbs(336),
322 0 : fHisZAbsErrors(0),
323 : // fDeltaYres3(336),
324 : //fDeltaZres3(336),
325 0 : fFitAside(new TVectorD(5)), // drift fit - A side
326 0 : fFitCside(new TVectorD(5)), // drift fit - C- side
327 0 : fFitACside(new TVectorD(6)), // drift fit - AC- side
328 0 : fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
329 0 : fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
330 0 : fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
331 0 : fNcuts(0), // number of cuts
332 0 : fBeamSectorOuter(336),
333 0 : fBeamSectorInner(336),
334 0 : fBeamOffsetYOuter(336),
335 0 : fBeamSlopeYOuter(336),
336 0 : fBeamOffsetYInner(336),
337 0 : fBeamSlopeYInner(336),
338 0 : fBeamOffsetZOuter(336),
339 0 : fBeamSlopeZOuter(336),
340 0 : fBeamOffsetZInner(336),
341 0 : fBeamSlopeZInner(336),
342 0 : fInverseSlopeZ(kTRUE),
343 0 : fUseFixedDriftV(0),
344 0 : fFixedFitAside0(0.0),
345 0 : fFixedFitAside1(1.0),
346 0 : fFixedFitCside0(0.0),
347 0 : fFixedFitCside1(1.0)
348 0 : {
349 0 : SetName(name);
350 0 : SetTitle(title);
351 : //
352 : // Constructor
353 : //
354 0 : fTracksEsdParam.SetOwner(kTRUE);
355 0 : for (Int_t i=0; i<336; i++) {
356 0 : fFitZ[i]=0;
357 0 : fCounter[i]=0; //! counter of usage
358 0 : fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
359 0 : fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
360 : }
361 0 : }
362 :
363 : AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
364 0 : AliTPCcalibBase(calibLaser),
365 0 : fESD(0),
366 0 : fESDfriend(0),
367 0 : fNtracks(0),
368 0 : fTracksMirror(336),
369 0 : fTracksEsd(336),
370 0 : fTracksEsdParam(336),
371 0 : fTracksTPC(336),
372 0 : fFullCalib(calibLaser.fFullCalib),
373 : //
374 0 : fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
375 0 : fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
376 0 : fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
377 0 : fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
378 0 : fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
379 0 : fSignals(((calibLaser.fSignals))), // array of dedx signals
380 : //
381 : //
382 0 : fHisLaser(0), // N dim histogram of laser
383 0 : fHisLaserPad(0), // N dim histogram of laser
384 0 : fHisLaserTime(0), // N dim histogram of laser
385 :
386 0 : fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
387 0 : fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
388 0 : fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
389 0 : fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
390 0 : fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
391 0 : fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
392 0 : fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
393 0 : fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
394 : //
395 : //
396 0 : fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
397 0 : fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
398 0 : fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
399 0 : fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
400 0 : fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
401 0 : fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
402 0 : fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
403 0 : fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
404 0 : fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
405 0 : fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
406 : //
407 : //
408 0 : fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
409 0 : fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
410 0 : fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
411 0 : fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
412 0 : fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
413 0 : fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
414 0 : fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
415 0 : fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
416 0 : fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
417 : //
418 : //
419 0 : fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
420 0 : fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
421 0 : fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
422 0 : fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
423 0 : fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
424 0 : fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
425 0 : fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
426 0 : fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
427 0 : fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
428 : //
429 : //
430 0 : fDeltaYres(((calibLaser.fDeltaYres))),
431 0 : fDeltaZres(((calibLaser.fDeltaZres))),
432 0 : fDeltaYres2(((calibLaser.fDeltaYres))),
433 0 : fDeltaZres2(((calibLaser.fDeltaZres))),
434 0 : fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
435 0 : fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
436 0 : fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
437 0 : fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
438 : // fDeltaYres3(((calibLaser.fDeltaYres))),
439 : //fDeltaZres3(((calibLaser.fDeltaZres))),
440 0 : fFitAside(new TVectorD(5)), // drift fit - A side
441 0 : fFitCside(new TVectorD(5)), // drift fit - C- side
442 0 : fFitACside(new TVectorD(6)), // drift fit - C- side
443 0 : fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
444 0 : fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
445 0 : fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
446 0 : fNcuts(0), // number of cuts
447 0 : fBeamSectorOuter(336),
448 0 : fBeamSectorInner(336),
449 0 : fBeamOffsetYOuter(336),
450 0 : fBeamSlopeYOuter(336),
451 0 : fBeamOffsetYInner(336),
452 0 : fBeamSlopeYInner(336),
453 0 : fBeamOffsetZOuter(336),
454 0 : fBeamSlopeZOuter(336),
455 0 : fBeamOffsetZInner(336),
456 0 : fBeamSlopeZInner(336),
457 0 : fInverseSlopeZ(calibLaser.fInverseSlopeZ),
458 0 : fUseFixedDriftV(calibLaser.fUseFixedDriftV),
459 0 : fFixedFitAside0(calibLaser.fFixedFitAside0),
460 0 : fFixedFitAside1(calibLaser.fFixedFitAside1),
461 0 : fFixedFitCside0(calibLaser.fFixedFitCside0),
462 0 : fFixedFitCside1(calibLaser.fFixedFitCside1)
463 0 : {
464 : //
465 : // copy constructor
466 : //
467 0 : for (Int_t i=0; i<336; i++) {
468 0 : fFitZ[i]=0;
469 0 : fCounter[i]=0; //! counter of usage
470 0 : fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
471 0 : fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
472 : }
473 0 : }
474 :
475 :
476 :
477 : AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
478 : //
479 : // assgnment operator
480 : //
481 0 : if (this != &calibLaser) {
482 0 : new (this) AliTPCcalibLaser(calibLaser);
483 : }
484 0 : return *this;
485 :
486 0 : }
487 :
488 :
489 :
490 :
491 0 : AliTPCcalibLaser::~AliTPCcalibLaser() {
492 : //
493 : // destructor
494 : //
495 0 : if ( fHisNclIn){
496 0 : delete fHisLaser; //->
497 0 : delete fHisLaserPad; //->
498 0 : delete fHisLaserTime; //->
499 :
500 0 : delete fHisNclIn; //->Number of clusters inner
501 0 : delete fHisNclOut; //->Number of clusters outer
502 0 : delete fHisNclIO; //->Number of cluster inner outer
503 0 : delete fHisLclIn; //->Level arm inner
504 0 : delete fHisLclOut; //->Level arm outer
505 0 : delete fHisLclIO; //->Number of cluster inner outer
506 0 : delete fHisdEdx;
507 0 : delete fHisdZfit;
508 : //
509 : //
510 0 : delete fHisChi2YIn1; //->chi2 y inner - line
511 0 : delete fHisChi2YOut1; //->chi2 y inner - line
512 0 : delete fHisChi2YIn2; //->chi2 y inner - parabola
513 0 : delete fHisChi2YOut2; //->chi2 y inner - parabola
514 0 : delete fHisChi2YIO1; //->chi2 y IO - common
515 0 : delete fHisChi2ZIn1; //->chi2 z inner - line
516 0 : delete fHisChi2ZOut1; //->chi2 z inner - line
517 0 : delete fHisChi2ZIn2; //->chi2 z inner - parabola
518 0 : delete fHisChi2ZOut2; //->chi2 z inner - parabola
519 0 : delete fHisChi2ZIO1; //->chi2 z IO - common
520 : //
521 : //
522 0 : delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
523 0 : delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
524 0 : delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
525 0 : delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
526 0 : delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
527 0 : delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
528 0 : delete fHisPy2vP2In; //-> Curv P2inner - parabola
529 0 : delete fHisPy2vP2Out; //-> Curv P2outer - parabola
530 0 : delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
531 : //
532 : //
533 0 : delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
534 0 : delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
535 0 : delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
536 0 : delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
537 0 : delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
538 0 : delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
539 0 : delete fHisPz2vP2In; //-> Curv P2inner - parabola
540 0 : delete fHisPz2vP2Out; //-> Curv P2outer - parabola
541 0 : delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
542 : }
543 : //
544 : //
545 : //
546 0 : fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
547 0 : fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
548 0 : fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
549 0 : fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
550 0 : fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
551 0 : fSignals.SetOwner(); //->Array of dedx signals
552 :
553 0 : fDeltaZ.Delete(); //-> array of histograms of delta z for each track
554 0 : fDeltaP3.Delete(); //-> array of histograms of P3 for each track
555 0 : fDeltaP4.Delete(); //-> array of histograms of P4 for each track
556 0 : fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
557 0 : fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
558 0 : fSignals.Delete(); //->Array of dedx signals
559 :
560 0 : fDeltaYres.SetOwner();
561 0 : fDeltaYres.Delete();
562 0 : delete fHisYAbsErrors;
563 0 : fDeltaZres.SetOwner();
564 0 : fDeltaZres.Delete();
565 0 : delete fHisZAbsErrors;
566 0 : fDeltaYres2.SetOwner();
567 0 : fDeltaYres2.Delete();
568 0 : fDeltaZres2.SetOwner();
569 0 : fDeltaZres2.Delete();
570 :
571 0 : fDeltaYresAbs.SetOwner();
572 0 : fDeltaYresAbs.Delete();
573 0 : fDeltaZresAbs.SetOwner();
574 0 : fDeltaZresAbs.Delete();
575 0 : }
576 :
577 :
578 :
579 : void AliTPCcalibLaser::Process(AliESDEvent * event) {
580 : //
581 : //
582 : // Loop over tracks and call Process function
583 : //
584 : const Int_t kMinTracks=20;
585 : const Int_t kMinClusters=40;
586 :
587 0 : fESD = event;
588 0 : if (!fESD) {
589 0 : return;
590 : }
591 0 : fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
592 0 : if (!fESDfriend) {
593 0 : return;
594 : }
595 0 : if (fESDfriend->TestSkipBit()) return;
596 0 : if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
597 0 : AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
598 : //
599 : // find CE background if present
600 : //
601 0 : if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
602 0 : TH1D hisCE("hhisCE","hhisCE",100,-100,100);
603 0 : for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
604 0 : AliESDtrack *track=fESD->GetTrack(i);
605 0 : if (!track) continue;
606 0 : hisCE.Fill(track->GetZ());
607 0 : hisCE.Fill(track->GetZ()+2);
608 0 : hisCE.Fill(track->GetZ()-2);
609 0 : }
610 : //
611 : //
612 :
613 :
614 0 : fTracksTPC.Clear();
615 0 : fTracksEsd.Clear();
616 0 : fTracksEsdParam.Delete();
617 0 : for (Int_t id=0; id<336;id++) {
618 0 : fCounter[id]=0;
619 0 : fClusterCounter[id]=0;
620 0 : fClusterSatur[id]=0;
621 : }
622 : //
623 0 : Int_t n=fESD->GetNumberOfTracks();
624 : Int_t counter=0;
625 0 : for (Int_t i=0;i<n;++i) {
626 0 : AliESDtrack *track=fESD->GetTrack(i);
627 0 : if (!track) continue;
628 0 : AliESDfriendTrack *friendTrack=(AliESDfriendTrack*)track->GetFriendTrack();
629 0 : if (!friendTrack) continue;
630 0 : Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
631 0 : if (binC>336) continue; //remove CE background
632 : TObject *calibObject=0;
633 : AliTPCseed *seed=0;
634 0 : for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
635 0 : if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
636 0 : break;
637 0 : if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
638 : //filter CE tracks
639 0 : Int_t id = FindMirror(track,seed);
640 0 : if (id>=0) counter++;
641 0 : }
642 : //
643 0 : }
644 0 : fNtracks=counter;
645 0 : if (counter<kMinTracks) return;
646 :
647 : //FitDriftV();
648 0 : FitDriftV(0.2);
649 0 : if (!fFullCalib) return;
650 : static Bool_t init=kFALSE;
651 0 : if (!init){
652 0 : init = kTRUE; // way around for PROOF - to be investigated
653 0 : UpdateFitHistos();
654 : }
655 : //
656 0 : for (Int_t id=0; id<336; id++){
657 : //
658 0 : if (!fTracksEsdParam.At(id)) continue;
659 0 : if (fClusterSatur[id]>0.3) continue; // tracks in saturation
660 0 : DumpLaser(id);
661 0 : if ( AcceptLaser(id) && fFitZ[id]<0.5){
662 0 : RefitLaserJW(id);
663 0 : MakeDistHisto(id);
664 : }
665 : }
666 :
667 0 : }
668 :
669 : void AliTPCcalibLaser::MakeDistHisto(Int_t id){
670 : //
671 : //
672 : //
673 : // for (Int_t id=0; id<336; id++){
674 : //
675 : //
676 0 : if (!fTracksEsdParam.At(id)) return;
677 0 : if (!AcceptLaser(id)) return;
678 0 : Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
679 : //
680 : //
681 0 : TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
682 0 : if (!hisdz) UpdateFitHistos();
683 0 : hisdz = (TH1F*)fDeltaZ.At(id);
684 0 : TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
685 0 : TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
686 0 : TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
687 0 : TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
688 0 : TH1F * hisSignal = (TH1F*)fSignals.At(id);
689 : //
690 :
691 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
692 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
693 0 : AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
694 0 : if (!param) return;
695 0 : if (!ltrp) return;
696 0 : if (!track) return;
697 0 : Double_t xyz[3];
698 0 : Double_t pxyz[3];
699 0 : Double_t lxyz[3];
700 0 : Double_t lpxyz[3];
701 0 : param->GetXYZ(xyz);
702 0 : param->GetPxPyPz(pxyz);
703 0 : ltrp->GetXYZ(lxyz);
704 0 : ltrp->GetPxPyPz(lpxyz);
705 : //
706 0 : Float_t dz = param->GetZ()-ltrp->GetZ();
707 0 : Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
708 0 : Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
709 0 : if (hisdz) hisdz->Fill(dz);
710 0 : if (hisP3) hisP3->Fill(param->GetParameter()[3]);
711 0 : if (hisP4) hisP4->Fill(param->GetParameter()[4]);
712 0 : if (hisdphi) hisdphi->Fill(dphi);
713 0 : if (hisdphiP) hisdphiP->Fill(dphiP);
714 0 : if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
715 : // fill HisLaser
716 0 : xhis[0] = ltrp->GetId();
717 0 : xhis[1] = ltrp->GetSide();
718 0 : xhis[2] = ltrp->GetRod();
719 0 : xhis[3] = ltrp->GetBundle();
720 0 : xhis[4] = ltrp->GetBeam();
721 0 : xhis[5] = dphi;
722 0 : xhis[6] = fFitZ[id];
723 0 : xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
724 0 : xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
725 0 : xhis[9] = param->GetParameter()[4];
726 0 : xhis[10]= track->GetTPCNcls();
727 0 : xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
728 : // }
729 0 : fHisLaser->Fill(xhis);
730 : //
731 :
732 0 : }
733 :
734 : void AliTPCcalibLaser::FitDriftV(){
735 : //
736 : // Fit corrections to the drift velocity - linear approximation in the z and global y
737 : //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
738 : //
739 : /*
740 : Formulas:
741 :
742 : z = s* (z0 - vd*(t-t0))
743 :
744 : s - side -1 and +1
745 : t0 - time 0
746 : vd - nominal drift velocity
747 : zs - miscalibrated position
748 :
749 : zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
750 : vr - relative change of the drift velocity
751 : dzt - vd*dt
752 : dr = zz0-s*z
753 : ..
754 : ==>
755 : zs ~ z - s*vr*(z0-s*z)+s*dzt
756 : --------------------------------
757 : 1. Correction function vr constant:
758 :
759 :
760 : dz = zs-z = -s*vr *(z0-s*z)+s*dzt
761 : dzs/dl = dz/dl +s*s*vr*dz/dl
762 : d(dz/dl) = vr*dz/dl
763 : */
764 :
765 : const Float_t kZCut = 240; // remove the closest laser beam
766 : const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
767 : const Float_t kDistCut = 3; // distance sigma cut
768 : const Float_t kDistCutAbs = 0.25;
769 : const Float_t kMinClusters = 60; // minimal amount of the clusters
770 : const Float_t kMinSignal = 16; // minimal mean height of the signal
771 : const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
772 0 : static TLinearFitter fdriftA(3,"hyp2");
773 0 : static TLinearFitter fdriftC(3,"hyp2");
774 0 : static TLinearFitter fdriftAC(4,"hyp3");
775 0 : TVectorD fitA(3),fitC(3),fitAC(4);
776 :
777 0 : AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
778 0 : AliTPCParam * tpcparam = calib->GetParameters();
779 : //
780 0 : for (Int_t id=0; id<336; id++) fFitZ[id]=0;
781 :
782 0 : Float_t chi2A = 10;
783 0 : Float_t chi2C = 10;
784 0 : Float_t chi2AC = 10;
785 0 : Int_t npointsA=0;
786 0 : Int_t npointsC=0;
787 0 : Int_t npointsAC=0;
788 :
789 :
790 0 : for (Int_t iter=0; iter<3; iter++){
791 0 : fdriftA.ClearPoints();
792 0 : fdriftC.ClearPoints();
793 0 : fdriftAC.ClearPoints();
794 : //
795 0 : for (Int_t id=0; id<336; id++){
796 0 : if (!fTracksEsdParam.At(id)) continue;
797 0 : if (!AcceptLaser(id)) continue;
798 0 : if ( fClusterSatur[id]>kSaturCut) continue;
799 0 : if ( fClusterCounter[id]<kMinClusters) continue;
800 0 : AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
801 0 : if (track->GetTPCsignal()<kMinSignal) continue;
802 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
803 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
804 :
805 0 : Double_t xyz[3];
806 0 : Double_t pxyz[3];
807 0 : Double_t lxyz[3];
808 0 : Double_t lpxyz[3];
809 0 : param->GetXYZ(xyz);
810 0 : param->GetPxPyPz(pxyz);
811 0 : ltrp->GetXYZ(lxyz);
812 0 : ltrp->GetPxPyPz(lpxyz);
813 0 : if (TMath::Abs(lxyz[2])>kZCut) continue;
814 0 : Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
815 0 : if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
816 0 : if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
817 0 : if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
818 :
819 : // drift distance
820 0 : Double_t zlength = tpcparam->GetZLength(0);
821 0 : Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
822 0 : Double_t mdrift = zlength-TMath::Abs(xyz[2]);
823 0 : Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
824 0 : if (ltrp->GetSide()==0){
825 0 : fdriftA.AddPoint(xxx,mdrift,1);
826 : }else{
827 0 : fdriftC.AddPoint(xxx,mdrift,1);
828 : }
829 0 : Double_t xxx2[3] = { static_cast<Double_t>(ltrp->GetSide()), ldrift,lxyz[1]*ldrift/(zlength*250.)};
830 0 : fdriftAC.AddPoint(xxx2,mdrift,1);
831 0 : }
832 : //
833 0 : if (fdriftA.GetNpoints()>10){
834 : //
835 0 : fdriftA.Eval();
836 0 : if (iter==0) fdriftA.EvalRobust(0.7);
837 0 : else fdriftA.EvalRobust(0.8);
838 0 : fdriftA.GetParameters(fitA);
839 0 : npointsA= fdriftA.GetNpoints();
840 0 : chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
841 0 : if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
842 0 : if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
843 0 : (*fFitAside)[0] = fitA[0];
844 0 : (*fFitAside)[1] = fitA[1];
845 0 : (*fFitAside)[2] = fitA[2];
846 0 : (*fFitAside)[3] = fdriftA.GetNpoints();
847 0 : (*fFitAside)[4] = chi2A;
848 0 : }
849 : }
850 0 : if (fdriftC.GetNpoints()>10){
851 0 : fdriftC.Eval();
852 0 : if (iter==0) fdriftC.EvalRobust(0.7);
853 0 : else fdriftC.EvalRobust(0.8);
854 : //
855 0 : fdriftC.GetParameters(fitC);
856 0 : npointsC= fdriftC.GetNpoints();
857 0 : chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
858 0 : if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
859 0 : if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
860 0 : (*fFitCside)[0] = fitC[0];
861 0 : (*fFitCside)[1] = fitC[1];
862 0 : (*fFitCside)[2] = fitC[2];
863 0 : (*fFitCside)[3] = fdriftC.GetNpoints();
864 0 : (*fFitCside)[4] = chi2C;
865 0 : }
866 : }
867 :
868 0 : if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
869 0 : fdriftAC.Eval();
870 0 : if (iter==0) fdriftAC.EvalRobust(0.7);
871 0 : else fdriftAC.EvalRobust(0.8);
872 : //
873 0 : fdriftAC.GetParameters(fitAC);
874 0 : npointsAC= fdriftAC.GetNpoints();
875 0 : chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
876 0 : if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
877 : }
878 :
879 0 : for (Int_t id=0; id<336; id++){
880 0 : if (!fTracksEsdParam.At(id)) continue;
881 : //
882 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
883 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
884 0 : Double_t xyz[3];
885 0 : Double_t pxyz[3];
886 0 : Double_t lxyz[3];
887 0 : Double_t lpxyz[3];
888 0 : param->GetXYZ(xyz);
889 0 : param->GetPxPyPz(pxyz);
890 0 : ltrp->GetXYZ(lxyz);
891 0 : ltrp->GetPxPyPz(lpxyz);
892 0 : Double_t zlength = tpcparam->GetZLength(0);
893 0 : Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
894 0 : Double_t mdrift = zlength-TMath::Abs(xyz[2]);
895 :
896 : Float_t fz =0;
897 0 : if (ltrp->GetSide()==0){
898 0 : fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
899 0 : }else{
900 0 : fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
901 : }
902 0 : if (npointsAC>10){
903 0 : fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
904 0 : }
905 0 : fFitZ[id]=mdrift-fz;
906 0 : }
907 0 : if (fStreamLevel>0){
908 0 : TTreeSRedirector *cstream = GetDebugStreamer();
909 0 : TTimeStamp tstamp(fTime);
910 0 : Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
911 0 : Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
912 0 : Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
913 0 : Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
914 0 : Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
915 0 : Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
916 0 : TVectorD vecGoofie(20);
917 0 : AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
918 0 : if (goofieArray)
919 0 : for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
920 0 : AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
921 0 : if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
922 0 : }
923 :
924 0 : if (cstream){
925 0 : (*cstream)<<"driftv"<<
926 0 : "run="<<fRun<< // run number
927 0 : "event="<<fEvent<< // event number
928 0 : "time="<<fTime<< // time stamp of event
929 0 : "trigger="<<fTrigger<< // trigger
930 0 : "mag="<<fMagF<< // magnetic field
931 : // Environment values
932 0 : "press0="<<valuePressure0<<
933 0 : "press1="<<valuePressure1<<
934 0 : "pt0="<<ptrelative0<<
935 0 : "pt1="<<ptrelative1<<
936 0 : "temp0="<<temp0<<
937 0 : "temp1="<<temp1<<
938 0 : "vecGoofie.="<<&vecGoofie<<
939 : //
940 : //
941 0 : "iter="<<iter<<
942 0 : "driftA.="<<fFitAside<<
943 0 : "driftC.="<<fFitCside<<
944 0 : "driftAC.="<<fFitACside<<
945 0 : "chi2A="<<chi2A<<
946 0 : "chi2C="<<chi2C<<
947 0 : "chi2AC="<<chi2AC<<
948 0 : "nA="<<npointsA<<
949 0 : "nC="<<npointsC<<
950 0 : "nAC="<<npointsAC<<
951 : "\n";
952 : }
953 0 : }
954 : }
955 0 : }
956 : Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
957 : //
958 : // Fit corrections to the drift velocity - linear approximation in the z and global y
959 : //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
960 : //
961 : // Source of outlyers :
962 : // 0. Track in the saturation - postpeak
963 : // 1. gating grid close the part of the signal for first bundle
964 :
965 : // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
966 : // 1. Robust fit is used in the itteration number 0
967 : // only fraction of laser used
968 : // 2. Only the tracks close to the fit used in the second itteration
969 : /*
970 : Formulas:
971 :
972 : z = s* (z0 - vd*(t-t0))
973 :
974 : s - side -1 and +1
975 : t0 - time 0
976 : vd - nominal drift velocity
977 : zs - miscalibrated position
978 :
979 : zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
980 : vr - relative change of the drift velocity
981 : dzt - vd*dt
982 : dr = zz0-s*z
983 : ..
984 : ==>
985 : zs ~ z - s*vr*(z0-s*z)+s*dzt
986 : --------------------------------
987 : 1. Correction function vr constant:
988 :
989 :
990 : dz = zs-z = -s*vr *(z0-s*z)+s*dzt
991 : dzs/dl = dz/dl +s*s*vr*dz/dl
992 : d(dz/dl) = vr*dz/dl
993 : */
994 : const Int_t knLaser = 336; //n laser tracks
995 : const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
996 :
997 : const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
998 : const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
999 : const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
1000 : const Float_t kMinClusters = 40.; // minimal amount of the clusters
1001 : const Float_t kMinSignal = 2.5; // minimal mean height of the signal
1002 : const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
1003 : //
1004 0 : static TLinearFitter fdriftA(3,"hyp2");
1005 0 : static TLinearFitter fdriftC(3,"hyp2");
1006 0 : static TLinearFitter fdriftAC(4,"hyp3");
1007 0 : TVectorD fitA(3),fitC(3),fitAC(4);
1008 :
1009 0 : AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
1010 0 : AliTPCParam * tpcparam = calib->GetParameters();
1011 : //
1012 : // reset old data
1013 : //
1014 0 : for (Int_t id=0; id<336; id++) fFitZ[id]=0;
1015 0 : if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1016 0 : if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1017 0 : for (Int_t i=0;i<5; i++){
1018 0 : (*fFitCside)[i]=0;
1019 0 : (*fFitAside)[i]=0;
1020 : }
1021 : //
1022 : //
1023 0 : Float_t chi2A = 10;
1024 0 : Float_t chi2C = 10;
1025 0 : Float_t chi2AC = 10;
1026 0 : Int_t npointsA=0;
1027 0 : Int_t npointsC=0;
1028 0 : Int_t npointsAC=0;
1029 0 : Int_t nbA[4]={0,0,0,0};
1030 0 : Int_t nbC[4]={0,0,0,0};
1031 0 : TVectorD vecZM(336); // measured z potion of laser
1032 0 : TVectorD vecA(336); // accepted laser
1033 0 : TVectorD vecZF(336); // fitted position
1034 0 : TVectorD vecDz(336); // deltaZ
1035 0 : TVectorD vecZS(336); // surveyed position of laser
1036 : // additional variable to cut
1037 0 : TVectorD vecdEdx(336); // dEdx
1038 0 : TVectorD vecSy(336); // shape y
1039 0 : TVectorD vecSz(336); // shape z
1040 : //
1041 : //
1042 0 : for (Int_t id=0; id<336; id++){
1043 : Int_t reject=0;
1044 : AliTPCLaserTrack *ltrp =
1045 0 : (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1046 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1047 0 : AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1048 0 : vecZM(id)= (param==0) ? 0:param->GetZ();
1049 0 : vecZS(id)= ltrp->GetZ();
1050 0 : vecDz(id)= 0;
1051 0 : vecA(id)=1;
1052 0 : vecdEdx(id)=(seed)?seed->GetdEdx():0;
1053 0 : vecSy(id) =(seed)?seed->CookShape(1):0;
1054 0 : vecSz(id) =(seed)?seed->CookShape(2):0;
1055 : //
1056 0 : fFitZ[id]=0;
1057 0 : if (!fTracksEsdParam.At(id)) reject|=1;
1058 0 : if (!param) continue;
1059 0 : if (!AcceptLaser(id)) reject|=2;
1060 0 : if ( fClusterSatur[id]>kSaturCut) reject|=4;
1061 0 : if ( fClusterCounter[id]<kMinClusters) reject|=8;
1062 0 : vecA(id)=reject;
1063 0 : if (reject>0) continue;
1064 0 : if (ltrp->GetSide()==0){
1065 0 : npointsA++;
1066 0 : nbA[ltrp->GetBundle()]++;
1067 0 : }
1068 0 : if (ltrp->GetSide()==1){
1069 0 : npointsC++;
1070 0 : nbC[ltrp->GetBundle()]++;
1071 0 : }
1072 0 : }
1073 : //
1074 : // reject "bad events"
1075 : //
1076 : Bool_t isOK=kTRUE;
1077 : Int_t naA=0;
1078 : Int_t naC=0;
1079 0 : if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1080 0 : isOK=kFALSE;
1081 0 : for (Int_t i=0;i<4;i++){
1082 : //count accepted for all layers
1083 0 : if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1084 0 : if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1085 : }
1086 0 : if (naA<3 &&naC<3) isOK=kFALSE;
1087 0 : if (isOK==kFALSE) return kFALSE;
1088 :
1089 : //
1090 : //
1091 : //
1092 0 : for (Int_t iter=0; iter<2; iter++){
1093 0 : fdriftA.ClearPoints();
1094 0 : fdriftC.ClearPoints();
1095 0 : fdriftAC.ClearPoints();
1096 0 : npointsA=0; npointsC=0; npointsAC=0;
1097 : //
1098 0 : for (Int_t id=0; id<336; id++){
1099 0 : if (!fTracksEsdParam.At(id)) continue;
1100 0 : if (!AcceptLaser(id)) continue;
1101 0 : if ( fClusterSatur[id]>kSaturCut) continue;
1102 0 : if ( fClusterCounter[id]<kMinClusters) continue;
1103 0 : AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1104 0 : if (track->GetTPCsignal()<kMinSignal) continue;
1105 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1106 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1107 0 : Double_t xyz[3];
1108 0 : Double_t pxyz[3];
1109 0 : Double_t lxyz[3];
1110 0 : Double_t lpxyz[3];
1111 0 : param->GetXYZ(xyz);
1112 0 : param->GetPxPyPz(pxyz);
1113 0 : ltrp->GetXYZ(lxyz);
1114 0 : ltrp->GetPxPyPz(lpxyz);
1115 0 : Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1116 : //if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1117 0 : if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1118 0 : if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1119 :
1120 : // // drift distance
1121 : // Double_t zlength = tpcparam->GetZLength(0);
1122 : // Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1123 : // Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1124 : //
1125 0 : Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1126 0 : Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1127 0 : Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1128 :
1129 :
1130 0 : Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1131 0 : if (iter==0 &<rp->GetBundle()==0) continue;
1132 : // skip bundle 0 - can be wrong in the 0 iteration
1133 0 : if (ltrp->GetSide()==0){
1134 0 : fdriftA.AddPoint(xxx,mdrift,1);
1135 : }else{
1136 0 : fdriftC.AddPoint(xxx,mdrift,1);
1137 : }
1138 0 : Double_t xxx2[3] = { static_cast<Double_t>(ltrp->GetSide()), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1139 0 : fdriftAC.AddPoint(xxx2,mdrift,1);
1140 0 : }
1141 : //
1142 0 : if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1143 : //
1144 0 : fdriftA.Eval();
1145 : //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
1146 0 : npointsA= fdriftA.GetNpoints();
1147 0 : chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1148 0 : fdriftA.EvalRobust(kFraction[iter]);
1149 0 : fdriftA.GetParameters(fitA);
1150 0 : if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1151 0 : if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1152 0 : (*fFitAside)[0] = fitA[0];
1153 0 : (*fFitAside)[1] = fitA[1];
1154 0 : (*fFitAside)[2] = fitA[2];
1155 0 : (*fFitAside)[3] = fdriftA.GetNpoints();
1156 0 : (*fFitAside)[4] = chi2A;
1157 0 : }
1158 : }
1159 0 : if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1160 0 : fdriftC.Eval();
1161 : //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
1162 0 : npointsC= fdriftC.GetNpoints();
1163 0 : chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1164 0 : fdriftC.EvalRobust(kFraction[iter]);
1165 0 : fdriftC.GetParameters(fitC);
1166 0 : if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1167 0 : if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1168 0 : (*fFitCside)[0] = fitC[0];
1169 0 : (*fFitCside)[1] = fitC[1];
1170 0 : (*fFitCside)[2] = fitC[2];
1171 0 : (*fFitCside)[3] = fdriftC.GetNpoints();
1172 0 : (*fFitCside)[4] = chi2C;
1173 0 : }
1174 : }
1175 :
1176 0 : if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1177 0 : fdriftAC.Eval();
1178 : //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
1179 0 : npointsAC= fdriftAC.GetNpoints();
1180 0 : chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1181 0 : fdriftAC.EvalRobust(kFraction[iter]);
1182 0 : fdriftAC.GetParameters(fitAC);
1183 0 : if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1184 0 : (*fFitACside)[0] = fitAC[0];
1185 0 : (*fFitACside)[1] = fitAC[1];
1186 0 : (*fFitACside)[2] = fitAC[2];
1187 0 : (*fFitACside)[3] = fdriftAC.GetNpoints();
1188 0 : (*fFitACside)[4] = chi2AC;
1189 0 : }
1190 :
1191 0 : for (Int_t id=0; id<336; id++){
1192 0 : if (!fTracksEsdParam.At(id)) continue;
1193 : //
1194 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1195 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1196 0 : Double_t xyz[3];
1197 0 : Double_t pxyz[3];
1198 0 : Double_t lxyz[3];
1199 0 : Double_t lpxyz[3];
1200 0 : param->GetXYZ(xyz);
1201 0 : param->GetPxPyPz(pxyz);
1202 0 : ltrp->GetXYZ(lxyz);
1203 0 : ltrp->GetPxPyPz(lpxyz);
1204 : //Double_t zlength = tpcparam->GetZLength(0);
1205 : //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1206 : //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1207 0 : Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1208 0 : Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1209 0 : Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1210 :
1211 :
1212 : Float_t fz =0;
1213 0 : if (ltrp->GetSide()==0){
1214 0 : fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1215 0 : }else{
1216 0 : fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1217 : }
1218 0 : if (npointsAC>10){
1219 : //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1220 : }
1221 0 : fFitZ[id]=mdrift-fz;
1222 0 : vecZF[id]=fz;
1223 0 : vecDz[id]=mdrift-fz;
1224 0 : }
1225 0 : if (fStreamLevel>0){
1226 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1227 0 : TTimeStamp tstamp(fTime);
1228 0 : Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1229 0 : Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1230 0 : Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1231 0 : Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1232 0 : Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1233 0 : Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1234 0 : TVectorD vecGoofie(20);
1235 0 : AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1236 0 : if (goofieArray)
1237 0 : for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1238 0 : AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1239 0 : if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1240 0 : }
1241 :
1242 0 : if (cstream){
1243 0 : (*cstream)<<"driftvN"<<
1244 0 : "run="<<fRun<< // run number
1245 0 : "event="<<fEvent<< // event number
1246 0 : "time="<<fTime<< // time stamp of event
1247 0 : "trigger="<<fTrigger<< // trigger
1248 0 : "mag="<<fMagF<< // magnetic field
1249 : // Environment values
1250 0 : "press0="<<valuePressure0<<
1251 0 : "press1="<<valuePressure1<<
1252 0 : "pt0="<<ptrelative0<<
1253 0 : "pt1="<<ptrelative1<<
1254 0 : "temp0="<<temp0<<
1255 0 : "temp1="<<temp1<<
1256 0 : "vecGoofie.="<<&vecGoofie<<
1257 : //
1258 : //
1259 0 : "vecZM.="<<&vecZM<< // measured z position
1260 0 : "vecZS.="<<&vecZS<< // surveyed z position
1261 0 : "vecZF.="<<&vecZF<< // fitted z position
1262 0 : "vecDz.="<<&vecDz<< // fitted z position
1263 0 : "vecA.="<<&vecA<< // accept laser flag
1264 0 : "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
1265 0 : "vecSy.="<<&vecSy<< // shape y - to cut on
1266 0 : "vecSz.="<<&vecSz<< // shape z - to cut on
1267 : //
1268 0 : "iter="<<iter<<
1269 0 : "driftA.="<<fFitAside<<
1270 0 : "driftC.="<<fFitCside<<
1271 0 : "driftAC.="<<fFitACside<<
1272 0 : "chi2A="<<chi2A<<
1273 0 : "chi2C="<<chi2C<<
1274 0 : "chi2AC="<<chi2AC<<
1275 0 : "nA="<<npointsA<<
1276 0 : "nC="<<npointsC<<
1277 0 : "nAC="<<npointsAC<<
1278 : "\n";
1279 : /*
1280 : //
1281 : variables to check in debug mode:
1282 : //
1283 : chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1284 : chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1285 : chainDriftN->SetAlias("driftF","vecZF.fElements");
1286 : chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
1287 : TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1288 : TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1289 :
1290 : */
1291 : }
1292 0 : }
1293 : }
1294 0 : return kTRUE;
1295 0 : }
1296 :
1297 : Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1298 : //
1299 : // get distance between mirror and laser track
1300 : //
1301 : //
1302 0 : Double_t xyz[3];
1303 0 : Double_t lxyz[3];
1304 0 : param->GetXYZ(xyz);
1305 0 : ltrp->GetXYZ(lxyz);
1306 : //
1307 : //
1308 : Double_t dist = 0;
1309 : //radial distance
1310 0 : dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1311 : //
1312 : // z distance
1313 : // apply drift correction if already exist
1314 : //
1315 : Float_t dz = 0;
1316 0 : if (ltrp->GetSide()==0){
1317 0 : if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1318 : }else{
1319 0 : if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1320 : }
1321 0 : if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1322 0 : dist+=TMath::Abs(dz);
1323 : //
1324 : // phi dist - divergence on 50 cm
1325 : //
1326 0 : dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1327 0 : return dist;
1328 0 : }
1329 :
1330 :
1331 : Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
1332 : //
1333 : //
1334 : //
1335 : /*
1336 : TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1337 : TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1338 : TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1339 : TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1340 :
1341 : TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1342 : */
1343 0 : AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1344 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1345 : //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1346 0 : Double_t xyz[3];
1347 0 : Double_t lxyz[3];
1348 0 : param->GetXYZ(xyz);
1349 0 : ltrp->GetXYZ(lxyz);
1350 0 : if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1351 0 : if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
1352 0 : if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
1353 0 : if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
1354 : //
1355 : //
1356 :
1357 0 : return kTRUE;
1358 0 : }
1359 :
1360 : Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1361 : //
1362 : // Find corresponding mirror
1363 : // add the corresponding tracks
1364 :
1365 0 : if (!track->GetOuterParam()) return -1;
1366 :
1367 : Float_t kRadius0 = 252;
1368 : Float_t kRadius = 254.2;
1369 : Int_t countercl=0;
1370 : Float_t counterSatur=0;
1371 : Int_t csideA =0;
1372 : Int_t csideC =0;
1373 0 : for (Int_t irow=158;irow>-1;--irow) {
1374 0 : AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1375 0 : if (!c) continue;
1376 0 : Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1377 0 : Double_t pedgeX = TMath::Min((irow)*0.75, (kMaxRow-irow)*1.5);
1378 0 : if (pedgeY<3) continue;
1379 0 : if (pedgeX<3) continue;
1380 0 : countercl++;
1381 0 : if (c->GetDetector()%36<18) csideA++;
1382 0 : if (c->GetDetector()%36>=18) csideC++;
1383 0 : if (c->GetMax()>900) counterSatur++;
1384 0 : }
1385 0 : counterSatur/=(countercl+1);
1386 : //
1387 : //
1388 : //
1389 0 : if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
1390 :
1391 : Int_t side= 0;
1392 0 : if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1393 :
1394 :
1395 0 : AliExternalTrackParam param(*(track->GetOuterParam()));
1396 0 : AliTracker::PropagateTrackTo(¶m,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
1397 0 : AliTracker::PropagateTrackTo(¶m,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
1398 0 : AliTPCLaserTrack ltr;
1399 : AliTPCLaserTrack *ltrp=0x0;
1400 : // AliTPCLaserTrack *ltrpjw=0x0;
1401 : //
1402 0 : Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
1403 : // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
1404 : //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1405 :
1406 0 : if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1407 0 : ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1408 : else
1409 : ltrp=<r;
1410 :
1411 0 : if (id<0) return -1;
1412 0 : if (ltrp->GetSide()!=side) return -1;
1413 0 : fCounter[id]++;
1414 : //
1415 : //
1416 : //
1417 : //
1418 0 : if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1419 : //
1420 : //
1421 0 : Float_t radius=TMath::Abs(ltrp->GetX());
1422 0 : param.Rotate(ltrp->GetAlpha());
1423 0 : AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
1424 : //
1425 0 : if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1426 : Bool_t accept=kTRUE;
1427 : //
1428 : // choose closer track
1429 : //
1430 0 : AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1431 0 : if (param0){
1432 0 : Float_t dist0=GetDistance(param0,ltrp);
1433 0 : Float_t dist1=GetDistance(¶m,ltrp);
1434 0 : if (dist0<dist1) accept=kFALSE;
1435 0 : }
1436 :
1437 0 : if (accept){
1438 0 : fClusterCounter[id]=countercl;
1439 0 : fTracksEsdParam.AddAt(param.Clone(),id);
1440 0 : fTracksEsd.AddAt(track,id);
1441 0 : fTracksTPC.AddAt(seed,id);
1442 : }
1443 : return id;
1444 0 : }
1445 :
1446 :
1447 :
1448 : void AliTPCcalibLaser::DumpLaser(Int_t id) {
1449 : //
1450 : // Dump Laser info to the tree
1451 : //
1452 0 : AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
1453 0 : AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1454 0 : AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1455 : //
1456 : // Fast laser ID
1457 : //
1458 0 : Double_t xyz[3];
1459 0 : Double_t pxyz[3];
1460 0 : Double_t lxyz[3];
1461 0 : Double_t lpxyz[3];
1462 0 : param->GetXYZ(xyz);
1463 0 : param->GetPxPyPz(pxyz);
1464 0 : ltrp->GetXYZ(lxyz);
1465 0 : ltrp->GetPxPyPz(lpxyz);
1466 0 : Float_t dist3D = GetDistance(param,ltrp);
1467 0 : Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1468 0 : Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1469 :
1470 :
1471 0 : if (fStreamLevel>0){
1472 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1473 0 : Int_t time = fESD->GetTimeStamp();
1474 0 : Bool_t accept = AcceptLaser(id);
1475 0 : if (cstream){
1476 0 : (*cstream)<<"Track"<<
1477 : //
1478 0 : "run="<<fRun<< // run number
1479 0 : "event="<<fEvent<< // event number
1480 0 : "time="<<fTime<< // time stamp of event
1481 0 : "trigger="<<fTrigger<< // trigger
1482 0 : "mag="<<fMagF<< // magnetic field
1483 :
1484 0 : "id="<<id<<
1485 0 : "accept="<<accept<<
1486 0 : "driftA.="<<fFitAside<<
1487 0 : "driftC.="<<fFitCside<<
1488 0 : "time="<<time<<
1489 0 : "dist3D="<<dist3D<<
1490 0 : "dist0="<<dist0<<
1491 0 : "distphi="<<distphi<<
1492 : //
1493 : //
1494 0 : "counter="<<fCounter[id]<<
1495 0 : "clcounter="<<fClusterCounter[id]<<
1496 0 : "satur="<<fClusterSatur[id]<<
1497 0 : "fz="<<fFitZ[id]<<
1498 : //
1499 0 : "LTr.="<<ltrp<<
1500 0 : "Esd.="<<track<<
1501 0 : "Tr.="<<param<<
1502 0 : "x0="<<xyz[0]<<
1503 0 : "x1="<<xyz[1]<<
1504 0 : "x2="<<xyz[2]<<
1505 0 : "px0="<<pxyz[0]<<
1506 0 : "px1="<<pxyz[1]<<
1507 0 : "px2="<<pxyz[2]<<
1508 : //
1509 0 : "lx0="<<lxyz[0]<<
1510 0 : "lx1="<<lxyz[1]<<
1511 0 : "lx2="<<lxyz[2]<<
1512 0 : "lpx0="<<lpxyz[0]<<
1513 0 : "lpx1="<<lpxyz[1]<<
1514 0 : "lpx2="<<lpxyz[2]<<
1515 : "\n";
1516 0 : }
1517 0 : }
1518 0 : }
1519 :
1520 : void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1521 : //
1522 : // Refit the track with different tracklet models:
1523 : // 1. Per ROC using the kalman filter, different edge cuts
1524 : // 2. Per ROC linear in y and z
1525 : // 3. Per ROC quadratic in y and z
1526 : // 4. Per track offset for each sector, linear for each sector, common quadratic
1527 : // store x, y, z information for all models and the cluster to calculate the residuals
1528 : //
1529 :
1530 : // The clusters which do not fulfill given criteria are skipped
1531 : //
1532 : // Cluters removed from fit
1533 : // 1. Extended shape > kShapeCut
1534 : // 2. In saturation Max > kMax
1535 : // 3. Distance to edge < cutEdge
1536 : //
1537 : // Clusters not used for the calibration:
1538 : //
1539 : // 1. Extended shape > kShapeCut
1540 : // 2. In saturation Max > kMax
1541 :
1542 :
1543 0 : AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
1544 0 : AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1545 0 : AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1546 :
1547 0 : AliTPCclusterMI dummyCl;
1548 :
1549 : //two tracklets
1550 : Int_t kMaxTracklets=2;
1551 : Float_t kShapeCut = 1.3;
1552 : Float_t kRatioCut = 2.;
1553 :
1554 : Float_t kMax = 900.;
1555 :
1556 :
1557 : //=============================================//
1558 : // Linear Fitters for the Different Approaches //
1559 : //=============================================//
1560 : //linear fit model in y and z; inner - outer sector, combined with offset
1561 0 : static TLinearFitter fy1I(2,"hyp1");
1562 0 : static TLinearFitter fy1O(2,"hyp1");
1563 0 : static TLinearFitter fz1I(2,"hyp1");
1564 0 : static TLinearFitter fz1O(2,"hyp1");
1565 0 : static TLinearFitter fy1IO(3,"hyp2");
1566 0 : static TLinearFitter fz1IO(3,"hyp2");
1567 : //quadratic fit model in y and z; inner - sector
1568 0 : static TLinearFitter fy2I(3,"hyp2");
1569 0 : static TLinearFitter fy2O(3,"hyp2");
1570 0 : static TLinearFitter fz2I(3,"hyp2");
1571 0 : static TLinearFitter fz2O(3,"hyp2");
1572 : //common quadratic fit for IROC and OROC in y and z
1573 0 : static TLinearFitter fy4(5,"hyp4");
1574 0 : static TLinearFitter fz4(5,"hyp4");
1575 :
1576 :
1577 : //set standard cuts
1578 0 : if ( fNcuts==0 ){
1579 0 : fNcuts=1;
1580 0 : fEdgeXcuts[0]=4;
1581 0 : fEdgeYcuts[0]=3;
1582 0 : fNClCuts[0]=20;
1583 0 : }
1584 : //=============================//
1585 : // Loop over all Tracklet Cuts //
1586 : //=============================//
1587 0 : for (Int_t icut=0; icut<fNcuts; icut++){
1588 0 : Float_t xinMin = 2500, xinMax=-90;
1589 0 : Float_t xoutMin=2500, xoutMax=-90;
1590 0 : Float_t msigmaYIn=0;
1591 0 : Float_t msigmaYOut=0;
1592 0 : Float_t mqratioIn=0;
1593 0 : Float_t mqratioOut=0;
1594 :
1595 :
1596 0 : AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1597 : //cut parameters
1598 0 : Double_t edgeCutX = fEdgeXcuts[icut];
1599 0 : Double_t edgeCutY = fEdgeYcuts[icut];
1600 0 : Int_t nclCut = (Int_t)fNClCuts[icut];
1601 : //=========================//
1602 : // Parameters to calculate //
1603 : //=========================//
1604 : //fit parameter inner, outer tracklet and combined fit
1605 0 : TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1606 0 : TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1607 : //
1608 0 : TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1609 0 : TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1610 0 : TVectorD vecy4res(5),vecz4res(5);
1611 0 : TVectorD vecy1resIO(3),vecz1resIO(3);
1612 : // cluster and track positions for each row - used for residuals
1613 0 : TVectorD vecgX(kMaxRow); // global X
1614 0 : TVectorD vecgY(kMaxRow); // global Y
1615 0 : TVectorD vecgZ(kMaxRow); // global Z
1616 :
1617 0 : TVectorD vecX(kMaxRow); // x is the same for all (row center)
1618 0 : TVectorD vecYkalman(kMaxRow); // y from kalman fit
1619 0 : TVectorD vecZkalman(kMaxRow); // z from kalman fit
1620 0 : TVectorD vecY1(kMaxRow); // y from pol1 fit per ROC
1621 0 : TVectorD vecZ1(kMaxRow); // z from pol1 fit per ROC
1622 0 : TVectorD vecY1IO(kMaxRow); // y from pol1 fit per ROC
1623 0 : TVectorD vecZ1IO(kMaxRow); // z from pol1 fit per ROC
1624 0 : TVectorD vecY2(kMaxRow); // y from pol2 fit per ROC
1625 0 : TVectorD vecZ2(kMaxRow); // z from pol2 fit per ROC
1626 0 : TVectorD vecY4(kMaxRow); // y from sector fit
1627 0 : TVectorD vecZ4(kMaxRow); // z from sector fit
1628 0 : TVectorD vecClY(kMaxRow); // y cluster position
1629 0 : TVectorD vecClZ(kMaxRow); // z cluster position
1630 0 : TVectorD vecSec(kMaxRow); // sector for each row
1631 0 : TVectorD isReject(kMaxRow); // flag - cluster to be rejected
1632 : //chi2 of fits
1633 0 : Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
1634 0 : Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
1635 0 : Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
1636 0 : Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
1637 0 : Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
1638 0 : Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
1639 0 : Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
1640 0 : Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
1641 0 : Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
1642 0 : Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
1643 0 : Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
1644 0 : Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
1645 : //more
1646 0 : Int_t innerSector = -1; // number of inner sector
1647 0 : Int_t outerSector = -1; // number of outer sector
1648 0 : Int_t nclI=0; // number of clusters (inner)
1649 0 : Int_t nclO=0; // number of clusters (outer)
1650 : //=================================================//
1651 : // Perform the Kalman Fit using the Tracklet Class //
1652 : //=================================================//
1653 0 : AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1654 0 : TObjArray tracklets=
1655 0 : AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1656 0 : kFALSE,nclCut,kMaxTracklets);
1657 : // tracklet pointers
1658 0 : AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1659 0 : AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1660 : AliTPCTracklet *tr=0x0;
1661 0 : AliTPCTracklet dummy;
1662 : //continue if we didn't find a tracklet
1663 0 : if ( !trInner && !trOuter ) continue;
1664 : //================================================================================//
1665 : // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1666 : //================================================================================//
1667 0 : if ( trInner && trOuter ){
1668 0 : if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1669 0 : if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1670 : tr = trInner;
1671 : trInner=trOuter;
1672 : trOuter=tr;
1673 0 : AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1674 : }
1675 : } else {
1676 0 : if ( trInner ){
1677 0 : if ( !trInner->GetInner() ) continue;
1678 : trOuter=&dummy;
1679 0 : if ( trInner->GetSector()>35 ){
1680 : trOuter=trInner;
1681 : trInner=&dummy;
1682 0 : }
1683 : } else { //trOuter
1684 0 : if ( !trOuter->GetInner() ) continue;
1685 : trInner=&dummy;
1686 0 : if ( trOuter->GetSector()<36 ){
1687 : trInner=trOuter;
1688 : trOuter=&dummy;
1689 0 : }
1690 : }
1691 : }
1692 0 : innerSector = trInner->GetSector();
1693 0 : if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1694 0 : outerSector = trOuter->GetSector();
1695 0 : if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1696 :
1697 : // array of clusters
1698 0 : TClonesArray arrCl("AliTPCclusterMI",kMaxRow);
1699 0 : arrCl.ExpandCreateFast(kMaxRow);
1700 : //=======================================//
1701 : // fill fitters with cluster information //
1702 : //=======================================//
1703 0 : AliDebug(3,"Filing Cluster Information");
1704 :
1705 : //
1706 : //
1707 :
1708 0 : for (Int_t irow=158;irow>-1;--irow) {
1709 0 : AliTPCclusterMI *c=track->GetClusterPointer(irow);
1710 0 : AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1711 0 : cl=dummyCl;
1712 0 : vecX[irow] = 0;
1713 0 : vecClY[irow] = 0;
1714 0 : vecClZ[irow] = 0;
1715 : Float_t meanY=0, sumY=0;
1716 0 : for (Int_t drow=-1;drow<=1;drow++) {
1717 0 : if (irow+drow<0) continue;
1718 0 : if (irow+drow>158) continue;
1719 0 : AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1720 0 : if (!ccurrent) continue;
1721 0 : Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1722 0 : if ( roc!=innerSector && roc!=outerSector ) continue;
1723 0 : if (ccurrent->GetX()<10) continue;
1724 0 : if (ccurrent->GetY()==0) continue;
1725 0 : meanY+=ccurrent->GetY();
1726 0 : sumY++;
1727 0 : }
1728 0 : if (sumY>0) meanY/=sumY;
1729 :
1730 : //
1731 0 : vecSec[irow]=-1;
1732 0 : if (!c) continue;
1733 0 : Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1734 0 : Double_t pedgeX = TMath::Min((irow)*0.75, (kMaxRow-irow)*1.5);
1735 :
1736 : //
1737 0 : Int_t roc = static_cast<Int_t>(c->GetDetector());
1738 0 : if ( roc!=innerSector && roc!=outerSector ) continue;
1739 0 : vecSec[irow]=roc;
1740 : //store clusters in clones array
1741 0 : cl=*c;
1742 : //
1743 0 : if (c->GetMax()<4) continue; // noise cluster?
1744 0 : if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
1745 : //cluster position
1746 0 : vecX[irow] = c->GetX();
1747 0 : vecClY[irow] = c->GetY();
1748 0 : vecClZ[irow] = c->GetZ();
1749 : //
1750 : // Float_t gxyz[3];
1751 : // c->GetGlobalXYZ(gxyz);
1752 : // vecgX[irow] = gxyz[0];
1753 : // vecgY[irow] = gxyz[1];
1754 : // vecgZ[irow] = gxyz[2];
1755 : //
1756 0 : Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1757 0 : Double_t y = vecClY[irow];
1758 0 : Double_t z = vecClZ[irow];
1759 : //
1760 0 : Double_t x2[2]={x,x*x}; //linear and parabolic parameters
1761 0 : Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1762 0 : Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
1763 0 : if ( roc == innerSector ) {
1764 0 : x4[0]=1; //offset inner - outer sector
1765 0 : x4[1]=x; //slope parameter inner sector
1766 0 : xIO[0]=1;
1767 0 : } else {
1768 0 : x4[2]=x; //slope parameter outer sector
1769 : }
1770 0 : x4[3]=x*x; //common parabolic shape
1771 0 : if (pedgeX < fEdgeXcuts[icut]) continue;
1772 0 : if (pedgeY < fEdgeYcuts[icut]) continue;
1773 0 : if (c->GetMax()>900) continue; // cluster in saturation
1774 : //
1775 0 : if ( roc==innerSector ){
1776 0 : fy1I.AddPoint(x2,y);
1777 0 : fz1I.AddPoint(x2,z);
1778 0 : fy2I.AddPoint(x2,y);
1779 0 : fz2I.AddPoint(x2,z);
1780 0 : ++nclI;
1781 0 : if (c->GetX()<xinMin) xinMin=c->GetX();
1782 0 : if (c->GetX()>xinMax) xinMax=c->GetX();
1783 :
1784 0 : msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1785 0 : mqratioIn +=c->GetMax()/c->GetQ();
1786 :
1787 0 : }
1788 0 : if ( roc==outerSector ){
1789 0 : fy1O.AddPoint(x2,y);
1790 0 : fz1O.AddPoint(x2,z);
1791 0 : fy2O.AddPoint(x2,y);
1792 0 : fz2O.AddPoint(x2,z);
1793 0 : ++nclO;
1794 0 : if (c->GetX()<xoutMin) xoutMin=c->GetX();
1795 0 : if (c->GetX()>xoutMax) xoutMax=c->GetX();
1796 0 : msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1797 0 : mqratioOut +=c->GetMax()/c->GetQ();
1798 :
1799 0 : }
1800 0 : fy4.AddPoint(x4,y);
1801 0 : fz4.AddPoint(x4,z);
1802 0 : fy1IO.AddPoint(xIO,y);
1803 0 : fz1IO.AddPoint(xIO,z);
1804 0 : }
1805 0 : if (nclI>0) {
1806 0 : msigmaYIn/=nclI;
1807 0 : mqratioIn/=nclI;
1808 0 : }
1809 0 : if (nclO>0) {
1810 0 : msigmaYOut/=nclO;
1811 0 : mqratioOut/=nclO;
1812 0 : }
1813 : //======================================//
1814 : // Evaluate and retrieve fit parameters //
1815 : //======================================//
1816 0 : AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1817 : //inner sector
1818 0 : if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1819 0 : fy1I.Eval();
1820 0 : fz1I.Eval();
1821 0 : fy2I.Eval();
1822 0 : fz2I.Eval();
1823 0 : fy1I.GetParameters(vecy1resInner);
1824 0 : fz1I.GetParameters(vecz1resInner);
1825 0 : fy2I.GetParameters(vecy2resInner);
1826 0 : fz2I.GetParameters(vecz2resInner);
1827 0 : chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1828 0 : chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1829 0 : chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1830 0 : chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1831 0 : }
1832 : //outer sector
1833 0 : if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1834 0 : fy1O.Eval();
1835 0 : fz1O.Eval();
1836 0 : fy2O.Eval();
1837 0 : fz2O.Eval();
1838 0 : fy1O.GetParameters(vecy1resOuter);
1839 0 : fz1O.GetParameters(vecz1resOuter);
1840 0 : fy2O.GetParameters(vecy2resOuter);
1841 0 : fz2O.GetParameters(vecz2resOuter);
1842 0 : chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1843 0 : chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1844 0 : chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1845 0 : chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1846 0 : }
1847 : //combined hyp4 fit
1848 0 : if ( innerSector>0 && outerSector>0 ){
1849 0 : if (fy4.GetNpoints()>0) {
1850 0 : fy4.Eval();
1851 0 : fy4.GetParameters(vecy4res);
1852 0 : chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1853 0 : }
1854 0 : if (fz4.GetNpoints()>0) {
1855 0 : fz4.Eval();
1856 0 : fz4.GetParameters(vecz4res);
1857 0 : chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1858 0 : }
1859 0 : if (fy1IO.GetNpoints()>0) {
1860 0 : fy1IO.Eval();
1861 0 : fy1IO.GetParameters(vecy1resIO);
1862 0 : chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1863 0 : }
1864 0 : if (fz1IO.GetNpoints()>0) {
1865 0 : fz1IO.Eval();
1866 0 : fz1IO.GetParameters(vecz1resIO);
1867 0 : chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1868 0 : }
1869 : }
1870 : //clear points
1871 0 : fy4.ClearPoints(); fz4.ClearPoints();
1872 0 : fy1I.ClearPoints(); fy1O.ClearPoints();
1873 0 : fz1I.ClearPoints(); fz1O.ClearPoints();
1874 0 : fy2I.ClearPoints(); fy2O.ClearPoints();
1875 0 : fz2I.ClearPoints(); fz2O.ClearPoints();
1876 0 : fy1IO.ClearPoints(); fz1IO.ClearPoints();
1877 : //==============================//
1878 : // calculate tracklet positions //
1879 : //==============================//
1880 0 : AliDebug(4,"Calculate tracklet positions");
1881 0 : for (Int_t irow=158;irow>-1;--irow) {
1882 0 : isReject[irow]=0;
1883 0 : AliTPCclusterMI *c=track->GetClusterPointer(irow);
1884 0 : if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1885 0 : isReject[irow]+=1;
1886 0 : }
1887 :
1888 0 : if (!c) { //no cluster
1889 0 : isReject[irow]+=2;
1890 0 : }else{
1891 0 : if (c->GetMax()>kMax){ //saturation
1892 0 : isReject[irow]+=4;
1893 0 : }
1894 0 : if ( vecSec[irow] == outerSector ) { //extended shape
1895 0 : if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
1896 0 : if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
1897 : }else{
1898 0 : if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
1899 0 : if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
1900 : }
1901 : }
1902 :
1903 :
1904 :
1905 0 : if ( vecSec[irow]==-1 ) continue; //no cluster info
1906 0 : if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1907 : tr=&dummy;
1908 0 : Double_t x=vecX[irow];
1909 0 : Double_t xref=x-133.4;
1910 : //
1911 : Double_t yoffInner=0;
1912 : Double_t zoffInner=0;
1913 : Double_t yoffInner1=0;
1914 : Double_t zoffInner1=0;
1915 : Double_t yslopeInner=0;
1916 : Double_t yslopeOuter=0;
1917 : Double_t zslopeInner=0;
1918 : Double_t zslopeOuter=0;
1919 : //positions of hyperplane fits
1920 0 : if ( vecSec[irow] == outerSector ) {
1921 : tr=trOuter;
1922 0 : vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1923 0 : vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1924 0 : vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1925 0 : vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1926 0 : yslopeOuter=vecy4res[3];
1927 0 : zslopeOuter=vecz4res[3];
1928 0 : } else {
1929 : tr=trInner;
1930 0 : vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1931 0 : vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1932 0 : vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1933 0 : vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1934 0 : yoffInner=vecy4res[1];
1935 0 : zoffInner=vecz4res[1];
1936 0 : yoffInner1=vecy1resIO[1];
1937 0 : zoffInner1=vecz1resIO[1];
1938 0 : yslopeInner=vecy4res[2];
1939 0 : zslopeInner=vecz4res[2];
1940 : }
1941 0 : vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1942 0 : vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1943 0 : vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1944 0 : vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1945 : //positions of kalman fits
1946 0 : Double_t gxyz[3],xyz[3];
1947 : AliExternalTrackParam *param = 0x0;
1948 : //
1949 0 : param=tr->GetInner();
1950 0 : if (param){
1951 0 : param->GetXYZ(gxyz);
1952 0 : Float_t bz = AliTracker::GetBz(gxyz);
1953 0 : param->GetYAt(x, bz, xyz[1]);
1954 0 : param->GetZAt(x, bz, xyz[2]);
1955 0 : vecYkalman[irow]=xyz[1];
1956 0 : vecZkalman[irow]=xyz[2];
1957 0 : }
1958 : //
1959 : //
1960 : //
1961 :
1962 0 : }
1963 : //=====================================================================//
1964 : // write results from the different tracklet fits with debug streamers //
1965 : //=====================================================================//
1966 0 : if (fStreamLevel>4){
1967 0 : TTreeSRedirector *cstream = GetDebugStreamer();
1968 0 : if (cstream){
1969 0 : Float_t dedx = track->GetdEdx();
1970 0 : (*cstream)<<"FitModels"<<
1971 0 : "run="<<fRun<< // run number
1972 0 : "event="<<fEvent<< // event number
1973 0 : "time="<<fTime<< // time stamp of event
1974 0 : "trigger="<<fTrigger<< // trigger
1975 0 : "mag="<<fMagF<< // magnetic field
1976 : //
1977 0 : "cutNr=" << icut <<
1978 0 : "edgeCutX=" << edgeCutX <<
1979 0 : "edgeCutY=" << edgeCutY <<
1980 0 : "nclCut=" << nclCut <<
1981 0 : "innerSector="<< innerSector <<
1982 0 : "outerSector="<< outerSector <<
1983 0 : "dEdx=" << dedx <<
1984 0 : "LTr.=" << ltrp <<
1985 0 : "Tr.=" << extparam <<
1986 0 : "yPol1In.=" << &vecy1resInner <<
1987 0 : "zPol1In.=" << &vecz1resInner <<
1988 0 : "yPol1InOut.="<< &vecy1resIO <<
1989 0 : "zPol1InOut.="<< &vecz1resIO <<
1990 0 : "yPol2In.=" << &vecy2resInner <<
1991 0 : "zPol2In.=" << &vecz2resInner <<
1992 0 : "yPol1Out.=" << &vecy1resOuter <<
1993 0 : "zPol1Out.=" << &vecz1resOuter <<
1994 0 : "yPol2Out.=" << &vecy2resOuter <<
1995 0 : "zPol2Out.=" << &vecz2resOuter <<
1996 0 : "yInOut.=" << &vecy4res <<
1997 0 : "zInOut.=" << &vecz4res <<
1998 0 : "chi2y1In=" << chi2I1y <<
1999 0 : "chi2z1In=" << chi2I1z <<
2000 0 : "chi2y1InOut="<< chi2IO1y <<
2001 0 : "chi2z1InOut="<< chi2IO1z <<
2002 0 : "chi2y1Out=" << chi2O1y <<
2003 0 : "chi2z1Out=" << chi2O1z <<
2004 0 : "chi2y2In=" << chi2I2y <<
2005 0 : "chi2z2In=" << chi2I2z <<
2006 0 : "chi2y2Out=" << chi2O2y <<
2007 0 : "chi2z2Out=" << chi2O2z <<
2008 0 : "chi2yInOut=" << chi2IOy <<
2009 0 : "chi2zInOut=" << chi2IOz <<
2010 0 : "trletIn.=" << trInner <<
2011 0 : "trletOut.=" << trOuter <<
2012 0 : "nclI=" << nclI <<
2013 0 : "nclO=" << nclO <<
2014 0 : "xinMin=" << xinMin<<
2015 0 : "xinMax=" << xinMax<<
2016 0 : "xoutMin=" << xoutMin<<
2017 0 : "xoutMax=" << xoutMax<<
2018 0 : "msigmaYIn=" <<msigmaYIn<<
2019 0 : "msigmaYOut=" <<msigmaYOut<<
2020 0 : "mqratioIn=" <<mqratioIn<<
2021 0 : "mqratioOut=" << mqratioOut <<
2022 : "\n";
2023 0 : }
2024 0 : }
2025 :
2026 : // wirte residuals information
2027 0 : if (fStreamLevel>5){
2028 0 : TTreeSRedirector *cstream = GetDebugStreamer();
2029 0 : if (cstream){
2030 0 : Float_t dedx = track->GetdEdx();
2031 0 : (*cstream)<<"Residuals"<<
2032 0 : "run="<<fRun<< // run number
2033 0 : "event="<<fEvent<< // event number
2034 0 : "time="<<fTime<< // time stamp of event
2035 0 : "trigger="<<fTrigger<< // trigger
2036 0 : "mag="<<fMagF<< // magnetic field
2037 : //
2038 0 : "cutNr=" << icut <<
2039 0 : "edgeCutX=" << edgeCutX <<
2040 0 : "edgeCutY=" << edgeCutY <<
2041 0 : "nclCut=" << nclCut <<
2042 0 : "LTr.=" << ltrp <<
2043 0 : "Tr.=" << extparam<<
2044 0 : "dEdx=" << dedx <<
2045 0 : "Cl.=" << &arrCl <<
2046 0 : "vX.=" << &vecgX<< // global x
2047 0 : "vY.=" << &vecgY<< // global y
2048 0 : "vZ.=" << &vecgZ<< // global z
2049 0 : "TrX.=" << &vecX <<
2050 0 : "TrYpol1.=" << &vecY1 <<
2051 0 : "TrZpol1.=" << &vecZ1 <<
2052 0 : "TrYpol2.=" << &vecY2 <<
2053 0 : "TrZpol2.=" << &vecZ2 <<
2054 0 : "TrYpol1InOut.="<< &vecY1IO <<
2055 0 : "TrZpol1InOut.="<< &vecZ1IO <<
2056 0 : "TrYInOut.=" << &vecY4 <<
2057 0 : "TrZInOut.=" << &vecZ4 <<
2058 0 : "ClY.=" << &vecClY <<
2059 0 : "ClZ.=" << &vecClZ <<
2060 0 : "isReject.=" << &isReject<<
2061 0 : "sec.=" << &vecSec <<
2062 0 : "nclI=" << nclI <<
2063 0 : "nclO=" << nclO <<
2064 0 : "xinMin=" << xinMin<<
2065 0 : "xinMax=" << xinMax<<
2066 0 : "xoutMin=" << xoutMin<<
2067 0 : "xoutMax=" << xoutMax<<
2068 0 : "msigmaYIn=" <<msigmaYIn<<
2069 0 : "msigmaYOut=" <<msigmaYOut<<
2070 0 : "mqratioIn=" <<mqratioIn<<
2071 0 : "mqratioOut=" << mqratioOut <<
2072 0 : "yInOut.=" << &vecy4res <<
2073 0 : "zInOut.=" << &vecz4res <<
2074 : //chi2s
2075 0 : "chi2y1In=" << chi2I1y << //
2076 0 : "chi2z1In=" << chi2I1z <<
2077 0 : "chi2y1Out=" << chi2O1y <<
2078 0 : "chi2z1Out=" << chi2O1z <<
2079 0 : "chi2y1InOut="<< chi2IO1y <<
2080 0 : "chi2z1InOut="<< chi2IO1z <<
2081 0 : "chi2y2In=" << chi2I2y <<
2082 0 : "chi2z2In=" << chi2I2z <<
2083 0 : "chi2y2Out=" << chi2O2y <<
2084 0 : "chi2z2Out=" << chi2O2z <<
2085 0 : "chi2yInOut=" << chi2IOy <<
2086 0 : "chi2zInOut=" << chi2IOz <<
2087 : // fit parameters
2088 0 : "yPol1In.=" << &vecy1resInner <<
2089 0 : "zPol1In.=" << &vecz1resInner <<
2090 0 : "yPol2In.=" << &vecy2resInner <<
2091 0 : "zPol2In.=" << &vecz2resInner <<
2092 0 : "yPol1Out.=" << &vecy1resOuter <<
2093 0 : "zPol1Out.=" << &vecz1resOuter <<
2094 0 : "yPol1InOut.="<< &vecy1resIO <<
2095 0 : "zPol1InOut.="<< &vecz1resIO <<
2096 0 : "yPol2Out.=" << &vecy2resOuter <<
2097 0 : "zPol2Out.=" << &vecz2resOuter <<
2098 :
2099 : "\n";
2100 :
2101 0 : }
2102 0 : }
2103 : //==========================//
2104 : // Fill Residual Histograms //
2105 : //==========================//
2106 0 : if (!fHisNclIn) UpdateFitHistos();
2107 :
2108 0 : TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2109 0 : TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2110 0 : TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2111 0 : TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2112 0 : TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2113 0 : TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2114 : // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2115 : //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2116 : //
2117 0 : for (Int_t irow=158;irow>-1;--irow) {
2118 0 : if (vecSec[irow]==-1)continue; // no cluster info
2119 0 : if (isReject[irow]>0.5) continue; //
2120 0 : Double_t ycl = vecClY[irow];
2121 0 : Double_t yfit = vecY1[irow];
2122 0 : Double_t yfit2 = vecY2[irow];
2123 0 : Double_t x = vecX[irow];
2124 : Double_t yabsbeam = -1000;
2125 0 : if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2126 0 : yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2127 0 : else if(innerSector==fBeamSectorInner[id])
2128 0 : yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2129 :
2130 : // Double_t yfit3 = vecY2[irow];
2131 0 : Double_t zcl = vecClZ[irow];
2132 0 : Double_t zfit = vecZ1[irow];
2133 0 : Double_t zfit2 = vecZ2[irow];
2134 : //Double_t zfit3 = vecZ2[irow];
2135 :
2136 : // dz abs
2137 : // The expressions for zcorrected has been obtained by
2138 : // inverting the fits in the FitDriftV() method (ignoring the
2139 : // global y dependence for now):
2140 : // A side:
2141 : // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2142 : // =>
2143 : // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2144 : //
2145 : // C side:
2146 : // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
2147 : // =>
2148 : // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2149 :
2150 : Double_t dzabs = -1000;
2151 : Double_t zcorrected = -1000;
2152 0 : if (ltrp->GetSide()==0){
2153 0 : if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2154 : // ignore global y dependence for now
2155 : zcorrected = 0;
2156 0 : if(!fUseFixedDriftV)
2157 0 : zcorrected = (zcl + (*fFitAside)[0] -
2158 0 : (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2159 : else
2160 0 : zcorrected = (zcl + fFixedFitAside0 -
2161 0 : (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2162 : // zcorrected = zcl;
2163 0 : if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2164 0 : dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2165 0 : else if(innerSector==fBeamSectorInner[id])
2166 0 : dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2167 : }
2168 : } else {
2169 0 : if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2170 :
2171 0 : if(!fUseFixedDriftV)
2172 0 : zcorrected = (zcl - (*fFitCside)[0] +
2173 0 : (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2174 : else
2175 0 : zcorrected = (zcl - fFixedFitCside0 +
2176 0 : (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2177 :
2178 : // zcorrected = zcl;
2179 0 : if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2180 0 : dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2181 0 : else if(innerSector==fBeamSectorInner[id])
2182 0 : dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2183 : }
2184 : }
2185 :
2186 0 : if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2187 0 : if (profy){
2188 0 : profy->Fill(irow,ycl-yfit);
2189 0 : profy2->Fill(irow,ycl-yfit2);
2190 0 : if(yabsbeam<-100) {
2191 0 : fHisYAbsErrors->Fill(id);
2192 : // profyabs->Fill(irow,-0.99);
2193 : } else
2194 0 : profyabs->Fill(irow,ycl-yabsbeam);
2195 :
2196 : // profy3->Fill(irow,ycl-yfit3);
2197 : }
2198 0 : if (profz) {
2199 0 : profz->Fill(irow,zcl-zfit);
2200 0 : profz2->Fill(irow,zcl-zfit2);
2201 : //profz3->Fill(irow,zcl-zfit3);
2202 0 : if(dzabs<-100) {
2203 :
2204 0 : fHisZAbsErrors->Fill(id);
2205 : }else
2206 0 : profzabs->Fill(irow,dzabs);
2207 : }
2208 : }
2209 0 : }
2210 : //
2211 : //
2212 : // Fill laser fit histograms
2213 : //
2214 0 : Float_t dedx = track->GetdEdx();
2215 0 : if (nclI>20&&nclO>20){
2216 0 : fHisNclIn->Fill(id,nclI); //->Number of clusters inner
2217 0 : fHisNclOut->Fill(id,nclO); //->Number of clusters outer
2218 0 : fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
2219 : //
2220 0 : fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
2221 0 : fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
2222 0 : fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
2223 : //
2224 0 : fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2225 0 : fHisdZfit->Fill(id,fFitZ[id]);
2226 : //
2227 : //
2228 0 : fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
2229 0 : fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
2230 0 : fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
2231 0 : fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
2232 0 : fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
2233 :
2234 :
2235 0 : fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
2236 0 : fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
2237 0 : fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
2238 0 : fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
2239 0 : fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
2240 : //
2241 : //
2242 0 : fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
2243 0 : fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
2244 0 : fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
2245 : //
2246 0 : fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
2247 0 : fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
2248 0 : fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
2249 : //
2250 0 : fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
2251 0 : fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
2252 0 : fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
2253 0 : fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
2254 : //
2255 0 : fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
2256 0 : fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
2257 0 : fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
2258 0 : fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
2259 : }
2260 0 : if(nclI>20){
2261 0 : fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
2262 0 : fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
2263 : }
2264 : //
2265 0 : if (nclO>20){
2266 0 : fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
2267 0 : fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
2268 : }
2269 :
2270 0 : }
2271 : //
2272 : // Fill raw THnSparses
2273 : //
2274 0 : for (Int_t irow=0;irow<kMaxRow;irow++) {
2275 0 : AliTPCclusterMI *c=track->GetClusterPointer(irow);
2276 0 : if (!c) continue;
2277 0 : if (c->GetMax()>800) continue; // saturation cut
2278 : //if (TMath::Sqrt(TMath::Abs(c->GetSigmaY2()))>1) continue; // saturation cut
2279 : //
2280 0 : Double_t deltaY=c->GetY()-(*ltrp->GetVecLY())[irow];
2281 0 : Double_t deltaZ=c->GetZ()-(*ltrp->GetVecLZ())[irow];
2282 : //TString axisName[6]={"Delta","bin", "rms shape", "Q", "row","trackID"}
2283 0 : Double_t xyz[6]={0, 0, 0,TMath::Sqrt(c->GetMax()), static_cast<Double_t>(irow), static_cast<Double_t>(id)};
2284 0 : xyz[0]=deltaY;
2285 0 : xyz[1]=c->GetPad();
2286 0 : xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaY2()));
2287 0 : fHisLaserPad->Fill(xyz);
2288 0 : xyz[0]=deltaZ;
2289 0 : xyz[1]=c->GetTimeBin();
2290 0 : xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaZ2()));
2291 0 : fHisLaserTime->Fill(xyz);
2292 0 : }
2293 0 : }
2294 :
2295 :
2296 :
2297 : void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2298 : //
2299 : // Dump information about laser beams
2300 : // isOK variable indicates usability of the beam
2301 : // Beam is not usable if:
2302 : // a. No entries in range (krmsCut0)
2303 : // b. Big sperad (krmscut1)
2304 : // c. RMSto fit sigma bigger then (kmultiCut)
2305 : // d. Too big angular spread
2306 : //
2307 :
2308 : const Float_t krmsCut0=0.001;
2309 : const Float_t krmsCut1=0.16;
2310 : const Float_t kmultiCut=2;
2311 : const Float_t kcutP0=0.002;
2312 0 : AliMagF* magF= dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
2313 0 : Double_t xyz[3]={90,0,10}; // tmp. global position
2314 0 : Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical
2315 0 : Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local
2316 : //
2317 : AliTPCcalibLaser *laser = this;
2318 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2319 0 : TF1 fg("fg","gaus");
2320 : AliTPCParam * tpcparam = 0;
2321 : // start set up for absolute residuals analysis
2322 : //
2323 0 : AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
2324 0 : tpcparam = calib->GetParameters();
2325 0 : if (!tpcparam) tpcparam = new AliTPCParamSR;
2326 0 : tpcparam->Update();
2327 0 : AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2328 : Float_t current=0;
2329 0 : Float_t bfield = 0, bz=0;
2330 :
2331 0 : if (grp){
2332 0 : Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
2333 0 : current = grp->GetL3Current((AliGRPObject::Stats)0);
2334 0 : bfield = polarity*5*current/30000.;
2335 : bz = polarity*5*current/30000.;
2336 0 : printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2337 0 : }
2338 :
2339 0 : SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2340 0 : SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2341 0 : TLinearFitter lfabsyInner(2);
2342 0 : lfabsyInner.SetFormula("1 ++ x");
2343 0 : TLinearFitter lfabszInner(2);
2344 0 : lfabszInner.SetFormula("1 ++ x");
2345 :
2346 0 : TLinearFitter lfabsyOuter(2);
2347 0 : lfabsyOuter.SetFormula("1 ++ x");
2348 0 : TLinearFitter lfabszOuter(2);
2349 0 : lfabszOuter.SetFormula("1 ++ x");
2350 : // end set up for absolute residuals analysis
2351 :
2352 : //
2353 : //
2354 0 : for (Int_t id=0; id<336; id++){
2355 0 : Bool_t isOK=kTRUE;
2356 0 : TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
2357 0 : TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2358 0 : TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
2359 0 : TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
2360 0 : TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
2361 0 : TH1F * hisS = (TH1F*)laser->fSignals.At(id);
2362 : //if (!hisphi) continue;
2363 0 : Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2364 : //if (entries<minEntries) continue;
2365 : //
2366 0 : AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2367 0 : if (!ltrp) {
2368 0 : AliTPCLaserTrack::LoadTracks();
2369 0 : ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2370 0 : }
2371 0 : ltrp->UpdatePoints();
2372 0 : pcstream->GetFile()->cd();
2373 0 : if (hisphi) hisphi->Write();
2374 0 : if (hisphiP) hisphiP->Write();
2375 0 : if (hisZ) hisZ->Write();
2376 0 : if (hisP3) hisP3->Write();
2377 0 : if (hisP4) hisP4->Write();
2378 :
2379 0 : Float_t meanphi = hisphi->GetMean();
2380 0 : Float_t rmsphi = hisphi->GetRMS();
2381 : //
2382 0 : Float_t meanphiP = hisphiP->GetMean();
2383 0 : Float_t rmsphiP = hisphiP->GetRMS();
2384 0 : Float_t meanZ = hisZ->GetMean();
2385 0 : Float_t rmsZ = hisZ->GetRMS();
2386 0 : if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2387 0 : hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2388 0 : Double_t gphi1 = fg.GetParameter(1);
2389 0 : Double_t gphi2 = fg.GetParameter(2);
2390 0 : if (hisphiP->GetRMS()>0)
2391 0 : hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2392 0 : Double_t gphiP1 = fg.GetParameter(1);
2393 0 : Double_t gphiP2 = fg.GetParameter(2);
2394 : //
2395 0 : if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2396 0 : hisZ->Fit(&fg,"","");
2397 0 : Double_t gz1 = fg.GetParameter(1);
2398 0 : Double_t gz2 = fg.GetParameter(2);
2399 : //
2400 0 : if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2401 0 : hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2402 0 : Double_t gp31 = fg.GetParameter(1);
2403 0 : Double_t gp32 = fg.GetParameter(2);
2404 0 : Double_t meanp3 = hisP3->GetMean();
2405 0 : Double_t rmsp3 = hisP3->GetRMS();
2406 : //
2407 0 : if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2408 0 : hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2409 0 : Double_t gp41 = fg.GetParameter(1);
2410 0 : Double_t gp42 = fg.GetParameter(2);
2411 0 : Double_t meanp4 = hisP4->GetMean();
2412 0 : Double_t rmsp4 = hisP4->GetRMS();
2413 : //
2414 0 : Float_t meanS=hisS->GetMean();
2415 : //
2416 0 : Double_t lxyz[3];
2417 0 : Double_t lpxyz[3];
2418 0 : ltrp->GetXYZ(lxyz);
2419 0 : ltrp->GetPxPyPz(lpxyz);
2420 :
2421 0 : if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2422 0 : if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2423 0 : if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2424 0 : if (gphiP2>kcutP0) isOK=kFALSE;
2425 : //
2426 : //
2427 : //
2428 : TH1 *his =0;
2429 : //
2430 0 : his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2431 0 : Float_t mnclIn = his->GetMean();
2432 0 : delete his;
2433 0 : his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2434 0 : Float_t mnclOut = his->GetMean();
2435 0 : delete his;
2436 0 : his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2437 0 : Float_t mnclIO = his->GetMean();
2438 0 : delete his;
2439 0 : his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2440 0 : Float_t mLclIn = his->GetMean();
2441 0 : delete his;
2442 0 : his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2443 0 : Float_t mLclOut = his->GetMean();
2444 0 : delete his;
2445 0 : his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2446 0 : Float_t mLclIO = his->GetMean();
2447 0 : delete his;
2448 : //
2449 0 : his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2450 0 : Float_t mdEdx = his->GetMean();
2451 0 : delete his;
2452 : //
2453 : //
2454 : //
2455 : //
2456 0 : his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2457 0 : Float_t mChi2YIn1= his->GetMean();
2458 0 : delete his;
2459 0 : his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2460 0 : Float_t mChi2YOut1 = his->GetMean();
2461 0 : delete his;
2462 0 : his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2463 0 : Float_t mChi2YIn2 = his->GetMean();
2464 0 : delete his;
2465 0 : his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2466 0 : Float_t mChi2YOut2 = his->GetMean();
2467 0 : delete his;
2468 0 : his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2469 0 : Float_t mChi2YIO1 = his->GetMean();
2470 0 : delete his;
2471 0 : his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2472 0 : Float_t mChi2ZIn1 = his->GetMean();
2473 0 : delete his;
2474 0 : his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2475 0 : Float_t mChi2ZOut1 = his->GetMean();
2476 0 : delete his;
2477 0 : his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2478 0 : Float_t mChi2ZIn2 = his->GetMean();
2479 0 : delete his;
2480 0 : his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2481 0 : Float_t mChi2ZOut2 = his->GetMean();
2482 0 : delete his;
2483 0 : his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2484 0 : Float_t mChi2ZIO1 = his->GetMean();
2485 0 : delete his;
2486 : //
2487 : // fit res. histos
2488 : //
2489 0 : his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2490 0 : Float_t edZfit = his->GetEntries();
2491 0 : Float_t mdZfit = his->GetMean();
2492 0 : Float_t rdZfit = his->GetRMS();
2493 0 : delete his;
2494 :
2495 0 : his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2496 0 : Float_t ePy1vP0 = his->GetEntries();
2497 0 : Float_t mPy1vP0 = his->GetMean();
2498 0 : Float_t rPy1vP0 = his->GetRMS();
2499 0 : delete his;
2500 :
2501 0 : his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2502 0 : Float_t ePy2vP0 = his->GetEntries();
2503 0 : Float_t mPy2vP0 = his->GetMean();
2504 0 : Float_t rPy2vP0 = his->GetRMS();
2505 0 : delete his;
2506 :
2507 0 : his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2508 0 : Float_t ePy3vP0 = his->GetEntries();
2509 0 : Float_t mPy3vP0 = his->GetMean();
2510 0 : Float_t rPy3vP0 = his->GetRMS();
2511 0 : delete his;
2512 :
2513 0 : his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2514 0 : Float_t ePy1vP1 = his->GetEntries();
2515 0 : Float_t mPy1vP1 = his->GetMean();
2516 0 : Float_t rPy1vP1 = his->GetRMS();
2517 0 : delete his;
2518 :
2519 0 : his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2520 0 : Float_t ePy2vP1 = his->GetEntries();
2521 0 : Float_t mPy2vP1 = his->GetMean();
2522 0 : Float_t rPy2vP1 = his->GetRMS();
2523 0 : delete his;
2524 :
2525 0 : his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2526 0 : Float_t ePy3vP1 = his->GetEntries();
2527 0 : Float_t mPy3vP1 = his->GetMean();
2528 0 : Float_t rPy3vP1 = his->GetRMS();
2529 0 : delete his;
2530 :
2531 0 : his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2532 0 : Float_t ePy2vP2In = his->GetEntries();
2533 0 : Float_t mPy2vP2In = his->GetMean();
2534 0 : Float_t rPy2vP2In = his->GetRMS();
2535 0 : delete his;
2536 :
2537 0 : his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2538 0 : Float_t ePy2vP2Out = his->GetEntries();
2539 0 : Float_t mPy2vP2Out = his->GetMean();
2540 0 : Float_t rPy2vP2Out = his->GetRMS();
2541 0 : delete his;
2542 :
2543 0 : his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2544 0 : Float_t ePy3vP2IO = his->GetEntries();
2545 0 : Float_t mPy3vP2IO = his->GetMean();
2546 0 : Float_t rPy3vP2IO = his->GetRMS();
2547 0 : delete his;
2548 :
2549 : //
2550 : //
2551 0 : his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2552 0 : Float_t ePz1vP0 = his->GetEntries();
2553 0 : Float_t mPz1vP0 = his->GetMean();
2554 0 : Float_t rPz1vP0 = his->GetRMS();
2555 0 : delete his;
2556 :
2557 0 : his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2558 0 : Float_t ePz2vP0 = his->GetEntries();
2559 0 : Float_t mPz2vP0 = his->GetMean();
2560 0 : Float_t rPz2vP0 = his->GetRMS();
2561 0 : delete his;
2562 :
2563 0 : his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2564 0 : Float_t ePz3vP0 = his->GetEntries();
2565 0 : Float_t mPz3vP0 = his->GetMean();
2566 0 : Float_t rPz3vP0 = his->GetRMS();
2567 0 : delete his;
2568 :
2569 0 : his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2570 0 : Float_t ePz1vP1 = his->GetEntries();
2571 0 : Float_t mPz1vP1 = his->GetMean();
2572 0 : Float_t rPz1vP1 = his->GetRMS();
2573 0 : delete his;
2574 :
2575 0 : his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2576 0 : Float_t ePz2vP1 = his->GetEntries();
2577 0 : Float_t mPz2vP1 = his->GetMean();
2578 0 : Float_t rPz2vP1 = his->GetRMS();
2579 0 : delete his;
2580 :
2581 0 : his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2582 0 : Float_t ePz3vP1 = his->GetEntries();
2583 0 : Float_t mPz3vP1 = his->GetMean();
2584 0 : Float_t rPz3vP1 = his->GetRMS();
2585 0 : delete his;
2586 :
2587 0 : his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2588 0 : Float_t ePz2vP2In = his->GetEntries();
2589 0 : Float_t mPz2vP2In = his->GetMean();
2590 0 : Float_t rPz2vP2In = his->GetRMS();
2591 0 : delete his;
2592 :
2593 0 : his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2594 0 : Float_t ePz2vP2Out = his->GetEntries();
2595 0 : Float_t mPz2vP2Out = his->GetMean();
2596 0 : Float_t rPz2vP2Out = his->GetRMS();
2597 0 : delete his;
2598 :
2599 0 : his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2600 0 : Float_t ePz3vP2IO = his->GetEntries();
2601 0 : Float_t mPz3vP2IO = his->GetMean();
2602 0 : Float_t rPz3vP2IO = his->GetRMS();
2603 0 : delete his;
2604 :
2605 : // Fit absolute laser residuals
2606 0 : TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2607 0 : TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2608 :
2609 0 : Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2610 0 : Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2611 :
2612 0 : TVectorD vecX(kMaxRow); // X
2613 0 : TVectorD vecY(kMaxRow); // Y
2614 0 : TVectorD vecR(kMaxRow); // R
2615 0 : TVectorD vecDY(kMaxRow); // absolute residuals in Y
2616 0 : TVectorD vecDZ(kMaxRow); // absolute residuals in Z
2617 0 : TVectorD vecN(kMaxRow); // number of clusters
2618 0 : TVectorD vecEy(kMaxRow); //error y
2619 0 : TVectorD vecEz(kMaxRow); //error z
2620 0 : TVectorD vecPhi(kMaxRow); // local tangent
2621 0 : TVectorD vecPhiR(kMaxRow); // local tangent
2622 : // magnetic field integrals
2623 0 : TVectorD vecIBR(kMaxRow); // radial
2624 0 : TVectorD vecIBRPhi(kMaxRow); // r-phi
2625 0 : TVectorD vecIBLX(kMaxRow); // local x
2626 0 : TVectorD vecIBLY(kMaxRow); // local y
2627 0 : TVectorD vecIBGX(kMaxRow); // local x
2628 0 : TVectorD vecIBGY(kMaxRow); // local y
2629 0 : TVectorD vecIBZ(kMaxRow); // z
2630 : //
2631 0 : for (Int_t irow=0;irow<kMaxRow;irow++){
2632 0 : vecIBR[irow]=0;
2633 0 : vecIBRPhi[irow]=0;
2634 0 : vecIBLX[irow]=0;
2635 0 : vecIBLY[irow]=0;
2636 0 : vecIBGX[irow]=0;
2637 0 : vecIBGY[irow]=0;
2638 0 : vecIBZ[irow]=0;
2639 0 : Double_t gx =(*(ltrp->fVecGX))[irow];
2640 0 : Double_t gy =(*(ltrp->fVecGY))[irow];
2641 0 : Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
2642 0 : Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2643 0 : Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2644 0 : xyz[2]=(*(ltrp->fVecGZ))[irow];
2645 0 : xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2646 0 : xyz[1]=TMath::ATan2(gy,gx);
2647 0 : Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2648 0 : if (magF){
2649 0 : magF->GetTPCIntCyl(xyz,bxyz);
2650 0 : magF->GetTPCInt(gxyz,bgxyz);
2651 0 : vecIBR[irow]=bxyz[0];
2652 0 : vecIBRPhi[irow]=bxyz[1];
2653 : //
2654 0 : vecIBGX[irow]=bgxyz[0];
2655 0 : vecIBGY[irow]=bgxyz[1];
2656 : //
2657 0 : vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
2658 0 : vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2659 : //
2660 :
2661 0 : vecIBZ[irow]=bxyz[2];
2662 0 : }
2663 0 : }
2664 :
2665 :
2666 0 : lfabsyInner.ClearPoints();
2667 0 : lfabszInner.ClearPoints();
2668 0 : lfabsyOuter.ClearPoints();
2669 0 : lfabszOuter.ClearPoints();
2670 : // dummy fit values
2671 0 : Int_t nClInY = 0;
2672 0 : Float_t yAbsInOffset = -100;
2673 0 : Float_t yAbsInSlope = -100;
2674 0 : Float_t yAbsInDeltay = -100;
2675 0 : Int_t nClInZ = 0;
2676 0 : Float_t zAbsInOffset = -100;
2677 0 : Float_t zAbsInSlope = -100;
2678 0 : Float_t zAbsInDeltaz = -100;
2679 0 : Int_t nClOutY = 0;
2680 0 : Float_t yAbsOutOffset = -100;
2681 0 : Float_t yAbsOutSlope = -100;
2682 0 : Float_t yAbsOutDeltay = -100;
2683 0 : Int_t nClOutZ = 0;
2684 0 : Float_t zAbsOutOffset = -100;
2685 0 : Float_t zAbsOutSlope = -100;
2686 0 : Float_t zAbsOutDeltaz = -100;
2687 :
2688 0 : Float_t lasTanPhiLocIn = -100;
2689 0 : Float_t lasTanPhiLocOut = -100;
2690 :
2691 0 : if(histAbsY && histAbsY->GetEntries()>0) {
2692 :
2693 : Double_t rotAngOut = 10;
2694 : Double_t rotAngIn = 10;
2695 0 : if((secInner%36)!=(secOuter%36))
2696 0 : rotAngIn += 20; // 30 degrees
2697 :
2698 : // Calculate laser mirror X position in local frame
2699 : Double_t laserposOut =
2700 0 : TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2701 : Double_t laserposIn =
2702 0 : TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2703 :
2704 : // Calculate laser tan phi in local frame
2705 0 : lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2706 0 : if(lasTanPhiLocOut<0) {
2707 0 : lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2708 0 : lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2709 0 : } else {
2710 :
2711 0 : lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2712 0 : lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2713 : }
2714 :
2715 0 : lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2716 0 : lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2717 :
2718 0 : TProfile* yprof = histAbsY->ProfileX("yprof");
2719 0 : TProfile* zprof = histAbsZ->ProfileX("zprof");
2720 :
2721 0 : for(Int_t bin = 1; bin<=kMaxRow; bin++) {
2722 :
2723 0 : if(yprof->GetBinEntries(bin)<5&&
2724 0 : zprof->GetBinEntries(bin)<5) {
2725 : continue;
2726 : }
2727 :
2728 : // There is a problem in defining inner and outer sectors for
2729 : // the outer beams (0 and 6) where both sectors are OROCs. To
2730 : // make sure there is no overlap row 94 to 99 are cutted.
2731 0 : if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2732 : continue;
2733 :
2734 0 : Int_t row = (bin-1);
2735 0 : if(row>62)
2736 0 : row -= 63;
2737 :
2738 : Bool_t isOuter = kTRUE;
2739 0 : Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2740 :
2741 0 : if(bin<=62 ||
2742 0 : (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2743 :
2744 : isOuter = kFALSE;
2745 0 : sector = TMath::Nint(fBeamSectorInner[id]);
2746 0 : }
2747 :
2748 :
2749 0 : Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2750 0 : vecN[bin-1] =yprof->GetBinEntries(bin);
2751 0 : vecEy[bin-1]=yprof->GetBinError(bin);
2752 0 : vecEz[bin-1]=zprof->GetBinError(bin);
2753 0 : vecX[bin-1] = x;
2754 0 : vecDY[bin-1] = yprof->GetBinContent(bin);
2755 0 : vecDZ[bin-1] = zprof->GetBinContent(bin);
2756 0 : if (bin>0&&bin<kMaxRow){
2757 : //
2758 : //truncated mean - skip first and the last pad row
2759 : //
2760 0 : Int_t firstBin=TMath::Max(bin-5,0);
2761 0 : Int_t lastBin =TMath::Min(bin+5,158);
2762 0 : histAbsY->GetXaxis()->SetRangeUser(firstBin,lastBin);
2763 0 : histAbsY->GetYaxis()->SetRangeUser(-2,2);
2764 0 : vecEy[bin-1]=histAbsY->GetRMS(2);
2765 0 : vecDY[bin-1]=histAbsY->GetMean(2);
2766 0 : histAbsY->GetXaxis()->SetRangeUser(firstBin+2,lastBin-2);//use+-2 bins
2767 0 : histAbsY->GetYaxis()->SetRangeUser(vecDY[bin-1]-4*vecEy[bin-1],vecDY[bin-1]+4*vecEy[bin-1]);
2768 0 : if (yprof->GetBinEntries(bin-1)>0) vecEy[bin-1]=histAbsY->GetRMS(2)/TMath::Sqrt(yprof->GetBinEntries(bin-1));
2769 0 : vecDY[bin-1]=histAbsY->GetMean(2);
2770 0 : }
2771 :
2772 0 : if(!isOuter) { // inner
2773 0 : vecPhi[bin-1]=lasTanPhiLocIn;
2774 : // Calculate local y from residual and database
2775 0 : Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2776 0 : + vecDY[bin-1];
2777 0 : vecY[bin-1] = y;
2778 0 : Double_t r = TMath::Sqrt(x*x + y*y);
2779 0 : vecR[bin-1] = r;
2780 : // Find angle between laser vector and R vector
2781 : // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2782 0 : Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2783 0 : cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2784 0 : vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2785 0 : if(lasTanPhiLocIn<0)
2786 0 : vecPhiR[bin-1]*=-1; // to have the same sign
2787 :
2788 0 : if(yprof->GetBinEntries(bin)>=10) {
2789 0 : lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2790 0 : TMath::Max(yprof->GetBinError(bin), 0.001));
2791 : }
2792 0 : if(zprof->GetBinEntries(bin)>=10) {
2793 0 : lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2794 0 : TMath::Max(zprof->GetBinError(bin), 0.001));
2795 : }
2796 0 : } else { // outer
2797 0 : vecPhi[bin-1]=lasTanPhiLocOut;
2798 : // Calculate local y from residual and database
2799 0 : Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2800 0 : + vecDY[bin-1];
2801 0 : vecY[bin-1] = y;
2802 0 : Double_t r = TMath::Sqrt(x*x + y*y);
2803 0 : vecR[bin-1] = r;
2804 :
2805 0 : Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2806 0 : cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2807 0 : vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2808 0 : if(lasTanPhiLocOut<0)
2809 0 : vecPhiR[bin-1]*=-1; // to have the same sign
2810 :
2811 0 : if(yprof->GetBinEntries(bin)>=10) {
2812 0 : lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2813 0 : TMath::Max(yprof->GetBinError(bin), 0.001));
2814 : }
2815 0 : if(zprof->GetBinEntries(bin)>=10) {
2816 0 : lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2817 0 : TMath::Max(zprof->GetBinError(bin), 0.001));
2818 : }
2819 : }
2820 : // global position
2821 :
2822 0 : }
2823 :
2824 0 : delete yprof; delete zprof;
2825 :
2826 :
2827 : // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2828 0 : nClInY = lfabsyInner.GetNpoints();
2829 0 : if(lfabsyInner.GetNpoints()>10) {
2830 0 : lfabsyInner.EvalRobust(0.95);
2831 :
2832 0 : TVectorD result(2);
2833 0 : lfabsyInner.GetParameters(result);
2834 0 : yAbsInOffset = result[0];
2835 0 : yAbsInSlope = result[1];
2836 0 : yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2837 0 : }
2838 0 : nClInZ = lfabszInner.GetNpoints();
2839 0 : if(lfabszInner.GetNpoints()>10) {
2840 0 : lfabszInner.EvalRobust(0.95);
2841 :
2842 0 : TVectorD result(2);
2843 0 : lfabszInner.GetParameters(result);
2844 0 : zAbsInOffset = result[0];
2845 0 : zAbsInSlope = result[1];
2846 0 : zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2847 0 : }
2848 0 : nClOutY = lfabsyOuter.GetNpoints();
2849 0 : if(lfabsyOuter.GetNpoints()>10) {
2850 0 : lfabsyOuter.EvalRobust(0.95);
2851 :
2852 0 : TVectorD result(2);
2853 0 : lfabsyOuter.GetParameters(result);
2854 0 : yAbsOutOffset = result[0];
2855 0 : yAbsOutSlope = result[1];
2856 0 : yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2857 0 : }
2858 0 : nClOutZ = lfabszOuter.GetNpoints();
2859 0 : if(lfabszOuter.GetNpoints()>10) {
2860 0 : lfabszOuter.EvalRobust(0.95);
2861 :
2862 0 : TVectorD result(2);
2863 0 : lfabszOuter.GetParameters(result);
2864 0 : zAbsOutOffset = result[0];
2865 0 : zAbsOutSlope = result[1];
2866 0 : zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2867 0 : }
2868 0 : }
2869 :
2870 :
2871 : Int_t itime=-1;
2872 0 : Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2873 0 : Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2874 0 : Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2875 0 : Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2876 : //
2877 0 : Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2878 0 : Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2879 : //
2880 0 : Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2881 0 : Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2882 :
2883 : //
2884 0 : (*pcstream)<<"Mean"<<
2885 0 : "run="<<run<< //
2886 : //voltages
2887 0 : "UcIA="<<coverIA<<
2888 0 : "UcIC="<<coverIC<<
2889 0 : "UcOA="<<coverOA<<
2890 0 : "UcOC="<<coverOC<<
2891 0 : "UsA="<<skirtA<<
2892 0 : "UsC="<<skirtC<<
2893 0 : "UggA="<<ggOffA<<
2894 0 : "UggC="<<ggOffC<<
2895 : //
2896 0 : "isOK="<<isOK<< //
2897 0 : "id="<<id<< // track id
2898 0 : "entries="<<entries<< // number of entries
2899 0 : "bz="<<bfield<< // bfield
2900 0 : "LTr.="<<ltrp<< // refernece track
2901 : //
2902 0 : "mnclIn="<<mnclIn<< // mean number of clusters in inner
2903 0 : "mnclOut="<<mnclOut<< // mean number of clusters in outer
2904 0 : "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2905 0 : "mLclIn="<<mLclIn<< // mean number of clusters in inner
2906 0 : "mLclOut="<<mLclOut<< // mean number of clusters in outer
2907 0 : "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2908 0 : "mdEdx="<<mdEdx<< // mean dEdx
2909 0 : "edZfit="<<edZfit<< // entries z fit
2910 0 : "mdZfit="<<mdZfit<< // mean z fit
2911 0 : "rdZfit="<<rdZfit<< // RMS z fit
2912 : //
2913 : //
2914 0 : "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2915 0 : "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2916 0 : "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2917 0 : "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2918 0 : "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2919 0 : "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2920 0 : "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2921 0 : "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2922 0 : "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2923 0 : "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2924 : //
2925 : //
2926 0 : "ePy1vP0="<<ePy1vP0<<
2927 0 : "mPy1vP0="<<mPy1vP0<<
2928 0 : "rPy1vP0="<<rPy1vP0<<
2929 0 : "ePy2vP0="<<ePy2vP0<<
2930 0 : "mPy2vP0="<<mPy2vP0<<
2931 0 : "rPy2vP0="<<rPy2vP0<<
2932 0 : "ePy3vP0="<<ePy3vP0<<
2933 0 : "mPy3vP0="<<mPy3vP0<<
2934 0 : "rPy3vP0="<<rPy3vP0<<
2935 0 : "ePy1vP1="<<ePy1vP1<<
2936 0 : "mPy1vP1="<<mPy1vP1<<
2937 0 : "rPy1vP1="<<rPy1vP1<<
2938 0 : "ePy2vP1="<<ePy2vP1<<
2939 0 : "mPy2vP1="<<mPy2vP1<<
2940 0 : "rPy2vP1="<<rPy2vP1<<
2941 0 : "ePy3vP1="<<ePy3vP1<<
2942 0 : "mPy3vP1="<<mPy3vP1<<
2943 0 : "rPy3vP1="<<rPy3vP1<<
2944 0 : "ePy2vP2In="<<ePy2vP2In<<
2945 0 : "mPy2vP2In="<<mPy2vP2In<<
2946 0 : "rPy2vP2In="<<rPy2vP2In<<
2947 0 : "ePy2vP2Out="<<ePy2vP2Out<<
2948 0 : "mPy2vP2Out="<<mPy2vP2Out<<
2949 0 : "rPy2vP2Out="<<rPy2vP2Out<<
2950 0 : "ePy3vP2IO="<<ePy3vP2IO<<
2951 0 : "mPy3vP2IO="<<mPy3vP2IO<<
2952 0 : "rPy3vP2IO="<<rPy3vP2IO<<
2953 : //
2954 : //
2955 0 : "ePz1vP0="<<ePz1vP0<<
2956 0 : "mPz1vP0="<<mPz1vP0<<
2957 0 : "rPz1vP0="<<rPz1vP0<<
2958 0 : "ePz2vP0="<<ePz2vP0<<
2959 0 : "mPz2vP0="<<mPz2vP0<<
2960 0 : "rPz2vP0="<<rPz2vP0<<
2961 0 : "ePz3vP0="<<ePz3vP0<<
2962 0 : "mPz3vP0="<<mPz3vP0<<
2963 0 : "rPz3vP0="<<rPz3vP0<<
2964 0 : "ePz1vP1="<<ePz1vP1<<
2965 0 : "mPz1vP1="<<mPz1vP1<<
2966 0 : "rPz1vP1="<<rPz1vP1<<
2967 0 : "ePz2vP1="<<ePz2vP1<<
2968 0 : "mPz2vP1="<<mPz2vP1<<
2969 0 : "rPz2vP1="<<rPz2vP1<<
2970 0 : "ePz3vP1="<<ePz3vP1<<
2971 0 : "mPz3vP1="<<mPz3vP1<<
2972 0 : "rPz3vP1="<<rPz3vP1<<
2973 0 : "ePz2vP2In="<<ePz2vP2In<<
2974 0 : "mPz2vP2In="<<mPz2vP2In<<
2975 0 : "rPz2vP2In="<<rPz2vP2In<<
2976 0 : "ePz2vP2Out="<<ePz2vP2Out<<
2977 0 : "mPz2vP2Out="<<mPz2vP2Out<<
2978 0 : "rPz2vP2Out="<<rPz2vP2Out<<
2979 0 : "ePz3vP2IO="<<ePz3vP2IO<<
2980 0 : "mPz3vP2IO="<<mPz3vP2IO<<
2981 0 : "rPz3vP2IO="<<rPz3vP2IO<<
2982 : //
2983 : //
2984 : //
2985 0 : "lx0="<<lxyz[0]<< // reference x
2986 0 : "lx1="<<lxyz[1]<< // reference y
2987 0 : "lx2="<<lxyz[2]<< // refernece z
2988 0 : "lpx0="<<lpxyz[0]<< // reference x
2989 0 : "lpx1="<<lpxyz[1]<< // reference y
2990 0 : "lpx2="<<lpxyz[2]<< // refernece z
2991 : //
2992 0 : "msig="<<meanS<<
2993 : //
2994 0 : "mphi="<<meanphi<< //
2995 0 : "rmsphi="<<rmsphi<< //
2996 0 : "gphi1="<<gphi1<<
2997 0 : "gphi2="<<gphi2<<
2998 : //
2999 0 : "mphiP="<<meanphiP<< //
3000 0 : "rmsphiP="<<rmsphiP<< //
3001 0 : "gphiP1="<<gphiP1<<
3002 0 : "gphiP2="<<gphiP2<<
3003 : //
3004 0 : "meanZ="<<meanZ<<
3005 0 : "rmsZ="<<rmsZ<<
3006 0 : "gz1="<<gz1<<
3007 0 : "gz2="<<gz2<<
3008 : //
3009 0 : "gp31="<<gp31<< //gaus mean - tgl
3010 0 : "gp32="<<gp32<< //gaus rms -tgl
3011 0 : "meanp3="<<meanp3<< //mean - tgl
3012 0 : "rmsp3="<<rmsp3<< //rms -tgl
3013 0 : "gp41="<<gp41<< //gaus mean - P4
3014 0 : "gp42="<<gp42<< //gaus rms - P4
3015 0 : "meanp4="<<meanp4<< //mean - P4
3016 0 : "rmsp4="<<rmsp4<< //rms - P4
3017 : // Parameters from abs res analysis
3018 0 : "SecIn="<<secInner<< // inner sector
3019 0 : "SecOut="<<secOuter<< // outer sector
3020 0 : "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
3021 0 : "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
3022 0 : "ibr.="<<&vecIBR<< // radial filed integral
3023 0 : "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
3024 0 : "ibr.="<<&vecIBR<< // radial filed integral
3025 0 : "ibz.="<<&vecIBZ<< // z filed integral
3026 : //
3027 0 : "iblx.="<<&vecIBLX<< // local bx integral
3028 0 : "ibly.="<<&vecIBLY<< // local by integral
3029 0 : "ibgx.="<<&vecIBGX<< // global bx integral
3030 0 : "ibgy.="<<&vecIBGY<< // global by integral
3031 : //
3032 0 : "X.="<<&vecX<< // local x
3033 0 : "Y.="<<&vecY<< // local y
3034 0 : "R.="<<&vecR<< // radius
3035 0 : "dY.="<<&vecDY<< // abs y residuals
3036 0 : "dZ.="<<&vecDZ<< // abs z residuals
3037 0 : "eY.="<<&vecEy<< // error of y residuals
3038 0 : "eZ.="<<&vecEz<< // error of z residuals
3039 0 : "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
3040 0 : "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
3041 0 : "nCl.="<<&vecN<< // number of clusters
3042 : //
3043 0 : "nClInY="<<nClInY<< // Number of clusters for inner y
3044 0 : "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
3045 0 : "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
3046 0 : "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
3047 0 : "nClInZ="<<nClInZ<< // Number of clusters for inner z
3048 0 : "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
3049 0 : "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
3050 0 : "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
3051 : //
3052 0 : "nClOutY="<<nClOutY<< // Number of clusters for outer y
3053 0 : "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
3054 0 : "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
3055 0 : "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
3056 0 : "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
3057 0 : "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
3058 0 : "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
3059 0 : "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
3060 : //
3061 : "\n";
3062 0 : }
3063 0 : delete pcstream;
3064 : /*
3065 : Browse the content
3066 : TFile fmean("laserMean.root");
3067 :
3068 :
3069 : */
3070 :
3071 :
3072 0 : }
3073 :
3074 :
3075 :
3076 : void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3077 : //
3078 : //
3079 : //
3080 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3081 0 : TFile * f = pcstream->GetFile();
3082 0 : f->mkdir("dirphi");
3083 0 : f->mkdir("dirphiP");
3084 0 : f->mkdir("dirZ");
3085 0 : TF1 fp("p1","pol1");
3086 : //
3087 : //
3088 0 : char cut[1000];
3089 0 : char grname[1000];
3090 0 : char grnamefull[1000];
3091 :
3092 0 : Double_t mphi[100];
3093 0 : Double_t mphiP[100];
3094 0 : Double_t smphi[100];
3095 0 : Double_t smphiP[100];
3096 0 : Double_t mZ[100];
3097 0 : Double_t smZ[100];
3098 0 : Double_t bz[100];
3099 0 : Double_t sbz[100];
3100 : // fit parameters
3101 0 : Double_t pphi[3];
3102 0 : Double_t pphiP[3];
3103 0 : Double_t pmZ[3];
3104 :
3105 : //
3106 0 : for (Int_t id=0; id<336; id++){
3107 : // id =205;
3108 0 : snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
3109 0 : Int_t entries = chain->Draw("bz",cut,"goff");
3110 0 : if (entries<3) continue;
3111 : AliTPCLaserTrack *ltrp = 0;
3112 0 : if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3113 0 : ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3114 0 : Double_t lxyz[3];
3115 0 : Double_t lpxyz[3];
3116 0 : ltrp->GetXYZ(lxyz);
3117 0 : ltrp->GetPxPyPz(lpxyz);
3118 :
3119 0 : chain->Draw("bz",cut,"goff");
3120 0 : memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3121 0 : chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3122 0 : memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3123 : //
3124 0 : chain->Draw("gphi1",cut,"goff");
3125 0 : memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3126 0 : chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3127 0 : memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3128 : //
3129 0 : chain->Draw("gphiP1",cut,"goff");
3130 0 : memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3131 0 : chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3132 0 : memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3133 : //
3134 0 : chain->Draw("gz1",cut,"goff");
3135 0 : memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3136 0 : chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3137 0 : memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3138 : //
3139 : //
3140 0 : snprintf(grnamefull,1000,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3141 0 : ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3142 : // store data
3143 : // phi
3144 0 : f->cd("dirphi");
3145 0 : Float_t phiB0 =0;
3146 0 : Float_t phiPB0=0;
3147 0 : Float_t zB0=0;
3148 0 : for (Int_t ientry=0; ientry<entries; ientry++){
3149 0 : if (TMath::Abs(bz[ientry])<0.05){
3150 0 : phiB0 = mphi[ientry];
3151 0 : phiPB0 = mphiP[ientry];
3152 0 : zB0 = mZ[ientry];
3153 0 : }
3154 : }
3155 0 : TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3156 0 : grphi->Draw("a*");
3157 0 : grphi->Fit(&fp);
3158 0 : pphi[0] = fp.GetParameter(0); // offset
3159 0 : pphi[1] = fp.GetParameter(1); // slope
3160 0 : pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3161 0 : snprintf(grname,1000,"phi_id%d",id);
3162 0 : grphi->SetName(grname); grphi->SetTitle(grnamefull);
3163 0 : grphi->GetXaxis()->SetTitle("b_{z} (T)");
3164 0 : grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3165 0 : grphi->SetMaximum(1.2);
3166 0 : grphi->SetMinimum(-1.2);
3167 0 : grphi->Draw("a*");
3168 :
3169 0 : grphi->Write();
3170 0 : gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3171 : // phiP
3172 0 : f->cd("dirphiP)");
3173 0 : TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3174 0 : grphiP->Draw("a*");
3175 0 : grphiP->Fit(&fp);
3176 0 : pphiP[0] = fp.GetParameter(0); // offset
3177 0 : pphiP[1] = fp.GetParameter(1); // slope
3178 0 : pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3179 0 : snprintf(grname,1000,"phiP_id%d",id);
3180 0 : grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3181 0 : grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3182 0 : grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3183 0 : grphiP->SetMaximum(pphiP[0]+0.005);
3184 0 : grphiP->SetMinimum(pphiP[0]-0.005);
3185 :
3186 0 : gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3187 0 : grphiP->Write();
3188 : //
3189 : //Z
3190 0 : f->cd("dirZ");
3191 0 : TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3192 0 : grmZ->Draw("a*");
3193 0 : grmZ->Fit(&fp);
3194 0 : pmZ[0] = fp.GetParameter(0); // offset
3195 0 : pmZ[1] = fp.GetParameter(1); // slope
3196 0 : pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3197 0 : snprintf(grname,1000,"mZ_id%d",id);
3198 0 : grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3199 0 : grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3200 0 : grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3201 :
3202 0 : gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3203 0 : grmZ->Write();
3204 : //
3205 : // P4
3206 : //
3207 :
3208 0 : for (Int_t ientry=0; ientry<entries; ientry++){
3209 0 : (*pcstream)<<"Mean"<<
3210 0 : "id="<<id<<
3211 0 : "LTr.="<<ltrp<<
3212 0 : "entries="<<entries<<
3213 0 : "bz="<<bz[ientry]<<
3214 0 : "lx0="<<lxyz[0]<< // reference x
3215 0 : "lx1="<<lxyz[1]<< // reference y
3216 0 : "lx2="<<lxyz[2]<< // refernece z
3217 0 : "lpx0="<<lpxyz[0]<< // reference x
3218 0 : "lpx1="<<lpxyz[1]<< // reference y
3219 0 : "lpx2="<<lpxyz[2]<< // refernece z
3220 : //values
3221 0 : "phiB0="<<phiB0<< // position shift at 0 field
3222 0 : "phiPB0="<<phiPB0<< // angular shift at 0 field
3223 0 : "zB0="<<zB0<< // z shift for 0 field
3224 : //
3225 0 : "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3226 0 : "pphi0="<<pphi[0]<< // offset
3227 0 : "pphi1="<<pphi[1]<< // slope
3228 0 : "pphi2="<<pphi[2]<< // norm chi2
3229 : //
3230 0 : "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3231 0 : "pphiP0="<<pphiP[0]<< // offset
3232 0 : "pphiP1="<<pphiP[1]<< // slope
3233 0 : "pphiP2="<<pphiP[2]<< // norm chi2
3234 : //
3235 0 : "gz1="<<mZ[ientry]<<
3236 0 : "pmZ0="<<pmZ[0]<< // offset
3237 0 : "pmZ1="<<pmZ[1]<< // slope
3238 0 : "pmZ2="<<pmZ[2]<< // norm chi2
3239 : "\n";
3240 : }
3241 0 : }
3242 :
3243 0 : delete pcstream;
3244 :
3245 0 : }
3246 :
3247 :
3248 : void AliTPCcalibLaser::Analyze(){
3249 : //
3250 : //
3251 : //
3252 0 : }
3253 :
3254 :
3255 : Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3256 :
3257 0 : TIterator* iter = li->MakeIterator();
3258 : AliTPCcalibLaser* cal = 0;
3259 : static Int_t counter0=0;
3260 0 : while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3261 0 : if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3262 : // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3263 0 : return -1;
3264 : }
3265 0 : AliDebug(5,Form("Marging number %d\n", counter0));
3266 0 : counter0++;
3267 : //
3268 0 : MergeFitHistos(cal);
3269 : TH1F *h=0x0;
3270 : TH1F *hm=0x0;
3271 : TH2F *h2=0x0;
3272 : TH2F *h2m=0x0;
3273 : // TProfile *hp=0x0;
3274 : //TProfile *hpm=0x0;
3275 :
3276 0 : for (Int_t id=0; id<336; id++){
3277 : // merge fDeltaZ histograms
3278 0 : hm = (TH1F*)cal->fDeltaZ.At(id);
3279 0 : h = (TH1F*)fDeltaZ.At(id);
3280 0 : if (!h &&hm &&hm->GetEntries()>0) {
3281 0 : h=(TH1F*)hm->Clone();
3282 0 : h->SetDirectory(0);
3283 0 : fDeltaZ.AddAt(h,id);
3284 0 : }
3285 0 : if (h && hm) h->Add(hm);
3286 :
3287 : // merge fP3 histograms
3288 0 : hm = (TH1F*)cal->fDeltaP3.At(id);
3289 0 : h = (TH1F*)fDeltaP3.At(id);
3290 0 : if (!h&&hm &&hm->GetEntries()>0) {
3291 0 : h=(TH1F*)hm->Clone();
3292 0 : h->SetDirectory(0);
3293 0 : fDeltaP3.AddAt(h,id);
3294 0 : }
3295 0 : if (h && hm) h->Add(hm);
3296 : // merge fP4 histograms
3297 0 : hm = (TH1F*)cal->fDeltaP4.At(id);
3298 0 : h = (TH1F*)fDeltaP4.At(id);
3299 0 : if (!h &&hm &&hm->GetEntries()>0) {
3300 0 : h=(TH1F*)hm->Clone();
3301 0 : h->SetDirectory(0);
3302 0 : fDeltaP4.AddAt(h,id);
3303 0 : }
3304 0 : if (h&&hm) h->Add(hm);
3305 :
3306 : //
3307 : // merge fDeltaPhi histograms
3308 0 : hm = (TH1F*)cal->fDeltaPhi.At(id);
3309 0 : h = (TH1F*)fDeltaPhi.At(id);
3310 0 : if (!h &&hm &&hm->GetEntries()>0) {
3311 0 : h= (TH1F*)hm->Clone();
3312 0 : h->SetDirectory(0);
3313 0 : fDeltaPhi.AddAt(h,id);
3314 0 : }
3315 0 : if (h&&hm) h->Add(hm);
3316 : // merge fDeltaPhiP histograms
3317 0 : hm = (TH1F*)cal->fDeltaPhiP.At(id);
3318 0 : h = (TH1F*)fDeltaPhiP.At(id);
3319 0 : if (!h&&hm &&hm->GetEntries()>0) {
3320 0 : h=(TH1F*)hm->Clone();
3321 0 : h->SetDirectory(0);
3322 0 : fDeltaPhiP.AddAt(h,id);
3323 0 : }
3324 0 : if (h&&hm) h->Add(hm);
3325 : // merge fSignals histograms
3326 0 : hm = (TH1F*)cal->fSignals.At(id);
3327 0 : h = (TH1F*)fSignals.At(id);
3328 0 : if (!h&&hm &&hm->GetEntries()>0) {
3329 0 : h=(TH1F*)hm->Clone();
3330 0 : h->SetDirectory(0);
3331 0 : fSignals.AddAt(h,id);
3332 0 : }
3333 0 : if (h&&hm) h->Add(hm);
3334 : //
3335 : //
3336 : // merge ProfileY histograms -0
3337 0 : h2m = (TH2F*)cal->fDeltaYres.At(id);
3338 0 : h2 = (TH2F*)fDeltaYres.At(id);
3339 0 : if (h2m&&h2) h2->Add(h2m);
3340 : //
3341 0 : h2m = (TH2F*)cal->fDeltaZres.At(id);
3342 0 : h2 = (TH2F*)fDeltaZres.At(id);
3343 0 : if (h2m&&h2) h2->Add(h2m);
3344 : // merge ProfileY histograms - 2
3345 0 : h2m = (TH2F*)cal->fDeltaYres2.At(id);
3346 0 : h2 = (TH2F*)fDeltaYres2.At(id);
3347 0 : if (h2m&&h2) h2->Add(h2m);
3348 : //
3349 0 : h2m = (TH2F*)cal->fDeltaZres2.At(id);
3350 0 : h2 = (TH2F*)fDeltaZres2.At(id);
3351 0 : if (h2m&&h2) h2->Add(h2m);
3352 :
3353 : // merge ProfileY histograms - abs
3354 0 : h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3355 0 : h2 = (TH2F*)fDeltaYresAbs.At(id);
3356 0 : if (h2m&&h2) h2->Add(h2m);
3357 0 : if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaYresAbs.AddAt(h2,id);}
3358 0 : h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3359 0 : h2 = (TH2F*)fDeltaZresAbs.At(id);
3360 0 : if (h2m&&h2) h2->Add(h2m);
3361 0 : if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaZresAbs.AddAt(h2,id);}
3362 : // merge ProfileY histograms - 3
3363 : //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3364 : //h2 = (TH2F*)fDeltaYres3.At(id);
3365 : //if (h2m) h2->Add(h2m);
3366 : //
3367 : //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3368 : //h2 = (TH2F*)fDeltaZres3.At(id);
3369 : //if (h2m) h->Add(h2m);
3370 : //
3371 : //
3372 : }
3373 : }
3374 0 : return 0;
3375 0 : }
3376 :
3377 : void AliTPCcalibLaser::MakeFitHistos(){
3378 : //
3379 : // Make a fit histograms
3380 : //
3381 : // Number of clusters
3382 : //
3383 : //TH2F *fHisNclIn; //->Number of clusters inner
3384 : //TH2F *fHisNclOut; //->Number of clusters outer
3385 : //TH2F *fHisNclIO; //->Number of cluster inner outer
3386 : // TH2F *fHisdEdx; //->dEdx histo
3387 0 : fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3388 0 : fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3389 0 : fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,kMaxRow,10,kMaxRow);
3390 : //
3391 0 : fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3392 0 : fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3393 0 : fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,kMaxRow,10,kMaxRow);
3394 : //
3395 0 : fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3396 0 : fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3397 :
3398 : //
3399 : // Chi2
3400 : //
3401 : // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3402 : // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3403 : // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3404 : // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3405 : // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3406 0 : fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3407 0 : fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3408 0 : fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3409 0 : fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3410 0 : fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3411 : // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3412 : // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3413 : // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3414 : // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3415 : // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3416 0 : fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3417 0 : fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3418 0 : fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3419 0 : fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3420 0 : fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3421 : //
3422 : // Fit
3423 : //
3424 : //
3425 : // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3426 : // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3427 : // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3428 : // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3429 : // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3430 : // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3431 : // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3432 : // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3433 : // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3434 0 : fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3435 0 : fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3436 0 : fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3437 0 : fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3438 0 : fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3439 0 : fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3440 0 : fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3441 0 : fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3442 0 : fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3443 : //
3444 : //
3445 : // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3446 : // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3447 : // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3448 : // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3449 : // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3450 : // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3451 : // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3452 : // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3453 : // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3454 0 : fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3455 0 : fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3456 0 : fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3457 0 : fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3458 0 : fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3459 0 : fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3460 0 : fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3461 0 : fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3462 0 : fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3463 :
3464 0 : fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3465 0 : fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3466 :
3467 0 : fHisNclIn->SetDirectory(0); //->Number of clusters inner
3468 0 : fHisNclOut->SetDirectory(0); //->Number of clusters outer
3469 0 : fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3470 0 : fHisLclIn->SetDirectory(0); //->Level arm inner
3471 0 : fHisLclOut->SetDirectory(0); //->Level arm outer
3472 0 : fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3473 0 : fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3474 0 : fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3475 : //
3476 : //
3477 0 : fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3478 0 : fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3479 0 : fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3480 0 : fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3481 0 : fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3482 0 : fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3483 0 : fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3484 0 : fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3485 0 : fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3486 0 : fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3487 : //
3488 : //
3489 0 : fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3490 0 : fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3491 0 : fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3492 0 : fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3493 0 : fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3494 0 : fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3495 0 : fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3496 0 : fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3497 0 : fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3498 : //
3499 : //
3500 0 : fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3501 0 : fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3502 0 : fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3503 0 : fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3504 0 : fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3505 0 : fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3506 0 : fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3507 0 : fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3508 0 : fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3509 :
3510 0 : fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3511 0 : fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3512 :
3513 :
3514 :
3515 : //
3516 : //
3517 : //
3518 0 : for (Int_t id=0; id<336;id++){
3519 0 : TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3520 0 : TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3521 : //TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3522 : TH2F *profy2 = 0;
3523 : TH2F *profz2 = 0;//(TH2F*)fDeltaZres2.UncheckedAt(id);
3524 : TH2F *profyabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
3525 : TH2F *profzabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
3526 : // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3527 : //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3528 0 : if (!profy){
3529 0 : profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3530 0 : profy->SetDirectory(0);
3531 0 : fDeltaYres.AddAt(profy,id);
3532 0 : profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3533 0 : profy2->SetDirectory(0);
3534 0 : fDeltaYres2.AddAt(profy2,id);
3535 0 : if(!fUseFixedDriftV)
3536 0 : profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3537 : else
3538 0 : 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
3539 0 : profyabs->SetDirectory(0);
3540 0 : fDeltaYresAbs.AddAt(profyabs,id);
3541 : //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3542 : //profy3->SetDirectory(0);
3543 : //fDeltaYres3.AddAt(profy3,id);
3544 0 : }
3545 0 : if (!profz){
3546 0 : profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3547 0 : profz->SetDirectory(0);
3548 0 : fDeltaZres.AddAt(profz,id);
3549 0 : profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3550 0 : profz2->SetDirectory(0);
3551 0 : fDeltaZres2.AddAt(profz2,id);
3552 0 : if(!fUseFixedDriftV)
3553 0 : profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3554 : else
3555 0 : 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
3556 0 : profzabs->SetDirectory(0);
3557 0 : fDeltaZresAbs.AddAt(profzabs,id);
3558 : //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3559 : //profz3->SetDirectory(0);
3560 : //fDeltaZres3.AddAt(profz3,id);
3561 0 : }
3562 : }
3563 : //
3564 : //
3565 0 : for (Int_t id=0; id<336;id++){
3566 0 : TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3567 : //TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3568 : TH1F * hisP3 = 0;
3569 : TH1F * hisP4 = 0;
3570 :
3571 : TH1F * hisdphi = 0;//(TH1F*)fDeltaPhi.At(id);
3572 : TH1F * hisdphiP = 0;//(TH1F*)fDeltaPhiP.At(id);
3573 : TH1F * hisSignal = 0; //(TH1F*)fSignals.At(id);
3574 :
3575 0 : if (!hisdz){
3576 0 : hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3577 0 : hisdz->SetDirectory(0);
3578 0 : fDeltaZ.AddAt(hisdz,id);
3579 :
3580 0 : hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3581 0 : hisP3->SetDirectory(0);
3582 0 : fDeltaP3.AddAt(hisP3,id);
3583 : //
3584 0 : hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3585 0 : hisP4->SetDirectory(0);
3586 0 : fDeltaP4.AddAt(hisP4,id);
3587 :
3588 : //
3589 0 : hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3590 0 : hisdphi->SetDirectory(0);
3591 0 : fDeltaPhi.AddAt(hisdphi,id);
3592 : //
3593 0 : hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3594 0 : hisdphiP->SetDirectory(0);
3595 0 : fDeltaPhiP.AddAt(hisdphiP,id);
3596 0 : hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3597 0 : hisSignal->SetDirectory(0);
3598 0 : fSignals.AddAt(hisSignal,id);
3599 0 : }
3600 : }
3601 :
3602 : //
3603 : // Make THnSparse
3604 : //
3605 : // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3606 0 : Int_t binsLaser[12]= {336, //id
3607 : 2, //side
3608 : 6, //rod
3609 : 4, //bundle
3610 : 7, //beam
3611 : 300, //dP0
3612 : 300, //dP1
3613 : 300, //dP2
3614 : 300, //dP3
3615 : 300, //dP4
3616 : 80, //ncl
3617 : 50}; //dEdx
3618 0 : Double_t xminLaser[12]= {0, //id
3619 : 0, //side
3620 : 0, //rod
3621 : 0, //bundle
3622 : 0, //beam
3623 : -1, //dP0
3624 : -1, //dP1
3625 : -0.01, //dP2
3626 : -0.01, //dP3
3627 : -0.1, //dP4
3628 : 0, //ncl
3629 : 0}; //sqrt dEdx
3630 0 : Double_t xmaxLaser[12]= {336, //id
3631 : 2, //side
3632 : 6, //rod
3633 : 4, //bundle
3634 : 7, //beam
3635 : 1, //dP0
3636 : 1, //dP1
3637 : 0.01, //dP2
3638 : 0.01, //dP3
3639 : 0.1, //dP4
3640 : 160, //ncl
3641 : 40}; //sqrt dEdx
3642 :
3643 0 : TString nameLaser[12]= {"id",
3644 0 : "side",
3645 0 : "rod",
3646 0 : "bundle",
3647 0 : "beam",
3648 0 : "dP0",
3649 0 : "dP1",
3650 0 : "dP2",
3651 0 : "dP3",
3652 0 : "dP4",
3653 0 : "ncl",
3654 0 : "sqrt dEdx"};
3655 0 : TString titleLaser[12]= {"id",
3656 0 : "side",
3657 0 : "rod",
3658 0 : "bundle",
3659 0 : "beam",
3660 0 : "#Delta_{P0}",
3661 0 : "#Delta_{P1}",
3662 0 : "#Delta_{P2}",
3663 0 : "#Delta_{P3}",
3664 0 : "#Delta_{P4}",
3665 0 : "N_{cl}",
3666 0 : "#sqrt{dEdx}"};
3667 0 : fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3668 0 : for (Int_t iaxis=1; iaxis<12; iaxis++){
3669 0 : fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3670 0 : fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3671 : }
3672 : //
3673 : // Delta Time bin
3674 : // Pad SigmaShape Q charge pad row trackID
3675 0 : Int_t binsRow[6]={200, 10000, 20, 30, kMaxRow, 336};
3676 0 : Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
3677 0 : Double_t axisMax[6]={ 1, 1000, 1, 30, kMaxRow, 336};
3678 0 : TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
3679 :
3680 0 : binsRow[1]=2000;
3681 0 : axisMin[1]=0;
3682 0 : axisMax[1]=200;
3683 0 : fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3684 : //
3685 0 : binsRow[0]=1000;
3686 0 : axisMin[0]=-20;
3687 0 : axisMax[0]=20;
3688 0 : binsRow[1]=10000;
3689 0 : axisMin[1]=0;
3690 0 : axisMax[1]=1000;
3691 : //
3692 0 : fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3693 : //
3694 0 : for (Int_t iaxis=0; iaxis<6; iaxis++){
3695 0 : fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
3696 0 : fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
3697 : }
3698 0 : }
3699 :
3700 : void AliTPCcalibLaser::UpdateFitHistos(){
3701 : //create the fit histos and set the beam parameters(needs OCDB access)
3702 0 : MakeFitHistos();
3703 0 : SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3704 0 : SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3705 0 : SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3706 0 : SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3707 0 : }
3708 :
3709 : void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3710 : //
3711 : // Merge content of histograms
3712 : //
3713 : // Only first histogram is checked - all other should be the same
3714 0 : if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3715 0 : if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
3716 0 : if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
3717 0 : if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
3718 0 : if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
3719 :
3720 0 : if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
3721 0 : if (!fHisNclIn) MakeFitHistos();
3722 0 : if (fHisNclIn->GetEntries()<1) MakeFitHistos();
3723 : //
3724 0 : fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3725 0 : fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3726 0 : fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3727 0 : fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3728 0 : fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3729 0 : fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3730 0 : fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3731 0 : fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3732 : //
3733 : //
3734 0 : fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3735 0 : fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3736 0 : fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3737 0 : fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3738 0 : fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3739 0 : fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3740 0 : fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3741 0 : fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3742 0 : fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3743 0 : fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3744 : //
3745 : //
3746 0 : fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3747 0 : fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3748 0 : fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3749 0 : fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3750 0 : fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3751 0 : fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3752 0 : fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3753 0 : fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3754 0 : fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3755 : //
3756 : //
3757 0 : fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3758 0 : fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3759 0 : fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3760 0 : fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3761 0 : fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3762 0 : fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3763 0 : fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3764 0 : fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3765 0 : fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
3766 0 : fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3767 0 : fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3768 : //
3769 : //
3770 : //
3771 :
3772 :
3773 :
3774 :
3775 0 : }
3776 :
3777 :
3778 :
3779 :
3780 : void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3781 : //
3782 : // Dump fit information - collect information from the streamers
3783 : //
3784 : /*
3785 : TChain * chainFit=0;
3786 : TChain * chainTrack=0;
3787 : TChain * chain=0;
3788 : //
3789 : gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3790 : gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3791 : AliXRDPROOFtoolkit tool;
3792 : chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3793 : chainTrack->Lookup();
3794 : chainTrack->SetProof(kTRUE);
3795 : chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3796 : chainDrift->Lookup();
3797 : chainDrift->SetProof(kTRUE);
3798 :
3799 : chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3800 : chain->Lookup();
3801 : chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3802 : chainFit->Lookup();
3803 : chainFit->SetProof(kTRUE);
3804 : chain->SetProof(kTRUE);
3805 : AliTPCLaserTrack::LoadTracks();
3806 : //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3807 :
3808 : */
3809 : //
3810 : // Fit cuts
3811 : //
3812 0 : TCut cutP3("abs(Tr.fP[3])<0.1");
3813 0 : TCut cutP4("abs(Tr.fP[4])<0.5");
3814 0 : TCut cutPx = cutP3+cutP4;
3815 0 : TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3816 0 : TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3817 0 : TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3818 0 : TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3819 : //
3820 0 : TCut cutdEdx("sqrt(dEdx)>3");
3821 0 : TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3822 0 : TCut cutN("nclO>20&&nclI>20");
3823 0 : TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3824 : //
3825 : // Cluster cuts
3826 : //
3827 0 : TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3828 0 : TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3829 0 : TCut cutClX("abs(Cl[].fX)>10");
3830 0 : TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3831 0 : TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3832 0 : TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3833 0 : TCut cutQ("sqrt(Cl[].fMax)>4");
3834 0 : TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3835 :
3836 :
3837 : TH1F * phisAl = 0;
3838 : TH1F * phisAccept = 0;
3839 : TH1F * phisOut = 0;
3840 : TProfile * pdEdx = 0;
3841 :
3842 : TProfile * pP0 = 0;
3843 : TProfile * pP1 = 0;
3844 : TProfile * pP2 = 0;
3845 : TProfile * pP3 = 0;
3846 : TProfile * pP4 = 0;
3847 : //
3848 : TProfile * pNclI = 0;
3849 : TProfile * pNclO = 0;
3850 : //
3851 : TProfile * pchi2YIn =0;
3852 : TProfile * pchi2ZIn =0;
3853 : TProfile * pchi2YOut =0;
3854 : TProfile * pchi2ZOut =0;
3855 : TProfile * pchi2YInOut =0;
3856 : TProfile * pchi2ZInOut =0;;
3857 : // laser counters
3858 0 : chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3859 0 : phisAl = (TH1F*)gROOT->FindObject("hisAl");
3860 0 : chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3861 0 : phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3862 0 : chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3863 0 : phisOut = (TH1F*)gROOT->FindObject("hisOut");
3864 : //
3865 0 : chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3866 0 : pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3867 : // track param
3868 : //
3869 0 : chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3870 0 : pP0 = (TProfile*)gROOT->FindObject("hP0");
3871 0 : chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3872 0 : pP1 = (TProfile*)gROOT->FindObject("hP1");
3873 0 : chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3874 0 : pP2 = (TProfile*)gROOT->FindObject("hP2");
3875 0 : chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3876 0 : pP3 = (TProfile*)gROOT->FindObject("hP3");
3877 0 : chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3878 0 : pP4 = (TProfile*)gROOT->FindObject("hP4");
3879 : //
3880 0 : chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3881 0 : pNclI = (TProfile*)gROOT->FindObject("hNclI");
3882 0 : chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3883 0 : pNclO = (TProfile*)gROOT->FindObject("hNclO");
3884 : //
3885 : //
3886 0 : chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3887 0 : pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3888 0 : chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3889 0 : pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3890 0 : chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3891 0 : pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3892 0 : chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3893 0 : pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3894 0 : chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3895 0 : pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3896 0 : chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3897 0 : pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3898 : //
3899 : // second derivatives
3900 : //
3901 0 : TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3902 0 : chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3903 0 : TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3904 0 : chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3905 0 : TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3906 0 : chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3907 : //
3908 0 : phisPy2In->FitSlicesY();
3909 0 : TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3910 0 : TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3911 0 : TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3912 : //
3913 0 : phisPy2Out->FitSlicesY();
3914 0 : TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3915 0 : TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3916 0 : TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3917 : //
3918 0 : phisPy2InOut->FitSlicesY();
3919 0 : TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3920 0 : TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3921 0 : TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3922 : //
3923 0 : TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3924 0 : chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3925 0 : TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3926 0 : chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3927 0 : TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3928 0 : chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3929 : //
3930 0 : phisPz2In->FitSlicesY();
3931 0 : TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3932 0 : TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3933 0 : TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3934 : //
3935 0 : phisPz2Out->FitSlicesY();
3936 0 : TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3937 0 : TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3938 0 : TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3939 : //
3940 0 : phisPz2InOut->FitSlicesY();
3941 0 : TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3942 0 : TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3943 0 : TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3944 : //
3945 : //
3946 : //
3947 :
3948 :
3949 : {
3950 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3951 0 : for (Int_t ilaser=0; ilaser<336; ilaser++){
3952 0 : Float_t all=phisAl->GetBinContent(ilaser+1);
3953 0 : Float_t accept=phisAccept->GetBinContent(ilaser+1);
3954 0 : Float_t out=phisOut->GetBinContent(ilaser+1);
3955 0 : Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3956 0 : Float_t mP0 = pP0->GetBinContent(ilaser+1);
3957 0 : Float_t mP1 = pP1->GetBinContent(ilaser+1);
3958 0 : Float_t mP2 = pP2->GetBinContent(ilaser+1);
3959 0 : Float_t mP3 = pP3->GetBinContent(ilaser+1);
3960 0 : Float_t mP4 = pP4->GetBinContent(ilaser+1);
3961 :
3962 :
3963 0 : Float_t nclI = pNclI->GetBinContent(ilaser+1);
3964 0 : Float_t nclO = pNclO->GetBinContent(ilaser+1);
3965 : //
3966 0 : Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3967 0 : Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3968 0 : Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3969 0 : Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3970 0 : Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3971 0 : Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3972 : //
3973 0 : Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3974 0 : Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3975 0 : Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3976 : //
3977 0 : Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3978 0 : Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3979 0 : Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3980 : //
3981 0 : Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3982 0 : Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3983 0 : Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3984 : //
3985 0 : Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3986 0 : Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3987 0 : Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3988 : //
3989 0 : Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3990 0 : Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3991 0 : Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3992 : //
3993 0 : Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3994 0 : Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3995 0 : Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3996 :
3997 0 : AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3998 0 : (*pcstream)<<"Scan"<<
3999 0 : "Run="<<id<<
4000 0 : "LTr.="<<ltrp<<
4001 0 : "all="<<all<<
4002 0 : "accept="<<accept<<
4003 0 : "out="<<out<<
4004 0 : "sdedx="<<sdedx<<
4005 0 : "mP0="<<mP0<<
4006 0 : "mP1="<<mP1<<
4007 0 : "mP2="<<mP2<<
4008 0 : "mP3="<<mP3<<
4009 0 : "mP4="<<mP4<<
4010 0 : "nclI="<<nclI<<
4011 0 : "nclO="<<nclO<<
4012 0 : "chi2YIn="<<chi2YIn<<
4013 0 : "chi2YOut="<<chi2YOut<<
4014 0 : "chi2YInOut="<<chi2YInOut<<
4015 0 : "chi2ZIn="<<chi2ZIn<<
4016 0 : "chi2ZOut="<<chi2ZOut<<
4017 0 : "chi2ZInOut="<<chi2ZInOut<<
4018 : //
4019 0 : "nPy2In="<<entriesPy2In<<
4020 0 : "mPy2In="<<meanPy2In<<
4021 0 : "sPy2In="<<sigmaPy2In<<
4022 : //
4023 0 : "nPy2Out="<<entriesPy2Out<<
4024 0 : "mPy2Out="<<meanPy2Out<<
4025 0 : "sPy2Out="<<sigmaPy2Out<<
4026 : //
4027 0 : "nPy2InOut="<<entriesPy2InOut<<
4028 0 : "mPy2InOut="<<meanPy2InOut<<
4029 0 : "sPy2InOut="<<sigmaPy2InOut<<
4030 : //
4031 0 : "nPz2In="<<entriesPz2In<<
4032 0 : "mPz2In="<<meanPz2In<<
4033 0 : "sPz2In="<<sigmaPz2In<<
4034 : //
4035 0 : "nPz2Out="<<entriesPz2Out<<
4036 0 : "mPz2Out="<<meanPz2Out<<
4037 0 : "sPz2Out="<<sigmaPz2Out<<
4038 : //
4039 0 : "nPz2InOut="<<entriesPz2InOut<<
4040 0 : "mPz2InOut="<<meanPz2InOut<<
4041 0 : "sPz2InOut="<<sigmaPz2InOut<<
4042 : "\n";
4043 0 : }
4044 :
4045 0 : delete pcstream;
4046 : }
4047 0 : }
4048 :
4049 : void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
4050 : TVectorD& meanSlope,
4051 : TVectorD& sectorArray,
4052 : Int_t option)
4053 : {
4054 : // This method should ideally go in AliTPCLaser
4055 : // option == 0 (pads outer - closest to beam)
4056 : // option == 1 (pads inner)
4057 : // option == 2 (time outer)
4058 : // option == 3 (time inner)
4059 : Int_t nFailures = 0;
4060 :
4061 0 : for(Int_t id = 0; id < 336; id++) {
4062 :
4063 0 : if (!AliTPCLaserTrack::GetTracks())
4064 0 : AliTPCLaserTrack::LoadTracks();
4065 : AliTPCLaserTrack *ltrp =
4066 0 : (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4067 :
4068 0 : AliExternalTrackParam trackParam(*ltrp);
4069 :
4070 : Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
4071 0 : if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
4072 0 : deltaangle = 30.0; // inner sector is 1 sector further away from mirror
4073 :
4074 0 : Double_t angle = trackParam.GetAlpha();
4075 0 : if(angle<0)
4076 0 : angle += 2*TMath::Pi();
4077 0 : if(trackParam.GetSnp()>0) // track points to sector "before"
4078 0 : angle -= deltaangle*TMath::DegToRad();
4079 : else // track points to sector "after"
4080 0 : angle += deltaangle*TMath::DegToRad();
4081 :
4082 0 : Bool_t success = trackParam.Rotate(angle);
4083 :
4084 0 : if(!success) {
4085 : // cout << "WARNING: Rotate failed for ID: " << id << endl;
4086 0 : nFailures++;
4087 0 : }
4088 :
4089 0 : angle *= TMath::RadToDeg();
4090 :
4091 0 : Int_t sector = TMath::Nint((angle-10.0)/20.0);
4092 0 : if(sector<0)
4093 0 : sector += 18;
4094 0 : else if(sector>=18)
4095 0 : sector -= 18;
4096 0 : if(ltrp->GetSide()==1) // C side
4097 0 : sector += 18;
4098 0 : if(option==0 || option==2)
4099 0 : sector += 36;
4100 0 : if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
4101 0 : sector += 36;
4102 :
4103 0 : sectorArray[id] = sector;
4104 :
4105 : const Double_t x0 = 0;
4106 :
4107 0 : Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
4108 0 : Double_t slopez = trackParam.GetTgl();
4109 : // One needs a factor sqrt(1+slopey**2) to take into account the
4110 : // longer path length
4111 0 : slopez *= TMath::Sqrt(1.0 + slopey*slopey);
4112 0 : if(fInverseSlopeZ) // wrong sign in database, should be fixed there
4113 0 : slopez *= -1;
4114 : // Double_t offsetz = trackParam.GetZ();
4115 0 : Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4116 0 : Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4117 0 : if(option==2 || option==3) {
4118 0 : meanOffset[id] = offsetz; meanSlope[id] = slopez;
4119 0 : } else {
4120 0 : meanOffset[id] = offsety; meanSlope[id] = slopey;
4121 : }
4122 0 : }
4123 :
4124 0 : if(nFailures>0)
4125 0 : AliWarning(Form("Rotate method failed %d times", nFailures));
4126 0 : }
4127 :
4128 :
4129 :
4130 : void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
4131 : //
4132 : //
4133 : //input="TPCLaserObjects.root"
4134 : //
4135 : // 0. OBJ: TAxis Delta
4136 : // 1. OBJ: TAxis bin
4137 : // 2. OBJ: TAxis rms shape
4138 : // 3. OBJ: TAxis sqrt(Q)
4139 : // 4. OBJ: TAxis row
4140 : // 5. OBJ: TAxis trackID
4141 :
4142 : const Double_t kSigma=4.;
4143 0 : TFile f(finput);
4144 0 : AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
4145 0 : THnSparse * hisPadInput = laserTPC->fHisLaserPad;
4146 0 : THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
4147 0 : TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
4148 0 : TVectorD meanY(kMaxRow), sigmaY(kMaxRow);
4149 0 : TVectorD meanZ(kMaxRow), sigmaZ(kMaxRow);
4150 0 : TVectorD meanPad(kMaxRow), sigmaPad(kMaxRow);
4151 0 : TVectorD meanTime(kMaxRow), sigmaTime(kMaxRow);
4152 0 : TVectorD meanDPad(kMaxRow), sigmaDPad(kMaxRow);
4153 0 : TVectorD meanDTime(kMaxRow), sigmaDTime(kMaxRow);
4154 0 : TVectorD meandEdx(kMaxRow), sigmadEdx(kMaxRow);
4155 0 : TVectorD meanSTime(kMaxRow), sigmaSTime(kMaxRow);
4156 0 : TVectorD meanSPad(kMaxRow), sigmaSPad(kMaxRow);
4157 0 : TVectorD entries(kMaxRow);
4158 : //
4159 0 : Int_t indexes[10]={0,1,2,3,4,5,6};
4160 : TH1 *his=0;
4161 0 : AliTPCLaserTrack::LoadTracks();
4162 : //
4163 0 : for (Int_t id=0; id<336; id++){ // llop over laser beams
4164 0 : printf("id=\t%d\n",id);
4165 : //
4166 0 : AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4167 : //
4168 0 : hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
4169 0 : hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
4170 : //
4171 0 : his=hisTimeInput->Projection(3);
4172 0 : Int_t firstBindEdx=his->FindFirstBinAbove(0);
4173 0 : Int_t lastBindEdx=his->FindLastBinAbove(0);
4174 0 : hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4175 0 : hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4176 0 : delete his;
4177 : //
4178 0 : his=hisTimeInput->Projection(1);
4179 : // Int_t firstBinTime=his->FindFirstBinAbove(0);
4180 : //Int_t lastBinTime=his->FindLastBinAbove(0);
4181 : //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
4182 0 : delete his;
4183 : //
4184 : //
4185 0 : his=hisTimeInput->Projection(2);
4186 : //Int_t firstBinZ=his->FindFirstBinAbove(0);
4187 : //Int_t lastBinZ=his->FindLastBinAbove(0);
4188 : //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
4189 0 : delete his;
4190 : //
4191 0 : his=hisPadInput->Projection(2);
4192 : // Int_t firstBinY=his->FindFirstBinAbove(0);
4193 : //Int_t lastBinY=his->FindLastBinAbove(0);
4194 : //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
4195 0 : delete his;
4196 : //
4197 : //
4198 : //
4199 0 : THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
4200 0 : THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
4201 : //
4202 : //
4203 0 : for (Int_t irow=0; irow<kMaxRow; irow++){
4204 0 : entries[irow]=0;
4205 0 : if ((*(ltrp->GetVecSec()))[irow] <0) continue;
4206 0 : if ((*(ltrp->GetVecLX()))[irow] <80) continue;
4207 :
4208 0 : hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
4209 0 : hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
4210 : //THnSparse *hisPad = hisPad0->Projection(4,indexes);
4211 : //THnSparse *hisTime = hisTime0->Projection(4,indexes);
4212 : THnSparse *hisPad = hisPad0;
4213 : THnSparse *hisTime = hisTime0;
4214 : //
4215 : // Get mean value of QA variables
4216 : //
4217 : // dEdx
4218 0 : his=hisTime->Projection(3);
4219 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4220 0 : meandEdx[irow] =his->GetMean();
4221 0 : sigmadEdx[irow]=his->GetRMS();
4222 : // Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
4223 : //Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
4224 : // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
4225 : //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
4226 0 : delete his;
4227 : //
4228 : // sigma Time
4229 : //
4230 0 : his=hisTime->Projection(2);
4231 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
4232 0 : meanSTime[irow] =his->GetMean();
4233 0 : sigmaSTime[irow]=his->GetRMS();
4234 : //Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4235 : //Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4236 : // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
4237 0 : delete his;
4238 : //
4239 : // sigma Pad
4240 0 : his=hisPad->Projection(2);
4241 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4242 0 : meanSPad[irow] =his->GetMean();
4243 0 : sigmaSPad[irow]=his->GetRMS();
4244 : // Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4245 : //Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
4246 : // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
4247 0 : delete his;
4248 : //
4249 : // apply selection on QA variables
4250 : //
4251 : //
4252 : //
4253 : // Y
4254 0 : his=hisPad->Projection(0);
4255 0 : entries[irow]=his->GetEntries();
4256 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4257 0 : meanY[irow] =his->GetMean();
4258 0 : sigmaY[irow]=his->GetRMS();
4259 0 : delete his;
4260 : // Z
4261 0 : his=hisTime->Projection(0);
4262 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4263 0 : meanZ[irow] =his->GetMean();
4264 0 : sigmaZ[irow]=his->GetRMS();
4265 0 : delete his;
4266 : // Pad
4267 0 : his=hisPad->Projection(1);
4268 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4269 0 : meanPad[irow] =his->GetMean();
4270 0 : meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
4271 0 : sigmaPad[irow]=his->GetRMS();
4272 0 : delete his;
4273 : // Time
4274 0 : his=hisTime->Projection(1);
4275 0 : his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4276 0 : meanTime[irow] = his->GetMean();
4277 0 : meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
4278 0 : sigmaTime[irow]=his->GetRMS();
4279 0 : delete his;
4280 : //
4281 : //delete hisTime;
4282 : //delete hisPad;
4283 0 : }
4284 : //
4285 : //
4286 : //
4287 0 : (*pcstream)<<"laserClusters"<<
4288 0 : "id="<<id<<
4289 0 : "run="<<run<<
4290 0 : "LTr.="<<ltrp<<
4291 : //
4292 0 : "entries.="<<&entries<<
4293 0 : "my.="<<&meanY<< //mean delta y
4294 0 : "rmsy.="<<&sigmaY<< //rms deltay
4295 0 : "mz.="<<&meanZ<< //mean deltaz
4296 0 : "rmsz.="<<&sigmaZ<< //rms z
4297 : //
4298 0 : "mPad.="<<&meanPad<< // mean pad
4299 0 : "mDPad.="<<&meanDPad<< // mead dpad
4300 0 : "rmsPad.="<<&sigmaPad<< // rms pad
4301 0 : "mTime.="<<&meanTime<<
4302 0 : "mDTime.="<<&meanTime<<
4303 0 : "rmsTime.="<<&sigmaTime<<
4304 : //
4305 0 : "mdEdx.="<<&meandEdx<< //mean dedx
4306 0 : "rmsdEdx.="<<&sigmadEdx<< //rms dedx
4307 0 : "mSPad.="<<&meanSPad<< //mean sigma pad
4308 0 : "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
4309 0 : "mSTime.="<<&meanSTime<< //mean sigma time
4310 0 : "rmsSTime.="<<&sigmaSTime<<
4311 : "\n";
4312 : //
4313 0 : delete hisPad0;
4314 0 : delete hisTime0;
4315 : }
4316 0 : delete pcstream;
4317 :
4318 : /*
4319 :
4320 : */
4321 0 : }
4322 :
4323 : void AliTPCcalibLaser::FitLaserClusters(Int_t run){
4324 : //
4325 : //
4326 : //input="TPCLaserObjects.root"
4327 : //Algorithm:
4328 : // 1. Select cluster candidates, remove outlyers
4329 : // edge clusters
4330 : // clusters with atypical spread (e.g due track overlaps)
4331 : // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
4332 : // 2. Fit the tracklets -per sector - in pad and time coordinate frame
4333 : // Remove outlyers
4334 : // Store info distance of track to pad, time center
4335 : // Fit the correction for distance to the center (sin,cos)
4336 : // 3. Do local fit
4337 : const Double_t kEpsilon=0.000001;
4338 : const Int_t kMinClusters=20;
4339 : const Double_t kEdgeCut=3;
4340 : const Double_t kDistCut=1.5; // cut distance to the ideal track
4341 : const Double_t kDistCutFit=0.5;
4342 : const Double_t kDistCutFitPad=0.25;
4343 : const Double_t kDistCutFitTime=0.25;
4344 : const Int_t kSmoothRow=5;
4345 0 : TFile f("hisLasers.root"); // Input file
4346 0 : TTree * treeInput=(TTree*)f.Get("laserClusters");
4347 0 : TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
4348 0 : TVectorD *vecN=0;
4349 0 : TVectorD *vecMY=0;
4350 0 : TVectorD *vecMZ=0;
4351 0 : TVectorD *vecPad=0;
4352 0 : TVectorD *vecTime=0;
4353 0 : TVectorD *vecSY=0;
4354 0 : TVectorD *vecSZ=0;
4355 0 : TVectorD *meandEdx=0;
4356 0 : TVectorD isOK(kMaxRow);
4357 0 : TVectorD fitPad(kMaxRow);
4358 0 : TVectorD fitTime(kMaxRow);
4359 0 : TVectorD fitPadLocal(kMaxRow);
4360 0 : TVectorD fitTimeLocal(kMaxRow);
4361 0 : TVectorD fitDPad(kMaxRow);
4362 0 : TVectorD fitDTime(kMaxRow);
4363 0 : TVectorD fitIPad(kMaxRow);
4364 0 : TVectorD fitITime(kMaxRow);
4365 0 : Double_t chi2PadIROC=0;
4366 0 : Double_t chi2PadOROC=0;
4367 : //
4368 0 : treeInput->SetBranchAddress("my.",&vecMY);
4369 0 : treeInput->SetBranchAddress("mz.",&vecMZ);
4370 0 : treeInput->SetBranchAddress("mPad.",&vecPad);
4371 0 : treeInput->SetBranchAddress("mTime.",&vecTime);
4372 0 : treeInput->SetBranchAddress("rmsy.",&vecSY);
4373 0 : treeInput->SetBranchAddress("rmsz.",&vecSZ);
4374 0 : treeInput->SetBranchAddress("entries.",&vecN);
4375 0 : treeInput->SetBranchAddress("mdEdx.",&meandEdx);
4376 :
4377 0 : AliTPCLaserTrack::LoadTracks();
4378 : //
4379 : //
4380 0 : TVectorD fitPadIROC(3), fitPadOROC(3);
4381 0 : TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
4382 0 : TVectorD fitTimeIROC(3), fitTimeOROC(3);
4383 0 : TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
4384 : //
4385 0 : AliTPCROC * roc = AliTPCROC::Instance();
4386 0 : Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
4387 :
4388 : //
4389 0 : for (Int_t id=0; id<336; id++){
4390 : //
4391 0 : treeInput->GetEntry(id);
4392 0 : AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4393 0 : Int_t medianEntries = TMath::Nint(TMath::Median(kMaxRow,vecN->GetMatrixArray()));
4394 0 : Double_t medianRMSY = TMath::Median(kMaxRow,vecSY->GetMatrixArray());
4395 0 : Double_t rmsRMSY = TMath::RMS(kMaxRow,vecSY->GetMatrixArray());
4396 0 : Double_t medianRMSZ = TMath::Median(kMaxRow,vecSZ->GetMatrixArray());
4397 0 : Double_t rmsRMSZ = TMath::RMS(kMaxRow,vecSZ->GetMatrixArray());
4398 0 : Double_t mdEdx = TMath::Median(kMaxRow,meandEdx->GetMatrixArray());
4399 0 : Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
4400 0 : Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
4401 0 : TLinearFitter fitterY(2,"pol1");
4402 0 : TLinearFitter fitterZ(2,"pol1");
4403 0 : TLinearFitter fitterPad(2,"pol1");
4404 0 : TLinearFitter fitterTime(2,"pol1");
4405 0 : TLinearFitter fitterPadSin(2,"hyp1");
4406 0 : TLinearFitter fitterTimeSin(3,"hyp2");
4407 : //
4408 : //
4409 0 : for (UInt_t irow=0; irow<kMaxRow; irow++){
4410 0 : fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
4411 0 : fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
4412 0 : Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
4413 0 : isOK[irow]=kFALSE;
4414 0 : fitPad[irow]=0;
4415 0 : fitTime[irow]=0;
4416 0 : Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
4417 0 : Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
4418 0 : (*vecPad)[irow]-=npads*0.5;
4419 : //
4420 0 : if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4421 0 : if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4422 : //
4423 0 : if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
4424 0 : if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
4425 0 : if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
4426 0 : if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
4427 0 : if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
4428 0 : if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
4429 0 : if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
4430 0 : Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
4431 0 : if (TMath::Abs(dEdge)<kEdgeCut) continue;
4432 0 : if (irow<roc->GetNRows(0)){
4433 0 : if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
4434 : }
4435 0 : if (irow>roc->GetNRows(0)){
4436 0 : if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
4437 : }
4438 :
4439 0 : isOK[irow]=kTRUE;
4440 0 : }
4441 : //
4442 : //fit OROC - get delta pad and delta time
4443 : //
4444 0 : fitterPad.ClearPoints();
4445 0 : fitterTime.ClearPoints();
4446 0 : fitterPadSin.ClearPoints();
4447 0 : fitterTimeSin.ClearPoints();
4448 0 : {for (Int_t irow=2; irow<157; irow++){
4449 0 : if (isOK[irow]<0.5) continue;
4450 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4451 0 : if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4452 0 : Double_t y=(*vecPad)[irow];
4453 0 : Double_t z=(*vecTime)[irow];
4454 0 : Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4455 0 : fitterPad.AddPoint(&x,y);
4456 0 : fitterTime.AddPoint(&x,z);
4457 0 : }}
4458 0 : chi2PadOROC=0;
4459 0 : if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4460 0 : fitterPad.Eval();
4461 0 : fitterTime.Eval();
4462 0 : chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4463 0 : for (Int_t irow=2; irow<157; irow++){
4464 0 : if (isOK[irow]<0.5) continue;
4465 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4466 0 : if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4467 0 : Double_t y=(*vecPad)[irow];
4468 0 : Double_t z=(*vecTime)[irow];
4469 0 : Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4470 0 : Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4471 0 : Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4472 0 : fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4473 0 : fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4474 0 : fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4475 0 : fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4476 0 : fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4477 0 : fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4478 0 : if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4479 0 : if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4480 0 : if (isOK[irow]>0){
4481 0 : Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
4482 0 : Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4483 0 : TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4484 0 : fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4485 0 : fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4486 0 : }
4487 0 : }
4488 0 : fitterPadSin.Eval();
4489 0 : fitterTimeSin.Eval();
4490 0 : fitterPadSin.FixParameter(0,0);
4491 0 : fitterTimeSin.FixParameter(0,0);
4492 0 : fitterPadSin.Eval();
4493 0 : fitterTimeSin.Eval();
4494 : //
4495 0 : fitterPad.GetParameters(fitPadOROC);
4496 0 : fitterTime.GetParameters(fitTimeOROC);
4497 0 : fitterPadSin.GetParameters(fitPadOROCSin);
4498 0 : fitterTimeSin.GetParameters(fitTimeOROCSin);
4499 : }
4500 : //
4501 : //
4502 : //fit IROC
4503 : //
4504 0 : fitterPad.ClearPoints();
4505 0 : fitterTime.ClearPoints();
4506 0 : fitterPadSin.ClearPoints();
4507 0 : fitterTimeSin.ClearPoints();
4508 0 : for (Int_t irow=2; irow<157; irow++){
4509 0 : if (isOK[irow]<0.5) continue;
4510 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4511 0 : if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4512 0 : Double_t y=(*vecPad)[irow];
4513 0 : Double_t z=(*vecTime)[irow];
4514 0 : Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4515 0 : fitterPad.AddPoint(&x,y);
4516 0 : fitterTime.AddPoint(&x,z);
4517 0 : }
4518 0 : chi2PadIROC=0;
4519 0 : if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4520 0 : fitterPad.Eval();
4521 0 : fitterTime.Eval();
4522 0 : chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4523 0 : for (Int_t irow=2; irow<157; irow++){
4524 0 : if (isOK[irow]<0.5) continue;
4525 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4526 0 : if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4527 0 : Double_t y=(*vecPad)[irow];
4528 0 : Double_t z=(*vecTime)[irow];
4529 0 : Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4530 0 : Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4531 0 : Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4532 0 : fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4533 0 : fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4534 0 : fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4535 0 : fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4536 0 : fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4537 0 : fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4538 0 : if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4539 0 : if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4540 0 : if (isOK[irow]>0.5){
4541 0 : Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
4542 0 : TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
4543 0 : Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4544 0 : TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4545 0 : fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4546 0 : fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4547 0 : }
4548 0 : }
4549 0 : fitterPadSin.Eval();
4550 0 : fitterTimeSin.Eval();
4551 0 : fitterPadSin.FixParameter(0,0);
4552 0 : fitterTimeSin.FixParameter(0,0);
4553 0 : fitterPadSin.Eval();
4554 0 : fitterTimeSin.Eval();
4555 0 : fitterPad.GetParameters(fitPadIROC);
4556 0 : fitterTime.GetParameters(fitTimeIROC);
4557 0 : fitterPadSin.GetParameters(fitPadIROCSin);
4558 0 : fitterTimeSin.GetParameters(fitTimeIROCSin);
4559 : }
4560 0 : for (Int_t irow=0; irow<kMaxRow; irow++){
4561 0 : if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
4562 0 : if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
4563 0 : if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
4564 0 : if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
4565 : }
4566 0 : for (Int_t irow=kSmoothRow/2; irow<kMaxRow-kSmoothRow/2; irow++){
4567 0 : fitPadLocal[irow]=0;
4568 0 : fitTimeLocal[irow]=0;
4569 0 : if (isOK[irow]<0.5) continue;
4570 0 : Int_t sector=(irow<Int_t(roc->GetNRows(0)))? sectorInner:sectorOuter;
4571 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
4572 : //
4573 0 : TLinearFitter fitterPadLocal(2,"pol1");
4574 0 : TLinearFitter fitterTimeLocal(2,"pol1");
4575 0 : Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
4576 0 : for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
4577 0 : Int_t jrow=irow+delta;
4578 0 : if (jrow<0) jrow=0;
4579 0 : if (jrow>=kMaxRow) jrow=kMaxRow-1;
4580 0 : if (isOK[jrow]<0.5) continue;
4581 0 : if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
4582 0 : Double_t y=(*vecPad)[jrow];
4583 0 : Double_t z=(*vecTime)[jrow];
4584 0 : Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
4585 0 : fitterPadLocal.AddPoint(&x,y);
4586 0 : fitterTimeLocal.AddPoint(&x,z);
4587 0 : }
4588 0 : if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
4589 0 : fitterPadLocal.Eval();
4590 0 : fitterTimeLocal.Eval();
4591 0 : fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
4592 0 : fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
4593 0 : }
4594 : //
4595 : //
4596 0 : (*pcstream)<<"fit"<<
4597 0 : "run="<<run<<
4598 0 : "id="<<id<<
4599 0 : "chi2PadIROC="<<chi2PadIROC<<
4600 0 : "chi2PadOROC="<<chi2PadOROC<<
4601 0 : "mdEdx="<<mdEdx<<
4602 0 : "LTr.="<<ltrp<<
4603 0 : "isOK.="<<&isOK<<
4604 : // mean measured-ideal values
4605 0 : "mY.="<<vecMY<<
4606 0 : "mZ.="<<vecMZ<<
4607 : // local coordinate fit
4608 0 : "mPad.="<<vecPad<<
4609 0 : "mTime.="<<vecTime<<
4610 0 : "fitPad.="<<&fitPad<<
4611 0 : "fitTime.="<<&fitTime<<
4612 0 : "fitPadLocal.="<<&fitPadLocal<<
4613 0 : "fitTimeLocal.="<<&fitTimeLocal<<
4614 0 : "fitDPad.="<<&fitDPad<<
4615 0 : "fitDTime.="<<&fitDTime<<
4616 0 : "fitIPad.="<<&fitIPad<<
4617 0 : "fitITime.="<<&fitITime<<
4618 : //
4619 0 : "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
4620 0 : "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
4621 0 : "fitPadOROC.="<<&fitPadOROC<<
4622 0 : "fitPadOROCSin.="<<&fitPadOROCSin<<
4623 : //
4624 0 : "fitTimeIROC.="<<&fitTimeIROC<<
4625 0 : "fitTimeIROCSin.="<<&fitTimeIROCSin<<
4626 0 : "fitTimeOROC.="<<&fitTimeOROC<<
4627 0 : "fitTimeOROCSin.="<<&fitTimeOROCSin<<
4628 : "\n";
4629 0 : }
4630 0 : delete pcstream;
4631 0 : }
|
