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