]> 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
correcting AliHLTGlobalTriggerDecision::Clear(): using Delete instead of Clear for...
[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);
0b736a46 2344 Double_t meanp3 = hisP3->GetMean();
2345 Double_t rmsp3 = hisP3->GetRMS();
e5c59218 2346 //
d3ce44cb 2347 if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
592a0c8f 2348 hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
e5c59218 2349 Double_t gp41 = fg.GetParameter(1);
2350 Double_t gp42 = fg.GetParameter(2);
0b736a46 2351 Double_t meanp4 = hisP4->GetMean();
2352 Double_t rmsp4 = hisP4->GetRMS();
e5c59218 2353 //
1fd56785 2354 Float_t meanS=hisS->GetMean();
2b35e8e6 2355 //
2356 Double_t lxyz[3];
2357 Double_t lpxyz[3];
2358 ltrp->GetXYZ(lxyz);
2359 ltrp->GetPxPyPz(lpxyz);
1fd56785 2360
2361 if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2362 if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2363 if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
2364 if (gphiP2>kcutP0) isOK=kFALSE;
be113f6e 2365 //
2366 //
2367 //
2368 TH1 *his =0;
2369 //
2370 his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2371 Float_t mnclIn = his->GetMean();
2372 delete his;
2373 his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2374 Float_t mnclOut = his->GetMean();
2375 delete his;
2376 his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2377 Float_t mnclIO = his->GetMean();
2378 delete his;
2379 his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2380 Float_t mLclIn = his->GetMean();
2381 delete his;
2382 his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2383 Float_t mLclOut = his->GetMean();
2384 delete his;
2385 his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2386 Float_t mLclIO = his->GetMean();
2387 delete his;
2388 //
2389 his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2390 Float_t mdEdx = his->GetMean();
2391 delete his;
2392 //
be113f6e 2393 //
2394 //
2395 //
2396 his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2397 Float_t mChi2YIn1= his->GetMean();
2398 delete his;
2399 his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
2400 Float_t mChi2YOut1 = his->GetMean();
2401 delete his;
2402 his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2403 Float_t mChi2YIn2 = his->GetMean();
2404 delete his;
2405 his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
2406 Float_t mChi2YOut2 = his->GetMean();
2407 delete his;
2408 his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
2409 Float_t mChi2YIO1 = his->GetMean();
2410 delete his;
2411 his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2412 Float_t mChi2ZIn1 = his->GetMean();
2413 delete his;
2414 his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
2415 Float_t mChi2ZOut1 = his->GetMean();
2416 delete his;
2417 his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2418 Float_t mChi2ZIn2 = his->GetMean();
2419 delete his;
2420 his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
2421 Float_t mChi2ZOut2 = his->GetMean();
2422 delete his;
2423 his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
2424 Float_t mChi2ZIO1 = his->GetMean();
2425 delete his;
2426 //
2427 // fit res. histos
2428 //
108953e9 2429 his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2430 Float_t edZfit = his->GetEntries();
2431 Float_t mdZfit = his->GetMean();
2432 Float_t rdZfit = his->GetRMS();
2433 delete his;
2434
be113f6e 2435 his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
2436 Float_t ePy1vP0 = his->GetEntries();
2437 Float_t mPy1vP0 = his->GetMean();
2438 Float_t rPy1vP0 = his->GetRMS();
2439 delete his;
2440
2441 his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
2442 Float_t ePy2vP0 = his->GetEntries();
2443 Float_t mPy2vP0 = his->GetMean();
2444 Float_t rPy2vP0 = his->GetRMS();
2445 delete his;
2446
2447 his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
2448 Float_t ePy3vP0 = his->GetEntries();
2449 Float_t mPy3vP0 = his->GetMean();
2450 Float_t rPy3vP0 = his->GetRMS();
2451 delete his;
2452
2453 his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
2454 Float_t ePy1vP1 = his->GetEntries();
2455 Float_t mPy1vP1 = his->GetMean();
2456 Float_t rPy1vP1 = his->GetRMS();
2457 delete his;
2458
2459 his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
2460 Float_t ePy2vP1 = his->GetEntries();
2461 Float_t mPy2vP1 = his->GetMean();
2462 Float_t rPy2vP1 = his->GetRMS();
2463 delete his;
2464
2465 his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
2466 Float_t ePy3vP1 = his->GetEntries();
2467 Float_t mPy3vP1 = his->GetMean();
2468 Float_t rPy3vP1 = his->GetRMS();
2469 delete his;
2470
2471 his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2472 Float_t ePy2vP2In = his->GetEntries();
2473 Float_t mPy2vP2In = his->GetMean();
2474 Float_t rPy2vP2In = his->GetRMS();
2475 delete his;
2476
2477 his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2478 Float_t ePy2vP2Out = his->GetEntries();
2479 Float_t mPy2vP2Out = his->GetMean();
2480 Float_t rPy2vP2Out = his->GetRMS();
2481 delete his;
2482
2483 his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2484 Float_t ePy3vP2IO = his->GetEntries();
2485 Float_t mPy3vP2IO = his->GetMean();
2486 Float_t rPy3vP2IO = his->GetRMS();
2487 delete his;
2488
2489 //
2490 //
2491 his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
2492 Float_t ePz1vP0 = his->GetEntries();
2493 Float_t mPz1vP0 = his->GetMean();
2494 Float_t rPz1vP0 = his->GetRMS();
2495 delete his;
2496
2497 his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
2498 Float_t ePz2vP0 = his->GetEntries();
2499 Float_t mPz2vP0 = his->GetMean();
2500 Float_t rPz2vP0 = his->GetRMS();
2501 delete his;
2502
2503 his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
2504 Float_t ePz3vP0 = his->GetEntries();
2505 Float_t mPz3vP0 = his->GetMean();
2506 Float_t rPz3vP0 = his->GetRMS();
2507 delete his;
2508
2509 his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
2510 Float_t ePz1vP1 = his->GetEntries();
2511 Float_t mPz1vP1 = his->GetMean();
2512 Float_t rPz1vP1 = his->GetRMS();
2513 delete his;
2514
2515 his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
2516 Float_t ePz2vP1 = his->GetEntries();
2517 Float_t mPz2vP1 = his->GetMean();
2518 Float_t rPz2vP1 = his->GetRMS();
2519 delete his;
2520
2521 his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
2522 Float_t ePz3vP1 = his->GetEntries();
2523 Float_t mPz3vP1 = his->GetMean();
2524 Float_t rPz3vP1 = his->GetRMS();
2525 delete his;
2526
2527 his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
2528 Float_t ePz2vP2In = his->GetEntries();
2529 Float_t mPz2vP2In = his->GetMean();
2530 Float_t rPz2vP2In = his->GetRMS();
2531 delete his;
2532
2533 his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
2534 Float_t ePz2vP2Out = his->GetEntries();
2535 Float_t mPz2vP2Out = his->GetMean();
2536 Float_t rPz2vP2Out = his->GetRMS();
2537 delete his;
2538
2539 his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
2540 Float_t ePz3vP2IO = his->GetEntries();
2541 Float_t mPz3vP2IO = his->GetMean();
2542 Float_t rPz3vP2IO = his->GetRMS();
2543 delete his;
f6220766 2544
2545 // Fit absolute laser residuals
2546 TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2547 TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2548
2549 Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2550 Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2551
2552 TVectorD vecX(159); // X
2553 TVectorD vecY(159); // Y
2554 TVectorD vecR(159); // R
2555 TVectorD vecDY(159); // absolute residuals in Y
2556 TVectorD vecDZ(159); // absolute residuals in Z
2557 TVectorD vecN(159); // number of clusters
2558 TVectorD vecEy(159); //error y
2559 TVectorD vecEz(159); //error z
2560 TVectorD vecPhi(159); // local tangent
2561 TVectorD vecPhiR(159); // local tangent
fb962ae3 2562 // magnetic field integrals
2563 TVectorD vecIBR(159); // radial
2564 TVectorD vecIBRPhi(159); // r-phi
b1c27e4f 2565 TVectorD vecIBLX(159); // local x
2566 TVectorD vecIBLY(159); // local y
2567 TVectorD vecIBGX(159); // local x
2568 TVectorD vecIBGY(159); // local y
fb962ae3 2569 TVectorD vecIBZ(159); // z
2570 //
2571 for (Int_t irow=0;irow<159;irow++){
2572 vecIBR[irow]=0;
2573 vecIBRPhi[irow]=0;
b1c27e4f 2574 vecIBLX[irow]=0;
2575 vecIBLY[irow]=0;
2576 vecIBGX[irow]=0;
2577 vecIBGY[irow]=0;
fb962ae3 2578 vecIBZ[irow]=0;
b1c27e4f 2579 Double_t gx =(*(ltrp->fVecGX))[irow];
2580 Double_t gy =(*(ltrp->fVecGY))[irow];
2581 Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
2582 Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2583 Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
fb962ae3 2584 xyz[2]=(*(ltrp->fVecGZ))[irow];
2585 xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2586 xyz[1]=TMath::ATan2(gy,gx);
b1c27e4f 2587 Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
fb962ae3 2588 if (magF){
2589 magF->GetTPCIntCyl(xyz,bxyz);
b1c27e4f 2590 magF->GetTPCInt(gxyz,bgxyz);
fb962ae3 2591 vecIBR[irow]=bxyz[0];
2592 vecIBRPhi[irow]=bxyz[1];
b1c27e4f 2593 //
2594 vecIBGX[irow]=bgxyz[0];
2595 vecIBGY[irow]=bgxyz[1];
2596 //
2597 vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
2598 vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2599 //
2600
fb962ae3 2601 vecIBZ[irow]=bxyz[2];
2602 }
2603 }
2604
f6220766 2605
2606 lfabsyInner.ClearPoints();
2607 lfabszInner.ClearPoints();
2608 lfabsyOuter.ClearPoints();
2609 lfabszOuter.ClearPoints();
2610 // dummy fit values
2611 Int_t nClInY = 0;
2612 Float_t yAbsInOffset = -100;
2613 Float_t yAbsInSlope = -100;
2614 Float_t yAbsInDeltay = -100;
2615 Int_t nClInZ = 0;
2616 Float_t zAbsInOffset = -100;
2617 Float_t zAbsInSlope = -100;
2618 Float_t zAbsInDeltaz = -100;
2619 Int_t nClOutY = 0;
2620 Float_t yAbsOutOffset = -100;
2621 Float_t yAbsOutSlope = -100;
2622 Float_t yAbsOutDeltay = -100;
2623 Int_t nClOutZ = 0;
2624 Float_t zAbsOutOffset = -100;
2625 Float_t zAbsOutSlope = -100;
2626 Float_t zAbsOutDeltaz = -100;
2627
2628 Float_t lasTanPhiLocIn = -100;
2629 Float_t lasTanPhiLocOut = -100;
2630
2631 if(histAbsY->GetEntries()>0) {
2632
2633 Double_t rotAngOut = 10;
2634 Double_t rotAngIn = 10;
2635 if((secInner%36)!=(secOuter%36))
2636 rotAngIn += 20; // 30 degrees
2637
2638 // Calculate laser mirror X position in local frame
2639 Double_t laserposOut =
2640 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2641 Double_t laserposIn =
2642 TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2643
2644 // Calculate laser tan phi in local frame
2645 lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
2646 if(lasTanPhiLocOut<0) {
2647 lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2648 lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2649 } else {
2650
2651 lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2652 lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2653 }
2654
2655 lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2656 lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2657
2658 TProfile* yprof = histAbsY->ProfileX("yprof");
2659 TProfile* zprof = histAbsZ->ProfileX("zprof");
2660
2661 for(Int_t bin = 1; bin<=159; bin++) {
2662
2663 if(yprof->GetBinEntries(bin)<10&&
2664 zprof->GetBinEntries(bin)<10) {
2665 continue;
2666 }
2667
2668 // There is a problem in defining inner and outer sectors for
2669 // the outer beams (0 and 6) where both sectors are OROCs. To
2670 // make sure there is no overlap row 94 to 99 are cutted.
2671 if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2672 continue;
2673
2674 Int_t row = (bin-1);
2675 if(row>62)
2676 row -= 63;
2677
2678 Bool_t isOuter = kTRUE;
2679 Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2680
2681 if(bin<=62 ||
2682 (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2683
2684 isOuter = kFALSE;
2685 sector = TMath::Nint(fBeamSectorInner[id]);
2686 }
2687
2688
2689 Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2690 vecN[bin-1] =yprof->GetBinEntries(bin);
2691 vecEy[bin-1]=yprof->GetBinError(bin);
2692 vecEz[bin-1]=zprof->GetBinError(bin);
2693 vecX[bin-1] = x;
2694 vecDY[bin-1] = yprof->GetBinContent(bin);
2695 vecDZ[bin-1] = zprof->GetBinContent(bin);
2696 if(!isOuter) { // inner
2697 vecPhi[bin-1]=lasTanPhiLocIn;
2698 // Calculate local y from residual and database
2699 Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2700 + vecDY[bin-1];
2701 vecY[bin-1] = y;
2702 Double_t r = TMath::Sqrt(x*x + y*y);
2703 vecR[bin-1] = r;
2704 // Find angle between laser vector and R vector
2705 // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2706 Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2707 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2708 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2709 if(lasTanPhiLocIn<0)
2710 vecPhiR[bin-1]*=-1; // to have the same sign
2711
2712 if(yprof->GetBinEntries(bin)>=10) {
2713 lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2714 TMath::Max(yprof->GetBinError(bin), 0.001));
2715 }
2716 if(zprof->GetBinEntries(bin)>=10) {
2717 lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2718 TMath::Max(zprof->GetBinError(bin), 0.001));
2719 }
2720 } else { // outer
2721 vecPhi[bin-1]=lasTanPhiLocOut;
2722 // Calculate local y from residual and database
2723 Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2724 + vecDY[bin-1];
2725 vecY[bin-1] = y;
2726 Double_t r = TMath::Sqrt(x*x + y*y);
2727 vecR[bin-1] = r;
2728
2729 Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2730 cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2731 vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2732 if(lasTanPhiLocOut<0)
2733 vecPhiR[bin-1]*=-1; // to have the same sign
2734
2735 if(yprof->GetBinEntries(bin)>=10) {
2736 lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2737 TMath::Max(yprof->GetBinError(bin), 0.001));
2738 }
2739 if(zprof->GetBinEntries(bin)>=10) {
2740 lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2741 TMath::Max(zprof->GetBinError(bin), 0.001));
2742 }
2743 }
fb962ae3 2744 // global position
2745
f6220766 2746 }
2747
2748 delete yprof; delete zprof;
2749
2750
2751 // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2752 nClInY = lfabsyInner.GetNpoints();
2753 if(lfabsyInner.GetNpoints()>10) {
2754 lfabsyInner.EvalRobust(0.95);
2755
2756 TVectorD result(2);
2757 lfabsyInner.GetParameters(result);
2758 yAbsInOffset = result[0];
2759 yAbsInSlope = result[1];
2760 yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2761 }
2762 nClInZ = lfabszInner.GetNpoints();
2763 if(lfabszInner.GetNpoints()>10) {
2764 lfabszInner.EvalRobust(0.95);
2765
2766 TVectorD result(2);
2767 lfabszInner.GetParameters(result);
2768 zAbsInOffset = result[0];
2769 zAbsInSlope = result[1];
2770 zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2771 }
2772 nClOutY = lfabsyOuter.GetNpoints();
2773 if(lfabsyOuter.GetNpoints()>10) {
2774 lfabsyOuter.EvalRobust(0.95);
2775
2776 TVectorD result(2);
2777 lfabsyOuter.GetParameters(result);
2778 yAbsOutOffset = result[0];
2779 yAbsOutSlope = result[1];
2780 yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2781 }
2782 nClOutZ = lfabszOuter.GetNpoints();
2783 if(lfabszOuter.GetNpoints()>10) {
2784 lfabszOuter.EvalRobust(0.95);
2785
2786 TVectorD result(2);
2787 lfabszOuter.GetParameters(result);
2788 zAbsOutOffset = result[0];
2789 zAbsOutSlope = result[1];
2790 zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2791 }
2792 }
2793
be113f6e 2794
f6220766 2795 Int_t itime=-1;
2796 Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2797 Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2798 Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2799 Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2800 //
2801 Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2802 Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2803 //
2804 Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2805 Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2806
2b35e8e6 2807 //
2808 (*pcstream)<<"Mean"<<
be113f6e 2809 "run="<<run<< //
f6220766 2810 //voltages
2811 "UcIA="<<coverIA<<
2812 "UcIC="<<coverIC<<
2813 "UcOA="<<coverOA<<
2814 "UcOC="<<coverOC<<
2815 "UsA="<<skirtA<<
2816 "UsC="<<skirtC<<
2817 "UggA="<<ggOffA<<
2818 "UggC="<<ggOffC<<
2819 //
be113f6e 2820 "isOK="<<isOK<< //
2821 "id="<<id<< // track id
2822 "entries="<<entries<< // number of entries
2823 "bz="<<bfield<< // bfield
2b35e8e6 2824 "LTr.="<<ltrp<< // refernece track
1fd56785 2825 //
be113f6e 2826 "mnclIn="<<mnclIn<< // mean number of clusters in inner
2827 "mnclOut="<<mnclOut<< // mean number of clusters in outer
2828 "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
2829 "mLclIn="<<mLclIn<< // mean number of clusters in inner
2830 "mLclOut="<<mLclOut<< // mean number of clusters in outer
2831 "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
2832 "mdEdx="<<mdEdx<< // mean dEdx
108953e9 2833 "edZfit="<<edZfit<< // entries z fit
be113f6e 2834 "mdZfit="<<mdZfit<< // mean z fit
108953e9 2835 "rdZfit="<<rdZfit<< // RMS z fit
be113f6e 2836 //
2837 //
2838 "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
2839 "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
2840 "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
2841 "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
2842 "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
2843 "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
2844 "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
2845 "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
2846 "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
2847 "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
2848 //
2849 //
2850 "ePy1vP0="<<ePy1vP0<<
2851 "mPy1vP0="<<mPy1vP0<<
2852 "rPy1vP0="<<rPy1vP0<<
2853 "ePy2vP0="<<ePy2vP0<<
2854 "mPy2vP0="<<mPy2vP0<<
2855 "rPy2vP0="<<rPy2vP0<<
2856 "ePy3vP0="<<ePy3vP0<<
2857 "mPy3vP0="<<mPy3vP0<<
2858 "rPy3vP0="<<rPy3vP0<<
2859 "ePy1vP1="<<ePy1vP1<<
2860 "mPy1vP1="<<mPy1vP1<<
2861 "rPy1vP1="<<rPy1vP1<<
2862 "ePy2vP1="<<ePy2vP1<<
2863 "mPy2vP1="<<mPy2vP1<<
2864 "rPy2vP1="<<rPy2vP1<<
2865 "ePy3vP1="<<ePy3vP1<<
2866 "mPy3vP1="<<mPy3vP1<<
2867 "rPy3vP1="<<rPy3vP1<<
2868 "ePy2vP2In="<<ePy2vP2In<<
2869 "mPy2vP2In="<<mPy2vP2In<<
2870 "rPy2vP2In="<<rPy2vP2In<<
2871 "ePy2vP2Out="<<ePy2vP2Out<<
2872 "mPy2vP2Out="<<mPy2vP2Out<<
2873 "rPy2vP2Out="<<rPy2vP2Out<<
2874 "ePy3vP2IO="<<ePy3vP2IO<<
2875 "mPy3vP2IO="<<mPy3vP2IO<<
2876 "rPy3vP2IO="<<rPy3vP2IO<<
2877 //
2878 //
2879 "ePz1vP0="<<ePz1vP0<<
2880 "mPz1vP0="<<mPz1vP0<<
2881 "rPz1vP0="<<rPz1vP0<<
2882 "ePz2vP0="<<ePz2vP0<<
2883 "mPz2vP0="<<mPz2vP0<<
2884 "rPz2vP0="<<rPz2vP0<<
2885 "ePz3vP0="<<ePz3vP0<<
2886 "mPz3vP0="<<mPz3vP0<<
2887 "rPz3vP0="<<rPz3vP0<<
2888 "ePz1vP1="<<ePz1vP1<<
2889 "mPz1vP1="<<mPz1vP1<<
2890 "rPz1vP1="<<rPz1vP1<<
2891 "ePz2vP1="<<ePz2vP1<<
2892 "mPz2vP1="<<mPz2vP1<<
2893 "rPz2vP1="<<rPz2vP1<<
2894 "ePz3vP1="<<ePz3vP1<<
2895 "mPz3vP1="<<mPz3vP1<<
2896 "rPz3vP1="<<rPz3vP1<<
2897 "ePz2vP2In="<<ePz2vP2In<<
2898 "mPz2vP2In="<<mPz2vP2In<<
2899 "rPz2vP2In="<<rPz2vP2In<<
2900 "ePz2vP2Out="<<ePz2vP2Out<<
2901 "mPz2vP2Out="<<mPz2vP2Out<<
2902 "rPz2vP2Out="<<rPz2vP2Out<<
2903 "ePz3vP2IO="<<ePz3vP2IO<<
2904 "mPz3vP2IO="<<mPz3vP2IO<<
2905 "rPz3vP2IO="<<rPz3vP2IO<<
2906 //
2907 //
2908 //
2909 "lx0="<<lxyz[0]<< // reference x
2910 "lx1="<<lxyz[1]<< // reference y
2911 "lx2="<<lxyz[2]<< // refernece z
2912 "lpx0="<<lpxyz[0]<< // reference x
2b35e8e6 2913 "lpx1="<<lpxyz[1]<< // reference y
263d466a 2914 "lpx2="<<lpxyz[2]<< // refernece z
1fd56785 2915 //
2916 "msig="<<meanS<<
2917 //
2b35e8e6 2918 "mphi="<<meanphi<< //
2919 "rmsphi="<<rmsphi<< //
2920 "gphi1="<<gphi1<<
2921 "gphi2="<<gphi2<<
1fd56785 2922 //
2b35e8e6 2923 "mphiP="<<meanphiP<< //
2924 "rmsphiP="<<rmsphiP<< //
2925 "gphiP1="<<gphiP1<<
2926 "gphiP2="<<gphiP2<<
1fd56785 2927 //
2b35e8e6 2928 "meanZ="<<meanZ<<
2929 "rmsZ="<<rmsZ<<
2930 "gz1="<<gz1<<
2931 "gz2="<<gz2<<
e5c59218 2932 //
2933 "gp31="<<gp31<< //gaus mean - tgl
2934 "gp32="<<gp32<< //gaus rms -tgl
0b736a46 2935 "meanp3="<<meanp3<< //mean - tgl
2936 "rmsp3="<<rmsp3<< //rms -tgl
e5c59218 2937 "gp41="<<gp41<< //gaus mean - P4
2938 "gp42="<<gp42<< //gaus rms - P4
0b736a46 2939 "meanp4="<<meanp4<< //mean - P4
2940 "rmsp4="<<rmsp4<< //rms - P4
f6220766 2941 // Parameters from abs res analysis
2942 "SecIn="<<secInner<< // inner sector
2943 "SecOut="<<secOuter<< // outer sector
2944 "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2945 "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
fb962ae3 2946 "ibr.="<<&vecIBR<< // radial filed integral
2947 "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
b1c27e4f 2948 "ibr.="<<&vecIBR<< // radial filed integral
2949 "ibz.="<<&vecIBZ<< // z filed integral
2950 //
2951 "iblx.="<<&vecIBLX<< // local bx integral
2952 "ibly.="<<&vecIBLY<< // local by integral
2953 "ibgx.="<<&vecIBGX<< // global bx integral
2954 "ibgy.="<<&vecIBGY<< // global by integral
2955 //
f6220766 2956 "X.="<<&vecX<< // local x
2957 "Y.="<<&vecY<< // local y
2958 "R.="<<&vecR<< // radius
2959 "dY.="<<&vecDY<< // abs y residuals
2960 "dZ.="<<&vecDZ<< // abs z residuals
2961 "eY.="<<&vecEy<< // error of y residuals
2962 "eZ.="<<&vecEz<< // error of z residuals
2963 "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
2964 "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
2965 "nCl.="<<&vecN<< // number of clusters
2966 //
2967 "nClInY="<<nClInY<< // Number of clusters for inner y
2968 "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2969 "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
2970 "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2971 "nClInZ="<<nClInZ<< // Number of clusters for inner z
2972 "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2973 "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
2974 "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2975 //
2976 "nClOutY="<<nClOutY<< // Number of clusters for outer y
2977 "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
2978 "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
2979 "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
2980 "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
2981 "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
2982 "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
2983 "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
2984 //
2b35e8e6 2985 "\n";
2986 }
2987 delete pcstream;
be113f6e 2988 /*
2989 Browse the content
2990 TFile fmean("laserMean.root");
2991
2992
2993 */
2994
2995
2b35e8e6 2996}
2997
1fd56785 2998
2999
592a0c8f 3000void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
1fd56785 3001 //
3002 //
3003 //
3004 TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3005 TFile * f = pcstream->GetFile();
3006 f->mkdir("dirphi");
3007 f->mkdir("dirphiP");
3008 f->mkdir("dirZ");
3009 TF1 fp("p1","pol1");
3010 //
3011 //
3012 char cut[1000];
3013 char grname[1000];
3014 char grnamefull[1000];
3015
3016 Double_t mphi[100];
3017 Double_t mphiP[100];
3018 Double_t smphi[100];
3019 Double_t smphiP[100];
3020 Double_t mZ[100];
3021 Double_t smZ[100];
3022 Double_t bz[100];
3023 Double_t sbz[100];
3024 // fit parameters
3025 Double_t pphi[3];
3026 Double_t pphiP[3];
3027 Double_t pmZ[3];
0cd8bb8d 3028
1fd56785 3029 //
3030 for (Int_t id=0; id<336; id++){
3031 // id =205;
592a0c8f 3032 sprintf(cut,"fId==%d&&%s",id,cutUser);
1fd56785 3033 Int_t entries = chain->Draw("bz",cut,"goff");
3034 if (entries<3) continue;
e5c59218 3035 AliTPCLaserTrack *ltrp = 0;
1fd56785 3036 if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3037 ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3038 Double_t lxyz[3];
3039 Double_t lpxyz[3];
3040 ltrp->GetXYZ(lxyz);
3041 ltrp->GetPxPyPz(lpxyz);
3042
3043 chain->Draw("bz",cut,"goff");
3044 memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3045 chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3046 memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3047 //
3048 chain->Draw("gphi1",cut,"goff");
3049 memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
592a0c8f 3050 chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
1fd56785 3051 memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3052 //
3053 chain->Draw("gphiP1",cut,"goff");
3054 memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
592a0c8f 3055 chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
1fd56785 3056 memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3057 //
3058 chain->Draw("gz1",cut,"goff");
3059 memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
592a0c8f 3060 chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
1fd56785 3061 memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3062 //
3063 //
3064 sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3065 ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3066 // store data
3067 // phi
3068 f->cd("dirphi");
592a0c8f 3069 Float_t phiB0 =0;
3070 Float_t phiPB0=0;
3071 Float_t zB0=0;
3072 for (Int_t ientry=0; ientry<entries; ientry++){
3073 if (TMath::Abs(bz[ientry])<0.05){
3074 phiB0 = mphi[ientry];
3075 phiPB0 = mphiP[ientry];
3076 zB0 = mZ[ientry];
3077 }
3078 }
1fd56785 3079 TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3080 grphi->Draw("a*");
3081 grphi->Fit(&fp);
3082 pphi[0] = fp.GetParameter(0); // offset
3083 pphi[1] = fp.GetParameter(1); // slope
3084 pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3085 sprintf(grname,"phi_id%d",id);
3086 grphi->SetName(grname); grphi->SetTitle(grnamefull);
3087 grphi->GetXaxis()->SetTitle("b_{z} (T)");
3088 grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
ae69f16f 3089 grphi->SetMaximum(1.2);
3090 grphi->SetMinimum(-1.2);
1fd56785 3091 grphi->Draw("a*");
3092
3093 grphi->Write();
3094 gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3095 // phiP
3096 f->cd("dirphiP)");
3097 TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3098 grphiP->Draw("a*");
3099 grphiP->Fit(&fp);
3100 pphiP[0] = fp.GetParameter(0); // offset
3101 pphiP[1] = fp.GetParameter(1); // slope
3102 pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3103 sprintf(grname,"phiP_id%d",id);
3104 grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
3105 grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3106 grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
ae69f16f 3107 grphiP->SetMaximum(pphiP[0]+0.005);
3108 grphiP->SetMinimum(pphiP[0]-0.005);
1fd56785 3109
3110 gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3111 grphiP->Write();
3112 //
3113 //Z
3114 f->cd("dirZ");
3115 TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3116 grmZ->Draw("a*");
3117 grmZ->Fit(&fp);
3118 pmZ[0] = fp.GetParameter(0); // offset
3119 pmZ[1] = fp.GetParameter(1); // slope
3120 pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
3121 sprintf(grname,"mZ_id%d",id);
3122 grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
3123 grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3124 grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3125
3126 gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3127 grmZ->Write();
592a0c8f 3128 //
3129 // P4
3130 //
1fd56785 3131
3132 for (Int_t ientry=0; ientry<entries; ientry++){
3133 (*pcstream)<<"Mean"<<
3134 "id="<<id<<
3135 "LTr.="<<ltrp<<
3136 "entries="<<entries<<
3137 "bz="<<bz[ientry]<<
3138 "lx0="<<lxyz[0]<< // reference x
3139 "lx1="<<lxyz[1]<< // reference y
3140 "lx2="<<lxyz[2]<< // refernece z
3141 "lpx0="<<lpxyz[0]<< // reference x
3142 "lpx1="<<lpxyz[1]<< // reference y
3143 "lpx2="<<lpxyz[2]<< // refernece z
3144 //values
592a0c8f 3145 "phiB0="<<phiB0<< // position shift at 0 field
3146 "phiPB0="<<phiPB0<< // angular shift at 0 field
3147 "zB0="<<zB0<< // z shift for 0 field
3148 //
1fd56785 3149 "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3150 "pphi0="<<pphi[0]<< // offset
592a0c8f 3151 "pphi1="<<pphi[1]<< // slope
1fd56785 3152 "pphi2="<<pphi[2]<< // norm chi2
3153 //
3154 "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3155 "pphiP0="<<pphiP[0]<< // offset
592a0c8f 3156 "pphiP1="<<pphiP[1]<< // slope
1fd56785 3157 "pphiP2="<<pphiP[2]<< // norm chi2
3158 //
3159 "gz1="<<mZ[ientry]<<
3160 "pmZ0="<<pmZ[0]<< // offset
592a0c8f 3161 "pmZ1="<<pmZ[1]<< // slope
1fd56785 3162 "pmZ2="<<pmZ[2]<< // norm chi2
3163 "\n";
3164 }
3165 }
3166
3167 delete pcstream;
3168
3169}
3170
3171
c6914c83 3172void AliTPCcalibLaser::Analyze(){
3173 //
3174 //
3175 //
3176}
3177
3178
c03e3250 3179Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
c6914c83 3180
c03e3250 3181 TIterator* iter = li->MakeIterator();
3182 AliTPCcalibLaser* cal = 0;
108953e9 3183 static Int_t counter0=0;
c03e3250 3184 while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3185 if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
5de4d13e 3186 // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
c03e3250 3187 return -1;
3188 }
108953e9 3189 printf("Marging number %d\n", counter0);
3190 counter0++;
c03e3250 3191 //
be113f6e 3192 MergeFitHistos(cal);
c03e3250 3193 TH1F *h=0x0;
3194 TH1F *hm=0x0;
108953e9 3195 TH2F *h2=0x0;
3196 TH2F *h2m=0x0;
3197 // TProfile *hp=0x0;
3198 //TProfile *hpm=0x0;
c03e3250 3199
3200 for (Int_t id=0; id<336; id++){
3201 // merge fDeltaZ histograms
263d466a 3202 hm = (TH1F*)cal->fDeltaZ.At(id);
c03e3250 3203 h = (TH1F*)fDeltaZ.At(id);
3204 if (!h) {
3205 h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
5732d662 3206 h->SetDirectory(0);
c03e3250 3207 fDeltaZ.AddAt(h,id);
3208 }
3209 if (hm) h->Add(hm);
e5c59218 3210 // merge fP3 histograms
3211 hm = (TH1F*)cal->fDeltaP3.At(id);
3212 h = (TH1F*)fDeltaP3.At(id);
3213 if (!h) {
3214 h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3215 h->SetDirectory(0);
3216 fDeltaP3.AddAt(h,id);
3217 }
3218 if (hm) h->Add(hm);
3219 // merge fP4 histograms
3220 hm = (TH1F*)cal->fDeltaP4.At(id);
3221 h = (TH1F*)fDeltaP4.At(id);
3222 if (!h) {
3223 h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3224 h->SetDirectory(0);
3225 fDeltaP4.AddAt(h,id);
3226 }
3227 if (hm) h->Add(hm);
3228
3229 //
c03e3250 3230 // merge fDeltaPhi histograms
263d466a 3231 hm = (TH1F*)cal->fDeltaPhi.At(id);
c03e3250 3232 h = (TH1F*)fDeltaPhi.At(id);
3233 if (!h) {
3234 h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
5732d662 3235 h->SetDirectory(0);
c03e3250 3236 fDeltaPhi.AddAt(h,id);
3237 }
3238 if (hm) h->Add(hm);
3239 // merge fDeltaPhiP histograms
263d466a 3240 hm = (TH1F*)cal->fDeltaPhiP.At(id);
c03e3250 3241 h = (TH1F*)fDeltaPhiP.At(id);
3242 if (!h) {
3243 h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
5732d662 3244 h->SetDirectory(0);
c03e3250 3245 fDeltaPhiP.AddAt(h,id);
3246 }
3247 if (hm) h->Add(hm);
3248 // merge fSignals histograms
263d466a 3249 hm = (TH1F*)cal->fSignals.At(id);
c03e3250 3250 h = (TH1F*)fSignals.At(id);
3251 if (!h) {
e5c59218 3252 h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
5732d662 3253 h->SetDirectory(0);
c03e3250 3254 fSignals.AddAt(h,id);
3255 }
263d466a 3256 if (hm) h->Add(hm);
95a0e09b 3257 //
3258 //
108953e9 3259 // merge ProfileY histograms -0
3260 h2m = (TH2F*)cal->fDeltaYres.At(id);
3261 h2 = (TH2F*)fDeltaYres.At(id);
5de4d13e 3262 if (h2m&&h2) h2->Add(h2m);
95a0e09b 3263 //
108953e9 3264 h2m = (TH2F*)cal->fDeltaZres.At(id);
3265 h2 = (TH2F*)fDeltaZres.At(id);
5de4d13e 3266 if (h2m&&h2) h2->Add(h2m);
108953e9 3267 // merge ProfileY histograms - 2
3268 h2m = (TH2F*)cal->fDeltaYres2.At(id);
3269 h2 = (TH2F*)fDeltaYres2.At(id);
5de4d13e 3270 if (h2m&&h2) h2->Add(h2m);
108953e9 3271 //
3272 h2m = (TH2F*)cal->fDeltaZres2.At(id);
3273 h2 = (TH2F*)fDeltaZres2.At(id);
5de4d13e 3274 if (h2m&&h2) h2->Add(h2m);
7af539c6 3275
3276 // merge ProfileY histograms - abs
3277 h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3278 h2 = (TH2F*)fDeltaYresAbs.At(id);
3279 if (h2m&&h2) h2->Add(h2m);
3280
3281 h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3282 h2 = (TH2F*)fDeltaZresAbs.At(id);
3283 if (h2m&&h2) h2->Add(h2m);
108953e9 3284 // merge ProfileY histograms - 3
3285 //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3286 //h2 = (TH2F*)fDeltaYres3.At(id);
3287 //if (h2m) h2->Add(h2m);
3288 //
3289 //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3290 //h2 = (TH2F*)fDeltaZres3.At(id);
3291 //if (h2m) h->Add(h2m);
95a0e09b 3292 //
3293 //
c03e3250 3294 }
3295 }
3296 return 0;
3297}
3298
be113f6e 3299void AliTPCcalibLaser::MakeFitHistos(){
3300 //
3301 // Make a fit histograms
3302 //
3303 // Number of clusters
3304 //
3305 //TH2F *fHisNclIn; //->Number of clusters inner
3306 //TH2F *fHisNclOut; //->Number of clusters outer
3307 //TH2F *fHisNclIO; //->Number of cluster inner outer
3308 // TH2F *fHisdEdx; //->dEdx histo
3309 fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3310 fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3311 fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3312 //
3313 fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3314 fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3315 fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3316 //
3317 fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
108953e9 3318 fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
be113f6e 3319
3320 //
3321 // Chi2
3322 //
3323 // TH2F *fHisChi2YIn1; //->chi2 y inner - line
3324 // TH2F *fHisChi2YOut1; //->chi2 y inner - line
3325 // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
3326 // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
3327 // TH2F *fHisChi2YIO1; //->chi2 y IO - common
3328 fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3329 fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3330 fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3331 fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3332 fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3333 // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
3334 // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
3335 // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
3336 // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
3337 // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
3338 fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3339 fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3340 fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3341 fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3342 fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3343 //
3344 // Fit
3345 //
3346 //
3347 // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
3348 // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
3349 // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
3350 // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
3351 // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
3352 // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
3353 // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
3354 // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
3355 // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
3356 fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
3357 fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
3358 fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
3359 fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
3360 fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
3361 fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
3362 fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3363 fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3364 fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3365 //
3366 //
3367 // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
3368 // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
3369 // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
3370 // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
3371 // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
3372 // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
3373 // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
3374 // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
3375 // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
3376 fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
3377 fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
3378 fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
3379 fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
3380 fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
3381 fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
3382 fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3383 fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3384 fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3385
41a1f702 3386 fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3387 fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3388
be113f6e 3389 fHisNclIn->SetDirectory(0); //->Number of clusters inner
3390 fHisNclOut->SetDirectory(0); //->Number of clusters outer
3391 fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
3392 fHisLclIn->SetDirectory(0); //->Level arm inner
3393 fHisLclOut->SetDirectory(0); //->Level arm outer
3394 fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
3395 fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
3396 fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
3397 //
3398 //
3399 fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
3400 fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
3401 fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
3402 fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
3403 fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
3404 fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
3405 fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
3406 fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
3407 fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
3408 fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
3409 //
3410 //
3411 fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
3412 fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
3413 fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
3414 fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
3415 fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
3416 fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
3417 fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3418 fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3419 fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3420 //
3421 //
3422 fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
3423 fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
3424 fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
3425 fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
3426 fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
3427 fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
3428 fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
3429 fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
3430 fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
3431
41a1f702 3432 fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3433 fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3434
3435
be113f6e 3436
108953e9 3437 //
3438 //
3439 //
3440 for (Int_t id=0; id<336;id++){
3441 TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3442 TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3443 TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3444 TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
7af539c6 3445 TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3446 TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
108953e9 3447 // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3448 //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3449 if (!profy){
3450 profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3451 profy->SetDirectory(0);
3452 fDeltaYres.AddAt(profy,id);
3453 profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3454 profy2->SetDirectory(0);
3455 fDeltaYres2.AddAt(profy2,id);
41a1f702 3456 if(!fUseFixedDriftV)
3457 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
3458 else
3459 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 3460 profyabs->SetDirectory(0);
3461 fDeltaYresAbs.AddAt(profyabs,id);
108953e9 3462 //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3463 //profy3->SetDirectory(0);
3464 //fDeltaYres3.AddAt(profy3,id);
3465 }
3466 if (!profz){
3467 profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3468 profz->SetDirectory(0);
3469 fDeltaZres.AddAt(profz,id);
3470 profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3471 profz2->SetDirectory(0);
3472 fDeltaZres2.AddAt(profz2,id);
41a1f702 3473 if(!fUseFixedDriftV)
3474 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
3475 else
3476 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 3477 profzabs->SetDirectory(0);
3478 fDeltaZresAbs.AddAt(profzabs,id);
108953e9 3479 //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3480 //profz3->SetDirectory(0);
3481 //fDeltaZres3.AddAt(profz3,id);
3482 }
3483 }
3484 //
3485 //
3486 for (Int_t id=0; id<336;id++){
3487 TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3488 TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3489 TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3490
3491 TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3492 TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3493 TH1F * hisSignal = (TH1F*)fSignals.At(id);
3494
3495 if (!hisdz){
3496 hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3497 hisdz->SetDirectory(0);
3498 fDeltaZ.AddAt(hisdz,id);
3499
3500 hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3501 hisP3->SetDirectory(0);
3502 fDeltaP3.AddAt(hisP3,id);
3503 //
3504 hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3505 hisP4->SetDirectory(0);
3506 fDeltaP4.AddAt(hisP4,id);
be113f6e 3507
108953e9 3508 //
3509 hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3510 hisdphi->SetDirectory(0);
3511 fDeltaPhi.AddAt(hisdphi,id);
3512 //
3513 hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3514 hisdphiP->SetDirectory(0);
3515 fDeltaPhiP.AddAt(hisdphiP,id);
3516 hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3517 hisSignal->SetDirectory(0);
3518 fSignals.AddAt(hisSignal,id);
3519 }
3520 }
7af539c6 3521
3522 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3523 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3524 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3525 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
d3ce44cb 3526 //
3527 // Make THnSparse
3528 //
3529 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3530 Int_t binsLaser[12]= {336, //id
3531 2, //side
3532 6, //rod
3533 4, //bundle
3534 7, //beam
3535 300, //dP0
3536 300, //dP1
3537 300, //dP2
3538 300, //dP3
3539 300, //dP4
3540 80, //ncl
3541 50}; //dEdx
3542 Double_t xminLaser[12]= {0, //id
3543 0, //side
3544 0, //rod
3545 0, //bundle
3546 0, //beam
3547 -1, //dP0
3548 -1, //dP1
3549 -0.01, //dP2
3550 -0.01, //dP3
3551 -0.1, //dP4
3552 0, //ncl
3553 0}; //sqrt dEdx
3554 Double_t xmaxLaser[12]= {336, //id
3555 2, //side
3556 6, //rod
3557 4, //bundle
3558 7, //beam
3559 1, //dP0
3560 1, //dP1
3561 0.01, //dP2
3562 0.01, //dP3
3563 0.1, //dP4
3564 160, //ncl
3565 40}; //sqrt dEdx
3566
3567 TString nameLaser[12]= {"id",
3568 "side",
3569 "rod",
3570 "bundle",
3571 "beam",
3572 "dP0",
3573 "dP1",
3574 "dP2",
3575 "dP3",
3576 "dP4",
3577 "ncl",
3578 "sqrt dEdx"};
3579 TString titleLaser[12]= {"id",
3580 "side",
3581 "rod",
3582 "bundle",
3583 "beam",
3584 "#Delta_{P0}",
3585 "#Delta_{P1}",
3586 "#Delta_{P2}",
3587 "#Delta_{P3}",
3588 "#Delta_{P4}",
3589 "N_{cl}",
3590 "#sqrt{dEdx}"};
3591 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3592 for (Int_t iaxis=1; iaxis<12; iaxis++){
3593 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3594 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3595 }
be113f6e 3596}
3597
3598void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3599 //
3600 // Merge content of histograms
3601 //
3602 // Only first histogram is checked - all other should be the same
d3ce44cb 3603 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3604
be113f6e 3605 if (!laser->fHisNclIn) return; // empty histograms
3606 if (!fHisNclIn) MakeFitHistos();
3607 //
be113f6e 3608 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3609 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3610 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3611 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3612 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3613 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3614 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3615 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3616 //
3617 //
3618 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3619 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3620 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3621 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3622 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3623 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3624 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3625 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3626 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3627 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3628 //
3629 //
3630 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3631 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3632 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3633 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3634 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3635 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3636 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3637 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3638 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3639 //
3640 //
3641 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3642 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3643 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3644 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3645 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3646 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3647 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3648 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3649 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
41a1f702 3650 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3651 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
be113f6e 3652 //
3653 //
3654 //
3655
3656
3657
3658
3659}
3660
3661
3662
3663
e5c59218 3664void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3665 //
3666 // Dump fit information - collect information from the streamers
3667 //
3668 /*
3669 TChain * chainFit=0;
3670 TChain * chainTrack=0;
3671 TChain * chain=0;
3672 //
3673 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3674 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3675 AliXRDPROOFtoolkit tool;
3676 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3677 chainTrack->Lookup();
3678 chainTrack->SetProof(kTRUE);
be113f6e 3679 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3680 chainDrift->Lookup();
3681 chainDrift->SetProof(kTRUE);
3682
e5c59218 3683 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3684 chain->Lookup();
3685 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3686 chainFit->Lookup();
3687 chainFit->SetProof(kTRUE);
3688 chain->SetProof(kTRUE);
3689 AliTPCLaserTrack::LoadTracks();
3690 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3691
3692 */
3693 //
3694 // Fit cuts
3695 //
3696 TCut cutP3("abs(Tr.fP[3])<0.1");
3697 TCut cutP4("abs(Tr.fP[4])<0.5");
3698 TCut cutPx = cutP3+cutP4;
3699 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3700 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3701 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3702 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3703 //
3704 TCut cutdEdx("sqrt(dEdx)>3");
3705 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3706 TCut cutN("nclO>20&&nclI>20");
be113f6e 3707 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
e5c59218 3708 //
3709 // Cluster cuts
3710 //
3711 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3712 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3713 TCut cutClX("abs(Cl[].fX)>10");
3714 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3715 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3716 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3717 TCut cutQ("sqrt(Cl[].fMax)>4");
3718 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3719
3720
3721 TH1F * phisAl = 0;
3722 TH1F * phisAccept = 0;
3723 TH1F * phisOut = 0;
3724 TProfile * pdEdx = 0;
3725
3726 TProfile * pP0 = 0;
3727 TProfile * pP1 = 0;
3728 TProfile * pP2 = 0;
3729 TProfile * pP3 = 0;
3730 TProfile * pP4 = 0;
3731 //
3732 TProfile * pNclI = 0;
3733 TProfile * pNclO = 0;
3734 //
3735 TProfile * pchi2YIn =0;
3736 TProfile * pchi2ZIn =0;
3737 TProfile * pchi2YOut =0;
3738 TProfile * pchi2ZOut =0;
3739 TProfile * pchi2YInOut =0;
3740 TProfile * pchi2ZInOut =0;;
3741 // laser counters
3742 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3743 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3744 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3745 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3746 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3747 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3748 //
3749 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3750 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3751 // track param
3752 //
3753 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3754 pP0 = (TProfile*)gROOT->FindObject("hP0");
3755 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3756 pP1 = (TProfile*)gROOT->FindObject("hP1");
3757 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3758 pP2 = (TProfile*)gROOT->FindObject("hP2");
3759 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3760 pP3 = (TProfile*)gROOT->FindObject("hP3");
3761 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3762 pP4 = (TProfile*)gROOT->FindObject("hP4");
3763 //
3764 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3765 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3766 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3767 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3768 //
3769 //
3770 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3771 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3772 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3773 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3774 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3775 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3776 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3777 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3778 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3779 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3780 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3781 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3782 //
3783 // second derivatives
3784 //
3785 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3786 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3787 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3788 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3789 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3790 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3791 //
3792 phisPy2In->FitSlicesY();
3793 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3794 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3795 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3796 //
3797 phisPy2Out->FitSlicesY();
3798 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3799 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3800 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3801 //
3802 phisPy2InOut->FitSlicesY();
3803 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3804 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3805 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3806 //
3807 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3808 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3809 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3810 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3811 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3812 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3813 //
3814 phisPz2In->FitSlicesY();
3815 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3816 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3817 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3818 //
3819 phisPz2Out->FitSlicesY();
3820 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3821 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3822 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3823 //
3824 phisPz2InOut->FitSlicesY();
3825 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3826 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3827 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3828 //
3829 //
3830 //
3831
3832
3833 {
3834 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3835 for (Int_t ilaser=0; ilaser<336; ilaser++){
3836 Float_t all=phisAl->GetBinContent(ilaser+1);
3837 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3838 Float_t out=phisOut->GetBinContent(ilaser+1);
3839 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3840 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3841 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3842 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3843 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3844 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3845
3846
3847 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3848 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3849 //
3850 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3851 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3852 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3853 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3854 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3855 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3856 //
3857 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3858 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3859 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3860 //
3861 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3862 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3863 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3864 //
3865 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3866 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3867 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3868 //
3869 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3870 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3871 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3872 //
3873 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3874 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3875 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3876 //
3877 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3878 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3879 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3880
3881 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3882 (*pcstream)<<"Scan"<<
3883 "Run="<<id<<
3884 "LTr.="<<ltrp<<
3885 "all="<<all<<
3886 "accept="<<accept<<
3887 "out="<<out<<
3888 "sdedx="<<sdedx<<
3889 "mP0="<<mP0<<
3890 "mP1="<<mP1<<
3891 "mP2="<<mP2<<
3892 "mP3="<<mP3<<
3893 "mP4="<<mP4<<
3894 "nclI="<<nclI<<
3895 "nclO="<<nclO<<
3896 "chi2YIn="<<chi2YIn<<
3897 "chi2YOut="<<chi2YOut<<
3898 "chi2YInOut="<<chi2YInOut<<
3899 "chi2ZIn="<<chi2ZIn<<
3900 "chi2ZOut="<<chi2ZOut<<
3901 "chi2ZInOut="<<chi2ZInOut<<
3902 //
3903 "nPy2In="<<entriesPy2In<<
3904 "mPy2In="<<meanPy2In<<
3905 "sPy2In="<<sigmaPy2In<<
3906 //
3907 "nPy2Out="<<entriesPy2Out<<
3908 "mPy2Out="<<meanPy2Out<<
3909 "sPy2Out="<<sigmaPy2Out<<
3910 //
3911 "nPy2InOut="<<entriesPy2InOut<<
3912 "mPy2InOut="<<meanPy2InOut<<
3913 "sPy2InOut="<<sigmaPy2InOut<<
3914 //
3915 "nPz2In="<<entriesPz2In<<
3916 "mPz2In="<<meanPz2In<<
3917 "sPz2In="<<sigmaPz2In<<
3918 //
3919 "nPz2Out="<<entriesPz2Out<<
3920 "mPz2Out="<<meanPz2Out<<
3921 "sPz2Out="<<sigmaPz2Out<<
3922 //
3923 "nPz2InOut="<<entriesPz2InOut<<
3924 "mPz2InOut="<<meanPz2InOut<<
3925 "sPz2InOut="<<sigmaPz2InOut<<
3926 "\n";
3927 }
3928
3929 delete pcstream;
3930 }
7af539c6 3931}
3932
3933void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3934 TVectorD& meanSlope,
3935 TVectorD& sectorArray,
3936 Int_t option)
3937{
3938 // This method should ideally go in AliTPCLaser
3939 // option == 0 (pads outer - closest to beam)
3940 // option == 1 (pads inner)
3941 // option == 2 (time outer)
3942 // option == 3 (time inner)
3943 Int_t nFailures = 0;
3944
3945 for(Int_t id = 0; id < 336; id++) {
3946
3947 if (!AliTPCLaserTrack::GetTracks())
3948 AliTPCLaserTrack::LoadTracks();
3949 AliTPCLaserTrack *ltrp =
3950 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3951
3952 AliExternalTrackParam trackParam(*ltrp);
3953
3954 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3955 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3956 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3957
3958 Double_t angle = trackParam.GetAlpha();
3959 if(angle<0)
3960 angle += 2*TMath::Pi();
3961 if(trackParam.GetSnp()>0) // track points to sector "before"
3962 angle -= deltaangle*TMath::DegToRad();
3963 else // track points to sector "after"
3964 angle += deltaangle*TMath::DegToRad();
3965
3966 Bool_t success = trackParam.Rotate(angle);
3967
3968 if(!success) {
3969 // cout << "WARNING: Rotate failed for ID: " << id << endl;
3970 nFailures++;
3971 }
3972
3973 angle *= TMath::RadToDeg();
3974
3975 Int_t sector = TMath::Nint((angle-10.0)/20.0);
3976 if(sector<0)
3977 sector += 18;
3978 else if(sector>=18)
3979 sector -= 18;
3980 if(ltrp->GetSide()==1) // C side
3981 sector += 18;
f6220766 3982 if(option==0 || option==2)
3983 sector += 36;
3984 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
3985 sector += 36;
3986
7af539c6 3987 sectorArray[id] = sector;
3988
3989 const Double_t x0 = 0;
3990
3991 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3992 Double_t slopez = trackParam.GetTgl();
3993 // One needs a factor sqrt(1+slopey**2) to take into account the
3994 // longer path length
3995 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3996 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3997 slopez *= -1;
3998 // Double_t offsetz = trackParam.GetZ();
3999 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4000 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4001 if(option==2 || option==3) {
4002 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4003 } else {
4004 meanOffset[id] = offsety; meanSlope[id] = slopey;
4005 }
4006 }
4007
4008 if(nFailures>0)
4009 AliWarning(Form("Rotate method failed %d times", nFailures));
4010}
4011
4012
e5c59218 4013 /*
4014 TFile f("vscan.root");
4015 */
4016
4017 /*
4018 pad binning effect
4019 chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
4020 //
4021 chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
4022 //
4023
4024chain->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);
4025
4026
4027chain->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)
4028
4029 */
4030
4031
4032
4033
4034
4035 /*
4036 // check edge effects
4037 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4038 //
4039 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4040
4041 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)
4042
4043
4044
4045 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4046 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
4047
4048*/
4049
4050
4051
4052
4053
4054 /*
4055 Edge y effect
4056
4057 dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
4058
4059
4060 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)
4061
4062 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)
4063
4064
4065
4066
4067
4068 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)
4069
4070 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)
4071
4072
4073
4074 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)
4075
4076 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)
4077
4078
4079
4080*/
4081
4082
4083/*
4084
4085chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4086
4087chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4088
4089
4090
4091chainFit->Draw("LTr.fId","nclI>10",100000)
4092
4093chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4094
4095chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4096
4097TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4098
4099*/
4100
e5c59218 4101
4102
c03e3250 4103
4104
263d466a 4105
c03e3250 4106/*
4107 gSystem->Load("libSTAT.so")
4108 TStatToolkit toolkit;
4109 Double_t chi2;
4110 TVectorD fitParam;
4111 TMatrixD covMatrix;
4112 Int_t npoints;
e5c59218 4113
ae69f16f 4114 TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
c03e3250 4115
4116
4117TString fstring="";
1fd56785 4118//
4119fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++"; //1
4120fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++"; //2
4121fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++"; //3
4122fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++"; //4
4123//
4124fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++" //5
4125fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //6
4126fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++" //7
4127fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" //8
4128//
4129fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++" //9
4130fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //10
4131fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++" //11
4132fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" //12
c5595838 4133
c5595838 4134
c5595838 4135
c5595838 4136
1fd56785 4137 TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c03e3250 4138
1fd56785 4139 treeT->SetAlias("fit",strq0->Data());
4140
c03e3250 4141
1fd56785 4142 TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
c5595838 4143
1fd56785 4144 treeT->SetAlias("fitP",strqP->Data());
c5595838 4145
4146
ae69f16f 4147 TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4148 treeT->SetAlias("fitD",strqDrift->Data());
c5595838 4149
c03e3250 4150
ae69f16f 4151treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4152{
4153for (Int_t i=0; i<6;i++){
4154treeT->SetLineColor(i+2);
4155treeT->SetMarkerSize(1);
4156treeT->SetMarkerStyle(22+i);
4157treeT->SetMarkerColor(i+2);
c03e3250 4158
ae69f16f 4159treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4160}
4161}
c03e3250 4162 */
68ff0583 4163
4164
4165
4166/*
4167 TTree * tree = (TTree*)f.Get("FitModels");
4168
4169 TEventList listLFit0("listLFit0","listLFit0");
4170 TEventList listLFit1("listLFit1","listLFit1");
68ff0583 4171 tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4172 tree->SetEventList(&listLFit0);
4173
4174
95a0e09b 4175
4176
4177 gSystem->Load("libSTAT.so")
4178 TStatToolkit toolkit;
4179 Double_t chi2;
4180 TVectorD fitParam;
4181 TMatrixD covMatrix;
4182 Int_t npoints;
4183
4184 chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4185 chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4186
4187
4188 TString fstring="";
4189 fstring+="cos(dp)++";
4190 fstring+="sin(dp)++";
4191 fstring+="cos(dt)++";
4192 fstring+="sin(dt)++";
4193
4194 TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4195
4196
e5c59218 4197
4198*/
4199
4200
4201
4202/*
4203 Edge effects
4204 //
4205 //
4206 //
4207 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4208 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4209 AliXRDPROOFtoolkit tool;
4210 TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4211 chainTrack->Lookup();
4212 chainTrack->SetProof(kTRUE);
4213
4214 TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4215 chain->Lookup();
4216 TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4217 chainFit->Lookup();
4218 chainFit->SetProof(kTRUE);
4219 chain->SetProof(kTRUE);
4220 //
4221 // Fit cuts
4222 //
4223 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4224 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4225 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4226 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4227 //
4228 TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4229 TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4230 TCut cutN("nclO>20&&nclI>20");
4231 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4232 //
4233 // Cluster cuts
4234 //
4235 TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4236 TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4237 TCut cutClX("abs(Cl.fX)>10");
4238 TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4239 TCut cutCl=cutClY+cutClZ+cutClX;
4240
4241
4242 // check edge effects
4243 chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4244 //
4245 chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4246
4247 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)
4248
4249
4250
4251 chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4252 chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
4253
68ff0583 4254*/
e5c59218 4255