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