]> 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
1. JoinTrees Index - Possible to specify the branch without prefix ()
[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
994c9f5b 652 UpdateFitHistos();
880c3382 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);
994c9f5b 681 if (!hisdz) UpdateFitHistos();
108953e9 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 //==========================//
994c9f5b 2105 if (!fHisNclIn) UpdateFitHistos();
108953e9 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
994c9f5b 3376void AliTPCcalibLaser::MakeFitHistos(){
be113f6e 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
d3ce44cb 3601 //
3602 // Make THnSparse
3603 //
3604 // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
3605 Int_t binsLaser[12]= {336, //id
3606 2, //side
3607 6, //rod
3608 4, //bundle
3609 7, //beam
3610 300, //dP0
3611 300, //dP1
3612 300, //dP2
3613 300, //dP3
3614 300, //dP4
3615 80, //ncl
3616 50}; //dEdx
3617 Double_t xminLaser[12]= {0, //id
3618 0, //side
3619 0, //rod
3620 0, //bundle
3621 0, //beam
3622 -1, //dP0
3623 -1, //dP1
3624 -0.01, //dP2
3625 -0.01, //dP3
3626 -0.1, //dP4
3627 0, //ncl
3628 0}; //sqrt dEdx
3629 Double_t xmaxLaser[12]= {336, //id
3630 2, //side
3631 6, //rod
3632 4, //bundle
3633 7, //beam
3634 1, //dP0
3635 1, //dP1
3636 0.01, //dP2
3637 0.01, //dP3
3638 0.1, //dP4
3639 160, //ncl
3640 40}; //sqrt dEdx
3641
3642 TString nameLaser[12]= {"id",
3643 "side",
3644 "rod",
3645 "bundle",
3646 "beam",
3647 "dP0",
3648 "dP1",
3649 "dP2",
3650 "dP3",
3651 "dP4",
3652 "ncl",
3653 "sqrt dEdx"};
3654 TString titleLaser[12]= {"id",
3655 "side",
3656 "rod",
3657 "bundle",
3658 "beam",
3659 "#Delta_{P0}",
3660 "#Delta_{P1}",
3661 "#Delta_{P2}",
3662 "#Delta_{P3}",
3663 "#Delta_{P4}",
3664 "N_{cl}",
3665 "#sqrt{dEdx}"};
3666 fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3667 for (Int_t iaxis=1; iaxis<12; iaxis++){
3668 fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3669 fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3670 }
76c58ee2 3671 //
3672 // Delta Time bin
3673 // Pad SigmaShape Q charge pad row trackID
3674 Int_t binsRow[6]={200, 10000, 20, 30, 159, 336};
3675 Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
3676 Double_t axisMax[6]={ 1, 1000, 1, 30, 159, 336};
3677 TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
3678
3679 binsRow[1]=2000;
3680 axisMin[1]=0;
3681 axisMax[1]=200;
3682 fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3683 //
3684 binsRow[0]=1000;
3685 axisMin[0]=-20;
3686 axisMax[0]=20;
3687 binsRow[1]=10000;
3688 axisMin[1]=0;
3689 axisMax[1]=1000;
3690 //
3691 fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
3692 //
3693 for (Int_t iaxis=0; iaxis<6; iaxis++){
3694 fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
3695 fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
3696 }
be113f6e 3697}
3698
994c9f5b 3699void AliTPCcalibLaser::UpdateFitHistos(){
3700 //create the fit histos and set the beam parameters(needs OCDB access)
3701 MakeFitHistos();
3702 SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3703 SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3704 SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3705 SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3706}
3707
be113f6e 3708void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3709 //
3710 // Merge content of histograms
3711 //
3712 // Only first histogram is checked - all other should be the same
d3ce44cb 3713 if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
76c58ee2 3714 if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
3715 if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
3716 if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
3717 if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
3718
3719 if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
be113f6e 3720 if (!fHisNclIn) MakeFitHistos();
76c58ee2 3721 if (fHisNclIn->GetEntries()<1) MakeFitHistos();
be113f6e 3722 //
be113f6e 3723 fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
3724 fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
3725 fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
3726 fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
3727 fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
3728 fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
3729 fHisdEdx->Add(laser->fHisdEdx ); //->dedx
3730 fHisdZfit->Add(laser->fHisdZfit ); //->dedx
3731 //
3732 //
3733 fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
3734 fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
3735 fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
3736 fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
3737 fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
3738 fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
3739 fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
3740 fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
3741 fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
3742 fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
3743 //
3744 //
3745 fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
3746 fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
3747 fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
3748 fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
3749 fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
3750 fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
3751 fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
3752 fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
3753 fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
3754 //
3755 //
3756 fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
3757 fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
3758 fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
3759 fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
3760 fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
3761 fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
3762 fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
3763 fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
3764 fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
41a1f702 3765 fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3766 fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
be113f6e 3767 //
3768 //
3769 //
3770
3771
3772
3773
3774}
3775
3776
3777
3778
e5c59218 3779void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3780 //
3781 // Dump fit information - collect information from the streamers
3782 //
3783 /*
3784 TChain * chainFit=0;
3785 TChain * chainTrack=0;
3786 TChain * chain=0;
3787 //
3788 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3789 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3790 AliXRDPROOFtoolkit tool;
3791 chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3792 chainTrack->Lookup();
3793 chainTrack->SetProof(kTRUE);
be113f6e 3794 chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3795 chainDrift->Lookup();
3796 chainDrift->SetProof(kTRUE);
3797
e5c59218 3798 chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3799 chain->Lookup();
3800 chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3801 chainFit->Lookup();
3802 chainFit->SetProof(kTRUE);
3803 chain->SetProof(kTRUE);
3804 AliTPCLaserTrack::LoadTracks();
3805 //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3806
3807 */
3808 //
3809 // Fit cuts
3810 //
3811 TCut cutP3("abs(Tr.fP[3])<0.1");
3812 TCut cutP4("abs(Tr.fP[4])<0.5");
3813 TCut cutPx = cutP3+cutP4;
3814 TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3815 TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3816 TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3817 TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3818 //
3819 TCut cutdEdx("sqrt(dEdx)>3");
3820 TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3821 TCut cutN("nclO>20&&nclI>20");
be113f6e 3822 TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
e5c59218 3823 //
3824 // Cluster cuts
3825 //
3826 TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3827 TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3828 TCut cutClX("abs(Cl[].fX)>10");
3829 TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3830 TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3831 TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3832 TCut cutQ("sqrt(Cl[].fMax)>4");
3833 TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3834
3835
3836 TH1F * phisAl = 0;
3837 TH1F * phisAccept = 0;
3838 TH1F * phisOut = 0;
3839 TProfile * pdEdx = 0;
3840
3841 TProfile * pP0 = 0;
3842 TProfile * pP1 = 0;
3843 TProfile * pP2 = 0;
3844 TProfile * pP3 = 0;
3845 TProfile * pP4 = 0;
3846 //
3847 TProfile * pNclI = 0;
3848 TProfile * pNclO = 0;
3849 //
3850 TProfile * pchi2YIn =0;
3851 TProfile * pchi2ZIn =0;
3852 TProfile * pchi2YOut =0;
3853 TProfile * pchi2ZOut =0;
3854 TProfile * pchi2YInOut =0;
3855 TProfile * pchi2ZInOut =0;;
3856 // laser counters
3857 chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3858 phisAl = (TH1F*)gROOT->FindObject("hisAl");
3859 chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3860 phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3861 chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3862 phisOut = (TH1F*)gROOT->FindObject("hisOut");
3863 //
3864 chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3865 pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
3866 // track param
3867 //
3868 chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3869 pP0 = (TProfile*)gROOT->FindObject("hP0");
3870 chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3871 pP1 = (TProfile*)gROOT->FindObject("hP1");
3872 chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3873 pP2 = (TProfile*)gROOT->FindObject("hP2");
3874 chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3875 pP3 = (TProfile*)gROOT->FindObject("hP3");
3876 chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3877 pP4 = (TProfile*)gROOT->FindObject("hP4");
3878 //
3879 chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3880 pNclI = (TProfile*)gROOT->FindObject("hNclI");
3881 chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3882 pNclO = (TProfile*)gROOT->FindObject("hNclO");
3883 //
3884 //
3885 chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3886 pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
3887 chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3888 pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
3889 chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3890 pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
3891 chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3892 pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
3893 chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3894 pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
3895 chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3896 pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3897 //
3898 // second derivatives
3899 //
3900 TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3901 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3902 TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3903 chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3904 TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3905 chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3906 //
3907 phisPy2In->FitSlicesY();
3908 TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3909 TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3910 TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3911 //
3912 phisPy2Out->FitSlicesY();
3913 TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3914 TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3915 TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3916 //
3917 phisPy2InOut->FitSlicesY();
3918 TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3919 TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3920 TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3921 //
3922 TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3923 chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3924 TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3925 chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3926 TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3927 chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3928 //
3929 phisPz2In->FitSlicesY();
3930 TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3931 TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3932 TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3933 //
3934 phisPz2Out->FitSlicesY();
3935 TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3936 TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3937 TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3938 //
3939 phisPz2InOut->FitSlicesY();
3940 TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3941 TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3942 TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3943 //
3944 //
3945 //
3946
3947
3948 {
3949 TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3950 for (Int_t ilaser=0; ilaser<336; ilaser++){
3951 Float_t all=phisAl->GetBinContent(ilaser+1);
3952 Float_t accept=phisAccept->GetBinContent(ilaser+1);
3953 Float_t out=phisOut->GetBinContent(ilaser+1);
3954 Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3955 Float_t mP0 = pP0->GetBinContent(ilaser+1);
3956 Float_t mP1 = pP1->GetBinContent(ilaser+1);
3957 Float_t mP2 = pP2->GetBinContent(ilaser+1);
3958 Float_t mP3 = pP3->GetBinContent(ilaser+1);
3959 Float_t mP4 = pP4->GetBinContent(ilaser+1);
3960
3961
3962 Float_t nclI = pNclI->GetBinContent(ilaser+1);
3963 Float_t nclO = pNclO->GetBinContent(ilaser+1);
3964 //
3965 Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3966 Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3967 Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3968 Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3969 Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3970 Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3971 //
3972 Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
3973 Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
3974 Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
3975 //
3976 Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
3977 Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
3978 Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
3979 //
3980 Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
3981 Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
3982 Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
3983 //
3984 Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
3985 Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
3986 Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
3987 //
3988 Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
3989 Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
3990 Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
3991 //
3992 Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
3993 Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
3994 Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
3995
3996 AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3997 (*pcstream)<<"Scan"<<
3998 "Run="<<id<<
3999 "LTr.="<<ltrp<<
4000 "all="<<all<<
4001 "accept="<<accept<<
4002 "out="<<out<<
4003 "sdedx="<<sdedx<<
4004 "mP0="<<mP0<<
4005 "mP1="<<mP1<<
4006 "mP2="<<mP2<<
4007 "mP3="<<mP3<<
4008 "mP4="<<mP4<<
4009 "nclI="<<nclI<<
4010 "nclO="<<nclO<<
4011 "chi2YIn="<<chi2YIn<<
4012 "chi2YOut="<<chi2YOut<<
4013 "chi2YInOut="<<chi2YInOut<<
4014 "chi2ZIn="<<chi2ZIn<<
4015 "chi2ZOut="<<chi2ZOut<<
4016 "chi2ZInOut="<<chi2ZInOut<<
4017 //
4018 "nPy2In="<<entriesPy2In<<
4019 "mPy2In="<<meanPy2In<<
4020 "sPy2In="<<sigmaPy2In<<
4021 //
4022 "nPy2Out="<<entriesPy2Out<<
4023 "mPy2Out="<<meanPy2Out<<
4024 "sPy2Out="<<sigmaPy2Out<<
4025 //
4026 "nPy2InOut="<<entriesPy2InOut<<
4027 "mPy2InOut="<<meanPy2InOut<<
4028 "sPy2InOut="<<sigmaPy2InOut<<
4029 //
4030 "nPz2In="<<entriesPz2In<<
4031 "mPz2In="<<meanPz2In<<
4032 "sPz2In="<<sigmaPz2In<<
4033 //
4034 "nPz2Out="<<entriesPz2Out<<
4035 "mPz2Out="<<meanPz2Out<<
4036 "sPz2Out="<<sigmaPz2Out<<
4037 //
4038 "nPz2InOut="<<entriesPz2InOut<<
4039 "mPz2InOut="<<meanPz2InOut<<
4040 "sPz2InOut="<<sigmaPz2InOut<<
4041 "\n";
4042 }
4043
4044 delete pcstream;
4045 }
7af539c6 4046}
4047
4048void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
4049 TVectorD& meanSlope,
4050 TVectorD& sectorArray,
4051 Int_t option)
4052{
4053 // This method should ideally go in AliTPCLaser
4054 // option == 0 (pads outer - closest to beam)
4055 // option == 1 (pads inner)
4056 // option == 2 (time outer)
4057 // option == 3 (time inner)
4058 Int_t nFailures = 0;
4059
4060 for(Int_t id = 0; id < 336; id++) {
4061
4062 if (!AliTPCLaserTrack::GetTracks())
4063 AliTPCLaserTrack::LoadTracks();
4064 AliTPCLaserTrack *ltrp =
4065 (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4066
4067 AliExternalTrackParam trackParam(*ltrp);
4068
4069 Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
4070 if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
4071 deltaangle = 30.0; // inner sector is 1 sector further away from mirror
4072
4073 Double_t angle = trackParam.GetAlpha();
4074 if(angle<0)
4075 angle += 2*TMath::Pi();
4076 if(trackParam.GetSnp()>0) // track points to sector "before"
4077 angle -= deltaangle*TMath::DegToRad();
4078 else // track points to sector "after"
4079 angle += deltaangle*TMath::DegToRad();
4080
4081 Bool_t success = trackParam.Rotate(angle);
4082
4083 if(!success) {
4084 // cout << "WARNING: Rotate failed for ID: " << id << endl;
4085 nFailures++;
4086 }
4087
4088 angle *= TMath::RadToDeg();
4089
4090 Int_t sector = TMath::Nint((angle-10.0)/20.0);
4091 if(sector<0)
4092 sector += 18;
4093 else if(sector>=18)
4094 sector -= 18;
4095 if(ltrp->GetSide()==1) // C side
4096 sector += 18;
f6220766 4097 if(option==0 || option==2)
4098 sector += 36;
4099 if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
4100 sector += 36;
4101
7af539c6 4102 sectorArray[id] = sector;
4103
4104 const Double_t x0 = 0;
4105
4106 Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
4107 Double_t slopez = trackParam.GetTgl();
4108 // One needs a factor sqrt(1+slopey**2) to take into account the
4109 // longer path length
4110 slopez *= TMath::Sqrt(1.0 + slopey*slopey);
4111 if(fInverseSlopeZ) // wrong sign in database, should be fixed there
4112 slopez *= -1;
4113 // Double_t offsetz = trackParam.GetZ();
4114 Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4115 Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4116 if(option==2 || option==3) {
4117 meanOffset[id] = offsetz; meanSlope[id] = slopez;
4118 } else {
4119 meanOffset[id] = offsety; meanSlope[id] = slopey;
4120 }
4121 }
4122
4123 if(nFailures>0)
4124 AliWarning(Form("Rotate method failed %d times", nFailures));
4125}
4126
4127
e5c59218 4128
76c58ee2 4129void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
e5c59218 4130 //
4131 //
76c58ee2 4132 //input="TPCLaserObjects.root"
e5c59218 4133 //
76c58ee2 4134 // 0. OBJ: TAxis Delta
4135 // 1. OBJ: TAxis bin
4136 // 2. OBJ: TAxis rms shape
4137 // 3. OBJ: TAxis sqrt(Q)
4138 // 4. OBJ: TAxis row
4139 // 5. OBJ: TAxis trackID
e5c59218 4140
76c58ee2 4141 const Double_t kSigma=4.;
4142 TFile f(finput);
4143 AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
4144 THnSparse * hisPadInput = laserTPC->fHisLaserPad;
4145 THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
4146 TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
4147 TVectorD meanY(159), sigmaY(159);
4148 TVectorD meanZ(159), sigmaZ(159);
4149 TVectorD meanPad(159), sigmaPad(159);
4150 TVectorD meanTime(159), sigmaTime(159);
4151 TVectorD meanDPad(159), sigmaDPad(159);
4152 TVectorD meanDTime(159), sigmaDTime(159);
4153 TVectorD meandEdx(159), sigmadEdx(159);
4154 TVectorD meanSTime(159), sigmaSTime(159);
4155 TVectorD meanSPad(159), sigmaSPad(159);
4156 TVectorD entries(159);
e5c59218 4157 //
76c58ee2 4158 Int_t indexes[10]={0,1,2,3,4,5,6};
4159 TH1 *his=0;
4160 AliTPCLaserTrack::LoadTracks();
e5c59218 4161 //
76c58ee2 4162 for (Int_t id=0; id<336; id++){ // llop over laser beams
4163 printf("id=\t%d\n",id);
4164 //
4165 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4166 //
4167 hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
4168 hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
4169 //
4170 his=hisTimeInput->Projection(3);
4171 Int_t firstBindEdx=his->FindFirstBinAbove(0);
4172 Int_t lastBindEdx=his->FindLastBinAbove(0);
4173 hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4174 hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
4175 delete his;
4176 //
4177 his=hisTimeInput->Projection(1);
5b129619 4178 // Int_t firstBinTime=his->FindFirstBinAbove(0);
4179 //Int_t lastBinTime=his->FindLastBinAbove(0);
76c58ee2 4180 //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
4181 delete his;
4182 //
4183 //
4184 his=hisTimeInput->Projection(2);
5b129619 4185 //Int_t firstBinZ=his->FindFirstBinAbove(0);
4186 //Int_t lastBinZ=his->FindLastBinAbove(0);
76c58ee2 4187 //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
4188 delete his;
4189 //
4190 his=hisPadInput->Projection(2);
5b129619 4191 // Int_t firstBinY=his->FindFirstBinAbove(0);
4192 //Int_t lastBinY=his->FindLastBinAbove(0);
76c58ee2 4193 //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
4194 delete his;
4195 //
4196 //
4197 //
4198 THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
4199 THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
4200 //
4201 //
4202 for (Int_t irow=0; irow<159; irow++){
4203 entries[irow]=0;
4204 if ((*(ltrp->GetVecSec()))[irow] <0) continue;
4205 if ((*(ltrp->GetVecLX()))[irow] <80) continue;
4206
4207 hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
4208 hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
4209 //THnSparse *hisPad = hisPad0->Projection(4,indexes);
4210 //THnSparse *hisTime = hisTime0->Projection(4,indexes);
4211 THnSparse *hisPad = hisPad0;
4212 THnSparse *hisTime = hisTime0;
4213 //
4214 // Get mean value of QA variables
4215 //
4216 // dEdx
4217 his=hisTime->Projection(3);
4218 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4219 meandEdx[irow] =his->GetMean();
4220 sigmadEdx[irow]=his->GetRMS();
5b129619 4221 // Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
4222 //Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
76c58ee2 4223 // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
4224 //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
4225 delete his;
4226 //
4227 // sigma Time
4228 //
4229 his=hisTime->Projection(2);
4230 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
4231 meanSTime[irow] =his->GetMean();
4232 sigmaSTime[irow]=his->GetRMS();
5b129619 4233 //Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4234 //Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
76c58ee2 4235 // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
4236 delete his;
4237 //
4238 // sigma Pad
4239 his=hisPad->Projection(2);
4240 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4241 meanSPad[irow] =his->GetMean();
4242 sigmaSPad[irow]=his->GetRMS();
5b129619 4243 // Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
4244 //Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
76c58ee2 4245 // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
4246 delete his;
4247 //
4248 // apply selection on QA variables
4249 //
4250 //
4251 //
4252 // Y
4253 his=hisPad->Projection(0);
4254 entries[irow]=his->GetEntries();
4255 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4256 meanY[irow] =his->GetMean();
4257 sigmaY[irow]=his->GetRMS();
4258 delete his;
4259 // Z
4260 his=hisTime->Projection(0);
4261 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4262 meanZ[irow] =his->GetMean();
4263 sigmaZ[irow]=his->GetRMS();
4264 delete his;
4265 // Pad
4266 his=hisPad->Projection(1);
4267 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4268 meanPad[irow] =his->GetMean();
4269 meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
4270 sigmaPad[irow]=his->GetRMS();
4271 delete his;
4272 // Time
4273 his=hisTime->Projection(1);
4274 his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
4275 meanTime[irow] = his->GetMean();
4276 meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
4277 sigmaTime[irow]=his->GetRMS();
4278 delete his;
4279 //
4280 //delete hisTime;
4281 //delete hisPad;
4282 }
4283 //
4284 //
4285 //
4286 (*pcstream)<<"laserClusters"<<
4287 "id="<<id<<
4288 "run="<<run<<
4289 "LTr.="<<ltrp<<
4290 //
4291 "entries.="<<&entries<<
4292 "my.="<<&meanY<< //mean delta y
4293 "rmsy.="<<&sigmaY<< //rms deltay
4294 "mz.="<<&meanZ<< //mean deltaz
4295 "rmsz.="<<&sigmaZ<< //rms z
4296 //
4297 "mPad.="<<&meanPad<< // mean pad
4298 "mDPad.="<<&meanDPad<< // mead dpad
4299 "rmsPad.="<<&sigmaPad<< // rms pad
4300 "mTime.="<<&meanTime<<
4301 "mDTime.="<<&meanTime<<
4302 "rmsTime.="<<&sigmaTime<<
4303 //
4304 "mdEdx.="<<&meandEdx<< //mean dedx
4305 "rmsdEdx.="<<&sigmadEdx<< //rms dedx
4306 "mSPad.="<<&meanSPad<< //mean sigma pad
4307 "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
4308 "mSTime.="<<&meanSTime<< //mean sigma time
4309 "rmsSTime.="<<&sigmaSTime<<
4310 "\n";
4311 //
4312 delete hisPad0;
4313 delete hisTime0;
4314 }
4315 delete pcstream;
e5c59218 4316
76c58ee2 4317 /*
4318
4319 */
4320}
e5c59218 4321
76c58ee2 4322void AliTPCcalibLaser::FitLaserClusters(Int_t run){
4323 //
4324 //
4325 //input="TPCLaserObjects.root"
4326 //Algorithm:
4327 // 1. Select cluster candidates, remove outlyers
4328 // edge clusters
4329 // clusters with atypical spread (e.g due track overlaps)
4330 // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
4331 // 2. Fit the tracklets -per sector - in pad and time coordinate frame
4332 // Remove outlyers
4333 // Store info distance of track to pad, time center
4334 // Fit the correction for distance to the center (sin,cos)
4335 // 3. Do local fit
4336 const Double_t kEpsilon=0.000001;
4337 const Int_t kMinClusters=20;
4338 const Double_t kEdgeCut=3;
4339 const Double_t kDistCut=1.5; // cut distance to the ideal track
4340 const Double_t kDistCutFit=0.5;
4341 const Double_t kDistCutFitPad=0.25;
4342 const Double_t kDistCutFitTime=0.25;
5b129619 4343 const Int_t kSmoothRow=5;
76c58ee2 4344 TFile f("hisLasers.root"); // Input file
4345 TTree * treeInput=(TTree*)f.Get("laserClusters");
4346 TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
4347 TVectorD *vecN=0;
4348 TVectorD *vecMY=0;
4349 TVectorD *vecMZ=0;
4350 TVectorD *vecPad=0;
4351 TVectorD *vecTime=0;
4352 TVectorD *vecSY=0;
4353 TVectorD *vecSZ=0;
4354 TVectorD *meandEdx=0;
4355 TVectorD isOK(159);
4356 TVectorD fitPad(159);
4357 TVectorD fitTime(159);
4358 TVectorD fitPadLocal(159);
4359 TVectorD fitTimeLocal(159);
4360 TVectorD fitDPad(159);
4361 TVectorD fitDTime(159);
4362 TVectorD fitIPad(159);
4363 TVectorD fitITime(159);
4364 Double_t chi2PadIROC=0;
4365 Double_t chi2PadOROC=0;
4366 //
4367 treeInput->SetBranchAddress("my.",&vecMY);
4368 treeInput->SetBranchAddress("mz.",&vecMZ);
4369 treeInput->SetBranchAddress("mPad.",&vecPad);
4370 treeInput->SetBranchAddress("mTime.",&vecTime);
4371 treeInput->SetBranchAddress("rmsy.",&vecSY);
4372 treeInput->SetBranchAddress("rmsz.",&vecSZ);
4373 treeInput->SetBranchAddress("entries.",&vecN);
4374 treeInput->SetBranchAddress("mdEdx.",&meandEdx);
4375
4376 AliTPCLaserTrack::LoadTracks();
4377 //
4378 //
4379 TVectorD fitPadIROC(3), fitPadOROC(3);
4380 TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
4381 TVectorD fitTimeIROC(3), fitTimeOROC(3);
4382 TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
4383 //
4384 AliTPCROC * roc = AliTPCROC::Instance();
4385 Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
e5c59218 4386
76c58ee2 4387 //
4388 for (Int_t id=0; id<336; id++){
4389 //
4390 treeInput->GetEntry(id);
4391 AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
4392 Int_t medianEntries = TMath::Nint(TMath::Median(159,vecN->GetMatrixArray()));
4393 Double_t medianRMSY = TMath::Median(159,vecSY->GetMatrixArray());
4394 Double_t rmsRMSY = TMath::RMS(159,vecSY->GetMatrixArray());
4395 Double_t medianRMSZ = TMath::Median(159,vecSZ->GetMatrixArray());
4396 Double_t rmsRMSZ = TMath::RMS(159,vecSZ->GetMatrixArray());
4397 Double_t mdEdx = TMath::Median(159,meandEdx->GetMatrixArray());
4398 Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
4399 Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
4400 TLinearFitter fitterY(2,"pol1");
4401 TLinearFitter fitterZ(2,"pol1");
4402 TLinearFitter fitterPad(2,"pol1");
4403 TLinearFitter fitterTime(2,"pol1");
4404 TLinearFitter fitterPadSin(2,"hyp1");
4405 TLinearFitter fitterTimeSin(3,"hyp2");
4406 //
4407 //
4408 for (UInt_t irow=0; irow<159; irow++){
4409 fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
4410 fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
4411 Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
4412 isOK[irow]=kFALSE;
4413 fitPad[irow]=0;
4414 fitTime[irow]=0;
4415 Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
4416 Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
4417 (*vecPad)[irow]-=npads*0.5;
4418 //
4419 if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4420 if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4421 //
4422 if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
4423 if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
4424 if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
4425 if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
4426 if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
4427 if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
4428 if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
4429 Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
4430 if (TMath::Abs(dEdge)<kEdgeCut) continue;
4431 if (irow<roc->GetNRows(0)){
4432 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
4433 }
4434 if (irow>roc->GetNRows(0)){
4435 if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
4436 }
4437
4438 isOK[irow]=kTRUE;
4439 }
4440 //
4441 //fit OROC - get delta pad and delta time
4442 //
4443 fitterPad.ClearPoints();
4444 fitterTime.ClearPoints();
4445 fitterPadSin.ClearPoints();
4446 fitterTimeSin.ClearPoints();
4447 {for (Int_t irow=2; irow<157; irow++){
4448 if (isOK[irow]<0.5) continue;
4449 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4450 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4451 Double_t y=(*vecPad)[irow];
4452 Double_t z=(*vecTime)[irow];
4453 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4454 fitterPad.AddPoint(&x,y);
4455 fitterTime.AddPoint(&x,z);
4456 }}
4457 chi2PadOROC=0;
4458 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4459 fitterPad.Eval();
4460 fitterTime.Eval();
4461 chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4462 for (Int_t irow=2; irow<157; irow++){
4463 if (isOK[irow]<0.5) continue;
4464 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
4465 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4466 Double_t y=(*vecPad)[irow];
4467 Double_t z=(*vecTime)[irow];
4468 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4469 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4470 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4471 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4472 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4473 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4474 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4475 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4476 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4477 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4478 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4479 if (isOK[irow]>0){
4480 Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
4481 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4482 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4483 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4484 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4485 }
4486 }
4487 fitterPadSin.Eval();
4488 fitterTimeSin.Eval();
4489 fitterPadSin.FixParameter(0,0);
4490 fitterTimeSin.FixParameter(0,0);
4491 fitterPadSin.Eval();
4492 fitterTimeSin.Eval();
4493 //
4494 fitterPad.GetParameters(fitPadOROC);
4495 fitterTime.GetParameters(fitTimeOROC);
4496 fitterPadSin.GetParameters(fitPadOROCSin);
4497 fitterTimeSin.GetParameters(fitTimeOROCSin);
4498 }
4499 //
4500 //
4501 //fit IROC
4502 //
4503 fitterPad.ClearPoints();
4504 fitterTime.ClearPoints();
4505 fitterPadSin.ClearPoints();
4506 fitterTimeSin.ClearPoints();
4507 for (Int_t irow=2; irow<157; irow++){
4508 if (isOK[irow]<0.5) continue;
4509 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4510 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4511 Double_t y=(*vecPad)[irow];
4512 Double_t z=(*vecTime)[irow];
4513 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4514 fitterPad.AddPoint(&x,y);
4515 fitterTime.AddPoint(&x,z);
4516 }
4517 chi2PadIROC=0;
4518 if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
4519 fitterPad.Eval();
4520 fitterTime.Eval();
4521 chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
4522 for (Int_t irow=2; irow<157; irow++){
4523 if (isOK[irow]<0.5) continue;
4524 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
4525 if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
4526 Double_t y=(*vecPad)[irow];
4527 Double_t z=(*vecTime)[irow];
4528 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
4529 Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4530 Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4531 fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
4532 fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
4533 fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
4534 fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
4535 fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
4536 fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
4537 if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
4538 if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
4539 if (isOK[irow]>0.5){
4540 Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
4541 TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
4542 Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
4543 TMath::Cos(2*TMath::Pi()*fitITime[irow])};
4544 fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
4545 fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
4546 }
4547 }
4548 fitterPadSin.Eval();
4549 fitterTimeSin.Eval();
4550 fitterPadSin.FixParameter(0,0);
4551 fitterTimeSin.FixParameter(0,0);
4552 fitterPadSin.Eval();
4553 fitterTimeSin.Eval();
4554 fitterPad.GetParameters(fitPadIROC);
4555 fitterTime.GetParameters(fitTimeIROC);
4556 fitterPadSin.GetParameters(fitPadIROCSin);
4557 fitterTimeSin.GetParameters(fitTimeIROCSin);
4558 }
4559 for (Int_t irow=0; irow<159; irow++){
4560 if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
4561 if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
4562 if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
4563 if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
4564 }
4565 for (Int_t irow=kSmoothRow/2; irow<159-kSmoothRow/2; irow++){
4566 fitPadLocal[irow]=0;
4567 fitTimeLocal[irow]=0;
4568 if (isOK[irow]<0.5) continue;
5b129619 4569 Int_t sector=(irow<Int_t(roc->GetNRows(0)))? sectorInner:sectorOuter;
76c58ee2 4570 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
4571 //
4572 TLinearFitter fitterPadLocal(2,"pol1");
4573 TLinearFitter fitterTimeLocal(2,"pol1");
4574 Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
4575 for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
4576 Int_t jrow=irow+delta;
4577 if (jrow<0) jrow=0;
4578 if (jrow>159) jrow=159;
4579 if (isOK[jrow]<0.5) continue;
4580 if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
4581 Double_t y=(*vecPad)[jrow];
4582 Double_t z=(*vecTime)[jrow];
4583 Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
4584 fitterPadLocal.AddPoint(&x,y);
4585 fitterTimeLocal.AddPoint(&x,z);
4586 }
4587 if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
4588 fitterPadLocal.Eval();
4589 fitterTimeLocal.Eval();
4590 fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
4591 fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
4592 }
4593 //
4594 //
4595 (*pcstream)<<"fit"<<
4596 "run="<<run<<
4597 "id="<<id<<
4598 "chi2PadIROC="<<chi2PadIROC<<
4599 "chi2PadOROC="<<chi2PadOROC<<
4600 "mdEdx="<<mdEdx<<
4601 "LTr.="<<ltrp<<
4602 "isOK.="<<&isOK<<
4603 // mean measured-ideal values
4604 "mY.="<<vecMY<<
4605 "mZ.="<<vecMZ<<
4606 // local coordinate fit
4607 "mPad.="<<vecPad<<
4608 "mTime.="<<vecTime<<
4609 "fitPad.="<<&fitPad<<
4610 "fitTime.="<<&fitTime<<
4611 "fitPadLocal.="<<&fitPadLocal<<
4612 "fitTimeLocal.="<<&fitTimeLocal<<
4613 "fitDPad.="<<&fitDPad<<
4614 "fitDTime.="<<&fitDTime<<
4615 "fitIPad.="<<&fitIPad<<
4616 "fitITime.="<<&fitITime<<
4617 //
4618 "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
4619 "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
4620 "fitPadOROC.="<<&fitPadOROC<<
4621 "fitPadOROCSin.="<<&fitPadOROCSin<<
4622 //
4623 "fitTimeIROC.="<<&fitTimeIROC<<
4624 "fitTimeIROCSin.="<<&fitTimeIROCSin<<
4625 "fitTimeOROC.="<<&fitTimeOROC<<
4626 "fitTimeOROCSin.="<<&fitTimeOROCSin<<
4627 "\n";
4628 }
4629 delete pcstream;
4630}