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