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