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