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