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