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