Test for Coverity
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibLaser.cxx
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 /*
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)
25   //
26   // 3. The drift velocity and jitter is calculated event by event
27   //    (see function drift velocity)
28   //
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 
34   //
35   // To make laser scan the user interaction neccessary
36   //
37   .x ~/NimStyle.C
38   gSystem->Load("libANALYSIS");
39   gSystem->Load("libTPCcalib");
40   TFile fcalib("CalibObjectsTrain2.root");
41   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
42   laser->DumpMeanInfo(run)
43   TFile fmean("laserMean.root")
44   //
45   //  laser track clasification;
46   //
47   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
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;
51   TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
52
53   //
54   //
55   // Analyze  LASER scan 
56   //
57   
58   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
59   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
60   AliXRDPROOFtoolkit tool;
61   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
62   chain->Lookup();
63   AliTPCcalibLaser::DumpScanInfo(chain)
64   TFile fscan("laserScan.root")
65   TTree * treeT = (TTree*)fscan.Get("Mean")
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;
72   AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1) 
73   TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
74   chainDrift->Lookup();
75   TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
76   chainDriftN->Lookup();
77
78
79   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
80   chain->Lookup();
81   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
82   chainFit->Lookup();
83   TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
84   chainTrack->Lookup();
85
86 */
87
88
89
90 #include "TLinearFitter.h"
91 #include "AliTPCcalibLaser.h"
92 #include "AliExternalTrackParam.h"
93 #include "AliESDEvent.h"
94 #include "AliESDfriend.h"
95 #include "AliESDtrack.h"
96 #include "AliTPCTracklet.h"
97 #include "TH1D.h"
98 #include "TH1F.h"
99 #include "TProfile.h"
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"
108 #include "AliLog.h"
109 #include "TClonesArray.h"
110 #include "TPad.h"
111 #include "TSystem.h"
112 #include "TCut.h"
113 #include "TChain.h"
114 #include "TH2F.h"
115 #include "TStatToolkit.h"
116 #include "TROOT.h"
117 #include "TDatabasePDG.h"
118
119
120 #include "TTreeStream.h"
121 #include <iostream>
122 #include <sstream>
123 #include "AliTPCLaserTrack.h"
124 #include "AliTPCcalibDB.h"
125 #include "AliTPCParam.h"
126 #include "AliTPCParamSR.h"
127 #include "TTimeStamp.h"
128 #include "AliDCSSensorArray.h"
129 #include "AliDCSSensor.h"
130 #include "AliGRPObject.h"
131 #include "AliTPCROC.h"
132
133 using namespace std;
134
135 ClassImp(AliTPCcalibLaser)
136
137 AliTPCcalibLaser::AliTPCcalibLaser():
138   AliTPCcalibBase(),
139   fESD(0),
140   fESDfriend(0),
141   fNtracks(0),
142   fTracksMirror(336),
143   fTracksEsd(336),
144   fTracksEsdParam(336),
145   fTracksTPC(336),
146   fFullCalib(kTRUE),
147   fDeltaZ(336),
148   fDeltaP3(336),
149   fDeltaP4(336),
150   fDeltaPhi(336),
151   fDeltaPhiP(336),
152   fSignals(336),
153   //
154   fHisLaser(0),      //  N dim histogram of laser
155   fHisLaserPad(0),      //  N dim histogram of laser
156   fHisLaserTime(0),      //  N dim histogram of laser
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   //
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
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
208   //fDeltaYres3(336),   //->2D histo of residuals
209   //fDeltaZres3(336),   //->2D histo fo residuals
210   fFitAside(new TVectorD(5)),
211   fFitCside(new TVectorD(5)),      
212   fFitACside(new TVectorD(6)),      
213   fEdgeXcuts(3),    
214   fEdgeYcuts(3),    
215   fNClCuts(5),      
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),
227   fInverseSlopeZ(kTRUE),
228   fUseFixedDriftV(0),
229   fFixedFitAside0(0.0),
230   fFixedFitAside1(1.0),
231   fFixedFitCside0(0.0),
232   fFixedFitCside1(1.0)
233 {
234   //
235   // Constructor
236   //
237   fTracksEsdParam.SetOwner(kTRUE);
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   }
244 }
245
246 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
247   AliTPCcalibBase(),
248   fESD(0),
249   fESDfriend(0),
250   fNtracks(0),
251   fTracksMirror(336),
252   fTracksEsd(336),
253   fTracksEsdParam(336),
254   fTracksTPC(336),
255   fFullCalib(full),
256   //
257   fDeltaZ(336),          // array of histograms of delta z for each track
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
261   fDeltaPhiP(336),          // array of histograms of delta z for each track
262   fSignals(336),           // array of dedx signals
263   //
264   //
265   fHisLaser(0),      //  N dim histogram of laser
266   fHisLaserPad(0),      //  N dim histogram of laser
267   fHisLaserTime(0),      //  N dim histogram of laser
268
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   //
314   fDeltaYres(336),
315   fDeltaZres(336),  
316   fDeltaYres2(336),
317   fDeltaZres2(336),  
318   fDeltaYresAbs(336),
319   fHisYAbsErrors(0),
320   fDeltaZresAbs(336),
321   fHisZAbsErrors(0),
322   //  fDeltaYres3(336),
323   //fDeltaZres3(336),  
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
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
329   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
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),
341   fInverseSlopeZ(kTRUE),
342   fUseFixedDriftV(0),
343   fFixedFitAside0(0.0),
344   fFixedFitAside1(1.0),
345   fFixedFitCside0(0.0),
346   fFixedFitCside1(1.0)
347 {
348   SetName(name);
349   SetTitle(title);
350   //
351   // Constructor
352   //
353   fTracksEsdParam.SetOwner(kTRUE);
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   }
360 }
361
362 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
363   AliTPCcalibBase(calibLaser), 
364   fESD(0),
365   fESDfriend(0),
366   fNtracks(0),
367   fTracksMirror(336),
368   fTracksEsd(336),
369   fTracksEsdParam(336),
370   fTracksTPC(336),
371   fFullCalib(calibLaser.fFullCalib),
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   //
381   fHisLaser(0),      //  N dim histogram of laser
382   fHisLaserPad(0),      //  N dim histogram of laser
383   fHisLaserTime(0),      //  N dim histogram of laser
384
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))),  
431   fDeltaYres2(((calibLaser.fDeltaYres))),
432   fDeltaZres2(((calibLaser.fDeltaZres))),  
433   fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
434   fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
435   fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
436   fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
437   //  fDeltaYres3(((calibLaser.fDeltaYres))),
438   //fDeltaZres3(((calibLaser.fDeltaZres))),  
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
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
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),
456   fInverseSlopeZ(calibLaser.fInverseSlopeZ),
457   fUseFixedDriftV(calibLaser.fUseFixedDriftV),
458   fFixedFitAside0(calibLaser.fFixedFitAside0),
459   fFixedFitAside1(calibLaser.fFixedFitAside1),
460   fFixedFitCside0(calibLaser.fFixedFitCside0),
461   fFixedFitCside1(calibLaser.fFixedFitCside1)
462 {
463   //
464   // copy constructor
465   //
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   }
472 }
473
474
475
476 AliTPCcalibLaser & 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
490 AliTPCcalibLaser::~AliTPCcalibLaser() {
491   //
492   // destructor
493   //
494   if ( fHisNclIn){
495     delete fHisLaser;      //->
496     delete fHisLaserPad;      //->
497     delete fHisLaserTime;      //->
498
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   }
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();
561   delete fHisYAbsErrors;
562   fDeltaZres.SetOwner();
563   fDeltaZres.Delete();
564   delete fHisZAbsErrors;
565   fDeltaYres2.SetOwner();
566   fDeltaYres2.Delete();
567   fDeltaZres2.SetOwner();
568   fDeltaZres2.Delete();
569   
570   fDeltaYresAbs.SetOwner();
571   fDeltaYresAbs.Delete();
572   fDeltaZresAbs.SetOwner();
573   fDeltaZresAbs.Delete();
574 }
575
576
577
578 void AliTPCcalibLaser::Process(AliESDEvent * event) {
579   //
580   //
581   // Loop over tracks and call  Process function
582   //
583   const Int_t  kMinTracks=20;
584   const Int_t  kMinClusters=40;
585
586   fESD = event;
587   if (!fESD) {
588     return;
589   }
590   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
591   if (!fESDfriend) {
592     return;
593   }
594   if (fESDfriend->TestSkipBit()) return;
595   if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
596   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
597   //
598   // find CE background if present
599   //
600   if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
601   TH1D hisCE("hhisCE","hhisCE",100,-100,100);  
602   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
603     AliESDtrack *track=fESD->GetTrack(i);
604     if (!track) continue;
605     hisCE.Fill(track->GetZ());
606     hisCE.Fill(track->GetZ()+2);
607     hisCE.Fill(track->GetZ()-2);
608   }
609   //
610   //
611
612
613   fTracksTPC.Clear();
614   fTracksEsd.Clear();
615   fTracksEsdParam.Delete();
616   for (Int_t id=0; id<336;id++) {
617     fCounter[id]=0;
618     fClusterCounter[id]=0;
619     fClusterSatur[id]=0;
620   }
621   //
622   Int_t n=fESD->GetNumberOfTracks();
623   Int_t counter=0;
624   for (Int_t i=0;i<n;++i) {
625     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
626     if (!friendTrack) continue;
627     AliESDtrack *track=fESD->GetTrack(i);
628     if (!track) continue;
629     Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
630     if (binC>336) continue; //remove CE background
631     TObject *calibObject=0;
632     AliTPCseed *seed=0;
633     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
634       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
635         break;
636     if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
637       //filter CE tracks
638       Int_t id = FindMirror(track,seed);
639       if (id>=0) counter++;
640     }
641     //
642   } 
643   fNtracks=counter;
644   if (counter<kMinTracks) return;
645
646   //FitDriftV();
647   FitDriftV(0.2);
648   if (!fFullCalib) return;
649   static Bool_t init=kFALSE;
650   if (!init){
651     init = kTRUE;  // way around for PROOF - to be investigated
652     MakeFitHistos();
653   }
654   //
655   for (Int_t id=0; id<336; id++){    
656     //
657     if (!fTracksEsdParam.At(id)) continue;
658     if (fClusterSatur[id]>0.3)   continue;  // tracks in saturation
659     DumpLaser(id);
660     if ( AcceptLaser(id) && fFitZ[id]<0.5){
661       RefitLaserJW(id);    
662       MakeDistHisto(id);
663     }
664   }
665
666 }
667
668 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
669   //
670   //
671   //
672   //  for (Int_t id=0; id<336; id++){
673     //
674     //
675     if (!fTracksEsdParam.At(id)) return;
676     if (!AcceptLaser(id)) return;
677     Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
678     //
679     //
680     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
681     if (!hisdz) MakeFitHistos();
682     hisdz = (TH1F*)fDeltaZ.At(id);
683     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
684     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
685     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
686     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
687     TH1F * hisSignal = (TH1F*)fSignals.At(id);
688     //
689
690     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
691     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
692     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
693     if (!param) return;
694     if (!ltrp) return;
695     if (!track) return;
696     Double_t xyz[3];
697     Double_t pxyz[3];
698     Double_t lxyz[3];
699     Double_t lpxyz[3];
700     param->GetXYZ(xyz);
701     param->GetPxPyPz(pxyz);
702     ltrp->GetXYZ(lxyz);
703     ltrp->GetPxPyPz(lpxyz);
704     //
705     Float_t dz   = param->GetZ()-ltrp->GetZ();
706     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
707     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
708     if (hisdz) hisdz->Fill(dz);
709     if (hisP3) hisP3->Fill(param->GetParameter()[3]);
710     if (hisP4) hisP4->Fill(param->GetParameter()[4]);
711     if (hisdphi) hisdphi->Fill(dphi);
712     if (hisdphiP) hisdphiP->Fill(dphiP);
713     if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
714     // fill HisLaser
715     xhis[0] = ltrp->GetId();
716     xhis[1] = ltrp->GetSide();
717     xhis[2] = ltrp->GetRod();
718     xhis[3] = ltrp->GetBundle();
719     xhis[4] = ltrp->GetBeam();
720     xhis[5] = dphi;
721     xhis[6] = fFitZ[id];
722     xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
723     xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
724     xhis[9] = param->GetParameter()[4];
725     xhis[10]= track->GetTPCNcls();
726     xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
727     // }
728     fHisLaser->Fill(xhis);
729     //
730     
731 }
732
733 void AliTPCcalibLaser::FitDriftV(){
734   //
735   // Fit corrections to the drift velocity - linear approximation in the z and global y
736   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
737   // 
738   /*
739     Formulas:
740     
741     z  = s* (z0 - vd*(t-t0))
742     
743     s  - side -1 and +1 
744     t0 - time 0
745     vd - nominal drift velocity
746     zs - miscalibrated position
747     
748     zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
749     vr  - relative change of the drift velocity
750     dzt - vd*dt
751     dr  = zz0-s*z
752     ..
753     ==>
754     zs ~ z - s*vr*(z0-s*z)+s*dzt
755     --------------------------------
756     1. Correction function vr constant:
757     
758     
759     dz = zs-z = -s*vr *(z0-s*z)+s*dzt         
760     dzs/dl = dz/dl +s*s*vr*dz/dl 
761     d(dz/dl) = vr*dz/dl     
762   */
763
764   const Float_t kZCut        = 240;    // remove the closest laser beam
765   const Float_t kSaturCut    = 0.05;   // remove saturated lasers - cut on fraction of saturated 
766   const Float_t kDistCut     = 3;    // distance sigma cut
767   const Float_t kDistCutAbs  = 0.25;  
768   const Float_t kMinClusters = 60;  // minimal amount of the clusters
769   const Float_t kMinSignal   = 16;  // minimal mean height of the signal
770   const Float_t kChi2Cut     = 0.1; // chi2 cut to accept drift fit
771   static TLinearFitter fdriftA(3,"hyp2");
772   static TLinearFitter fdriftC(3,"hyp2");
773   static TLinearFitter fdriftAC(4,"hyp3");
774   TVectorD fitA(3),fitC(3),fitAC(4);
775   
776   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
777   AliTPCParam  * tpcparam    = calib->GetParameters();    
778   //
779   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
780
781   Float_t chi2A = 10;
782   Float_t chi2C = 10;
783   Float_t chi2AC = 10;
784   Int_t npointsA=0;
785   Int_t npointsC=0;
786   Int_t npointsAC=0;
787
788
789   for (Int_t iter=0; iter<3; iter++){
790     fdriftA.ClearPoints();
791     fdriftC.ClearPoints();
792     fdriftAC.ClearPoints();
793     //
794     for (Int_t id=0; id<336; id++){
795       if (!fTracksEsdParam.At(id)) continue;
796       if (!AcceptLaser(id)) continue;
797       if ( fClusterSatur[id]>kSaturCut)  continue;
798       if ( fClusterCounter[id]<kMinClusters)  continue;
799       AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
800       if (track->GetTPCsignal()<kMinSignal) continue;
801       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
802       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
803
804       Double_t xyz[3];
805       Double_t pxyz[3];
806       Double_t lxyz[3];
807       Double_t lpxyz[3];
808       param->GetXYZ(xyz);
809       param->GetPxPyPz(pxyz);
810       ltrp->GetXYZ(lxyz);
811       ltrp->GetPxPyPz(lpxyz);
812       if (TMath::Abs(lxyz[2])>kZCut) continue;
813       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
814       if (npointsAC>0) sz =TMath::Sqrt(chi2AC); 
815       if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
816       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
817
818       // drift distance
819       Double_t zlength =  tpcparam->GetZLength(0);
820       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
821       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
822       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
823       if (ltrp->GetSide()==0){
824         fdriftA.AddPoint(xxx,mdrift,1);
825       }else{
826         fdriftC.AddPoint(xxx,mdrift,1);
827       }
828       Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
829       fdriftAC.AddPoint(xxx2,mdrift,1);
830     }
831     //
832     if (fdriftA.GetNpoints()>10){
833       //
834       fdriftA.Eval();
835       if (iter==0) fdriftA.EvalRobust(0.7);
836       else fdriftA.EvalRobust(0.8);
837       fdriftA.GetParameters(fitA);
838       npointsA= fdriftA.GetNpoints();
839       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
840       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
841         if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
842         (*fFitAside)[0] = fitA[0];
843         (*fFitAside)[1] = fitA[1];
844         (*fFitAside)[2] = fitA[2];
845         (*fFitAside)[3] = fdriftA.GetNpoints();
846         (*fFitAside)[4] = chi2A;        
847       }
848     }
849     if (fdriftC.GetNpoints()>10){
850       fdriftC.Eval();
851       if (iter==0) fdriftC.EvalRobust(0.7);
852       else fdriftC.EvalRobust(0.8);
853       //
854       fdriftC.GetParameters(fitC);
855       npointsC= fdriftC.GetNpoints();
856       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
857       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) { 
858         if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
859         (*fFitCside)[0] = fitC[0];
860         (*fFitCside)[1] = fitC[1];
861         (*fFitCside)[2] = fitC[2];
862         (*fFitCside)[3] = fdriftC.GetNpoints();
863         (*fFitCside)[4] = chi2C;        
864       }
865     }
866
867     if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
868       fdriftAC.Eval();
869       if (iter==0) fdriftAC.EvalRobust(0.7);
870       else fdriftAC.EvalRobust(0.8);
871       //
872       fdriftAC.GetParameters(fitAC);
873       npointsAC= fdriftAC.GetNpoints();
874       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
875       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
876     }
877     
878     for (Int_t id=0; id<336; id++){
879       if (!fTracksEsdParam.At(id)) continue;
880       //
881       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
882       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
883       Double_t xyz[3];
884       Double_t pxyz[3];
885       Double_t lxyz[3];
886       Double_t lpxyz[3];
887       param->GetXYZ(xyz);
888       param->GetPxPyPz(pxyz);
889       ltrp->GetXYZ(lxyz);
890       ltrp->GetPxPyPz(lpxyz); 
891       Double_t zlength =  tpcparam->GetZLength(0);
892       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
893       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
894
895       Float_t fz =0;
896       if (ltrp->GetSide()==0){
897         fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
898       }else{
899         fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);        
900       }
901       if (npointsAC>10){
902         fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
903       }
904       fFitZ[id]=mdrift-fz;
905     }
906     if (fStreamLevel>0){
907       TTreeSRedirector *cstream = GetDebugStreamer();
908       TTimeStamp tstamp(fTime);
909       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
910       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
911       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
912       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
913       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
914       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
915       TVectorD vecGoofie(20);
916       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
917       if (goofieArray) 
918         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
919           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
920           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
921         }
922
923       if (cstream){
924         (*cstream)<<"driftv"<<
925           "run="<<fRun<<              //  run number
926           "event="<<fEvent<<          //  event number
927           "time="<<fTime<<            //  time stamp of event
928           "trigger="<<fTrigger<<      //  trigger
929           "mag="<<fMagF<<             //  magnetic field
930           // Environment values
931           "press0="<<valuePressure0<<
932           "press1="<<valuePressure1<<
933           "pt0="<<ptrelative0<<
934           "pt1="<<ptrelative1<<
935           "temp0="<<temp0<<
936           "temp1="<<temp1<<
937           "vecGoofie.="<<&vecGoofie<<
938           //
939           //
940           "iter="<<iter<<
941           "driftA.="<<fFitAside<<
942           "driftC.="<<fFitCside<<
943           "driftAC.="<<fFitACside<<
944           "chi2A="<<chi2A<<
945           "chi2C="<<chi2C<<
946           "chi2AC="<<chi2AC<<
947           "nA="<<npointsA<<
948           "nC="<<npointsC<<
949           "nAC="<<npointsAC<<
950           "\n";
951       }
952     }
953   }
954 }
955 Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
956   //
957   // Fit corrections to the drift velocity - linear approximation in the z and global y
958   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
959   // 
960   // Source of outlyers : 
961   // 0. Track in the saturation - postpeak
962   // 1. gating grid close the part of the signal for first bundle
963
964   // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
965   // 1. Robust fit is used in the itteration number 0
966   // only fraction of laser used
967   // 2. Only the tracks close to the fit used in the second itteration
968   /*
969     Formulas:
970     
971     z  = s* (z0 - vd*(t-t0))
972     
973     s  - side -1 and +1 
974     t0 - time 0
975     vd - nominal drift velocity
976     zs - miscalibrated position
977     
978     zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
979     vr  - relative change of the drift velocity
980     dzt - vd*dt
981     dr  = zz0-s*z
982     ..
983     ==>
984     zs ~ z - s*vr*(z0-s*z)+s*dzt
985     --------------------------------
986     1. Correction function vr constant:
987     
988     
989     dz = zs-z = -s*vr *(z0-s*z)+s*dzt         
990     dzs/dl = dz/dl +s*s*vr*dz/dl 
991     d(dz/dl) = vr*dz/dl     
992   */
993   const Int_t   knLaser      = 336;     //n laser tracks
994   const Float_t kFraction[2] = {0.70,0.95};    // robust fit fraction
995
996   const Float_t kSaturCut    = 0.05;    // remove saturated lasers - cut on fraction of saturated 
997   const Float_t kDistCut     = 3.;      // distance sigma cut - 3 sigma
998   const Float_t kDistCutAbs  = 1.;      // absolute cut 1 cm
999   const Float_t kMinClusters = 40.;      // minimal amount of the clusters
1000   const Float_t kMinSignal   = 2.5;      // minimal mean height of the signal
1001   const Float_t kChi2Cut     = 1.0;     // chi2 cut to accept drift fit
1002   //
1003   static TLinearFitter fdriftA(3,"hyp2");
1004   static TLinearFitter fdriftC(3,"hyp2");
1005   static TLinearFitter fdriftAC(4,"hyp3");
1006   TVectorD fitA(3),fitC(3),fitAC(4);
1007   
1008   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
1009   AliTPCParam  * tpcparam    = calib->GetParameters();    
1010   //
1011   // reset old data
1012   //
1013   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
1014   if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1015   if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1016   for (Int_t i=0;i<5; i++){
1017     (*fFitCside)[i]=0;
1018     (*fFitAside)[i]=0;
1019   }
1020   //
1021   //
1022   Float_t chi2A = 10;
1023   Float_t chi2C = 10;
1024   Float_t chi2AC = 10;
1025   Int_t npointsA=0;
1026   Int_t npointsC=0;
1027   Int_t npointsAC=0;
1028   Int_t nbA[4]={0,0,0,0};
1029   Int_t nbC[4]={0,0,0,0};
1030   TVectorD vecZM(336);      // measured z potion of laser
1031   TVectorD vecA(336);       // accepted laser
1032   TVectorD vecZF(336);      // fitted position  
1033   TVectorD vecDz(336);      // deltaZ
1034   TVectorD vecZS(336);      // surveyed position of laser
1035   // additional variable to cut
1036   TVectorD vecdEdx(336);    // dEdx
1037   TVectorD vecSy(336);      // shape y
1038   TVectorD vecSz(336);      // shape z
1039   //
1040   //
1041   for (Int_t id=0; id<336; id++){
1042     Int_t reject=0;
1043     AliTPCLaserTrack *ltrp =
1044       (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1045     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1046     AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1047     vecZM(id)= (param==0) ? 0:param->GetZ();
1048     vecZS(id)= ltrp->GetZ();
1049     vecDz(id)= 0;
1050     vecA(id)=1;
1051     vecdEdx(id)=(seed)?seed->GetdEdx():0;
1052     vecSy(id) =(seed)?seed->CookShape(1):0;
1053     vecSz(id) =(seed)?seed->CookShape(2):0;
1054     //
1055     fFitZ[id]=0;
1056     if (!fTracksEsdParam.At(id)) reject|=1;
1057     if (!param) continue; 
1058     if (!AcceptLaser(id))        reject|=2;
1059     if ( fClusterSatur[id]>kSaturCut)  reject|=4;
1060     if ( fClusterCounter[id]<kMinClusters)  reject|=8;
1061     vecA(id)=reject;
1062     if (reject>0) continue;
1063     if (ltrp->GetSide()==0){
1064       npointsA++;
1065       nbA[ltrp->GetBundle()]++;
1066     }
1067     if (ltrp->GetSide()==1){
1068       npointsC++;
1069       nbC[ltrp->GetBundle()]++;
1070     }
1071   }
1072   //
1073   // reject  "bad events"
1074   //
1075   Bool_t isOK=kTRUE;
1076   Int_t  naA=0;
1077   Int_t  naC=0;
1078   if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1079     isOK=kFALSE;
1080   for (Int_t i=0;i<4;i++){
1081     //count accepted for all layers
1082     if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1083     if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1084   }
1085   if (naA<3 &&naC<3) isOK=kFALSE;
1086   if (isOK==kFALSE) return kFALSE;
1087   
1088   //
1089   //
1090   //
1091   for (Int_t iter=0; iter<2; iter++){
1092     fdriftA.ClearPoints();
1093     fdriftC.ClearPoints();
1094     fdriftAC.ClearPoints(); 
1095     npointsA=0;  npointsC=0;  npointsAC=0;
1096     //
1097     for (Int_t id=0; id<336; id++){
1098       if (!fTracksEsdParam.At(id)) continue;
1099       if (!AcceptLaser(id)) continue;
1100       if ( fClusterSatur[id]>kSaturCut)  continue;
1101       if ( fClusterCounter[id]<kMinClusters)  continue;
1102       AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1103       if (track->GetTPCsignal()<kMinSignal) continue;
1104       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1105       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1106       Double_t xyz[3];
1107       Double_t pxyz[3];
1108       Double_t lxyz[3];
1109       Double_t lpxyz[3];
1110       param->GetXYZ(xyz);
1111       param->GetPxPyPz(pxyz);
1112       ltrp->GetXYZ(lxyz);
1113       ltrp->GetPxPyPz(lpxyz);
1114       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1115       //if (npointsAC>0) sz =TMath::Sqrt(chi2AC); 
1116       if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1117       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1118
1119 //       // drift distance
1120 //       Double_t zlength =  tpcparam->GetZLength(0);
1121 //       Double_t ldrift  = zlength-TMath::Abs(lxyz[2]);
1122 //       Double_t mdrift  = zlength-TMath::Abs(xyz[2]);
1123       //
1124       Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1125       Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1126       Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1127
1128
1129       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1130       if (iter==0 &&ltrp->GetBundle()==0) continue;  
1131       // skip bundle 0 - can be wrong in the 0 iteration
1132       if (ltrp->GetSide()==0){
1133         fdriftA.AddPoint(xxx,mdrift,1);
1134       }else{
1135         fdriftC.AddPoint(xxx,mdrift,1);
1136       }
1137       Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1138       fdriftAC.AddPoint(xxx2,mdrift,1);
1139     }
1140     //
1141     if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1142       //
1143       fdriftA.Eval();
1144       //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
1145       npointsA= fdriftA.GetNpoints();
1146       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1147       fdriftA.EvalRobust(kFraction[iter]);
1148       fdriftA.GetParameters(fitA);
1149       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1150         if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1151         (*fFitAside)[0] = fitA[0];
1152         (*fFitAside)[1] = fitA[1];
1153         (*fFitAside)[2] = fitA[2];
1154         (*fFitAside)[3] = fdriftA.GetNpoints();
1155         (*fFitAside)[4] = chi2A;        
1156       }
1157     }
1158     if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1159       fdriftC.Eval();
1160       //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
1161       npointsC= fdriftC.GetNpoints();
1162       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1163       fdriftC.EvalRobust(kFraction[iter]);
1164       fdriftC.GetParameters(fitC);
1165       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) { 
1166         if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1167         (*fFitCside)[0] = fitC[0];
1168         (*fFitCside)[1] = fitC[1];
1169         (*fFitCside)[2] = fitC[2];
1170         (*fFitCside)[3] = fdriftC.GetNpoints();
1171         (*fFitCside)[4] = chi2C;        
1172       }
1173     }
1174
1175     if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1176       fdriftAC.Eval();
1177       //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
1178       npointsAC= fdriftAC.GetNpoints();
1179       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1180       fdriftAC.EvalRobust(kFraction[iter]);
1181       fdriftAC.GetParameters(fitAC);
1182       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1183       (*fFitACside)[0] = fitAC[0];
1184       (*fFitACside)[1] = fitAC[1];
1185       (*fFitACside)[2] = fitAC[2];
1186       (*fFitACside)[3] = fdriftAC.GetNpoints();
1187       (*fFitACside)[4] = chi2AC;
1188     }
1189     
1190     for (Int_t id=0; id<336; id++){
1191       if (!fTracksEsdParam.At(id)) continue;
1192       //
1193       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1194       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1195       Double_t xyz[3];
1196       Double_t pxyz[3];
1197       Double_t lxyz[3];
1198       Double_t lpxyz[3];
1199       param->GetXYZ(xyz);
1200       param->GetPxPyPz(pxyz);
1201       ltrp->GetXYZ(lxyz);
1202       ltrp->GetPxPyPz(lpxyz); 
1203       //Double_t zlength =  tpcparam->GetZLength(0);
1204       //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1205       //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1206       Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1207       Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1208       Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1209
1210
1211       Float_t fz =0;
1212       if (ltrp->GetSide()==0){
1213         fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1214       }else{
1215         fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);        
1216       }
1217       if (npointsAC>10){
1218         //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1219       }
1220       fFitZ[id]=mdrift-fz;
1221       vecZF[id]=fz;
1222       vecDz[id]=mdrift-fz;
1223     }
1224     if (fStreamLevel>0){
1225       TTreeSRedirector *cstream = GetDebugStreamer();
1226       TTimeStamp tstamp(fTime);
1227       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1228       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1229       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1230       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1231       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1232       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1233       TVectorD vecGoofie(20);
1234       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1235       if (goofieArray) 
1236         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1237           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1238           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1239         }
1240
1241       if (cstream){
1242         (*cstream)<<"driftvN"<<
1243           "run="<<fRun<<              //  run number
1244           "event="<<fEvent<<          //  event number
1245           "time="<<fTime<<            //  time stamp of event
1246           "trigger="<<fTrigger<<      //  trigger
1247           "mag="<<fMagF<<             //  magnetic field
1248           // Environment values
1249           "press0="<<valuePressure0<<
1250           "press1="<<valuePressure1<<
1251           "pt0="<<ptrelative0<<
1252           "pt1="<<ptrelative1<<
1253           "temp0="<<temp0<<
1254           "temp1="<<temp1<<
1255           "vecGoofie.="<<&vecGoofie<<
1256           //
1257           //
1258           "vecZM.="<<&vecZM<<   // measured z position
1259           "vecZS.="<<&vecZS<<   // surveyed z position
1260           "vecZF.="<<&vecZF<<   // fitted   z position
1261           "vecDz.="<<&vecDz<<   // fitted   z position
1262           "vecA.="<<&vecA<<     // accept laser flag
1263           "vecdEdx.="<<&vecdEdx<<   // dEdx  - to cut on
1264           "vecSy.="<<&vecSy<<       // shape y - to cut on
1265           "vecSz.="<<&vecSz<<       // shape z - to cut on
1266           //
1267           "iter="<<iter<<
1268           "driftA.="<<fFitAside<<
1269           "driftC.="<<fFitCside<<
1270           "driftAC.="<<fFitACside<<
1271           "chi2A="<<chi2A<<
1272           "chi2C="<<chi2C<<
1273           "chi2AC="<<chi2AC<<
1274           "nA="<<npointsA<<
1275           "nC="<<npointsC<<
1276           "nAC="<<npointsAC<<
1277           "\n";
1278         /*
1279           //
1280           variables to check in debug mode:
1281           //
1282           chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1283           chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1284           chainDriftN->SetAlias("driftF","vecZF.fElements");      
1285           chainDriftN->SetAlias("deltaZ","driftF-driftM");  //deltaZ
1286           TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1287           TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1288           
1289          */
1290       }
1291     }
1292   }
1293   return kTRUE;
1294 }
1295
1296 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1297   //
1298   // get distance between mirror and laser track
1299   //
1300   //
1301   Double_t xyz[3];
1302   Double_t lxyz[3];
1303   param->GetXYZ(xyz);
1304   ltrp->GetXYZ(lxyz);
1305   //
1306   //
1307   Double_t dist = 0;
1308   //radial distance
1309   dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1310   //
1311   // z distance
1312   // apply drift correction if already exist
1313   //
1314   Float_t dz = 0;
1315   if (ltrp->GetSide()==0){
1316     if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1317   }else{
1318     if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1319   }
1320   if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1321   dist+=TMath::Abs(dz);
1322   //
1323   // phi dist - divergence on 50 cm
1324   //
1325   dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1326   return dist;
1327 }
1328
1329
1330 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
1331   //
1332   //
1333   //
1334   /*
1335   TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1336   TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1337   TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1338   TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1339
1340   TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1341   */
1342   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1343   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1344   //AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1345   Double_t xyz[3];
1346   Double_t lxyz[3];
1347   param->GetXYZ(xyz);
1348   ltrp->GetXYZ(lxyz);
1349   if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1350   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE;    // cutZ -P1
1351   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE;  // cut -P2
1352   if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE;   // cut Tl -P3
1353   //
1354   //
1355
1356   return kTRUE;
1357 }
1358
1359 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1360   //
1361   // Find corresponding mirror
1362   // add the corresponding tracks
1363  
1364   if (!track->GetOuterParam()) return -1;
1365
1366   Float_t kRadius0  = 252;
1367   Float_t kRadius   = 254.2;
1368   Int_t countercl=0;
1369   Float_t counterSatur=0;
1370   Int_t csideA =0;
1371   Int_t csideC =0;
1372   for (Int_t irow=158;irow>-1;--irow) {
1373     AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1374     if (!c) continue;
1375     Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1376     Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1377     if (pedgeY<3) continue;
1378     if (pedgeX<3) continue;
1379     countercl++;
1380     if (c->GetDetector()%36<18) csideA++;
1381     if (c->GetDetector()%36>=18) csideC++;
1382     if (c->GetMax()>900) counterSatur++;
1383   }
1384   counterSatur/=(countercl+1);
1385   //
1386   //
1387   //
1388   if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0;  // cross laser track can not happen
1389
1390   Int_t side= 0;
1391   if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1392
1393
1394   AliExternalTrackParam param(*(track->GetOuterParam()));
1395   AliTracker::PropagateTrackTo(&param,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
1396   AliTracker::PropagateTrackTo(&param,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
1397   AliTPCLaserTrack ltr;
1398   AliTPCLaserTrack *ltrp=0x0;
1399   //  AliTPCLaserTrack *ltrpjw=0x0;
1400   //
1401   Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param,side);
1402  // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
1403   //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1404
1405   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1406     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1407   else
1408     ltrp=&ltr;
1409
1410   if (id<0) return -1;
1411   if (ltrp->GetSide()!=side) return -1;
1412   fCounter[id]++;
1413   //
1414   //
1415   //
1416   //
1417   if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1418   //
1419   //
1420   Float_t radius=TMath::Abs(ltrp->GetX());
1421   param.Rotate(ltrp->GetAlpha());
1422   AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kFALSE);
1423   //
1424   if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1425   Bool_t accept=kTRUE;  
1426   //
1427   // choose closer track
1428   //
1429   AliExternalTrackParam * param0  = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1430   if (param0){
1431     Float_t dist0=GetDistance(param0,ltrp);
1432     Float_t dist1=GetDistance(&param,ltrp);
1433     if (dist0<dist1)    accept=kFALSE;       
1434   }
1435   
1436   if (accept){
1437     fClusterCounter[id]=countercl;
1438     fTracksEsdParam.AddAt(param.Clone(),id);
1439     fTracksEsd.AddAt(track,id);
1440     fTracksTPC.AddAt(seed,id);
1441   }
1442   return id;
1443 }
1444
1445
1446
1447 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1448   //
1449   //  Dump Laser info to the tree
1450   //
1451   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1452   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1453   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1454   //
1455   // Fast laser ID
1456   //
1457   Double_t xyz[3];
1458   Double_t pxyz[3];
1459   Double_t lxyz[3];
1460   Double_t lpxyz[3];
1461   param->GetXYZ(xyz);
1462   param->GetPxPyPz(pxyz);
1463   ltrp->GetXYZ(lxyz);
1464   ltrp->GetPxPyPz(lpxyz);
1465   Float_t dist3D   = GetDistance(param,ltrp);
1466   Float_t dist0    = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1467   Float_t distphi  = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1468   
1469   
1470   if (fStreamLevel>0){
1471     TTreeSRedirector *cstream = GetDebugStreamer();
1472     Int_t time = fESD->GetTimeStamp();
1473     Bool_t accept = AcceptLaser(id);
1474     if (cstream){
1475       (*cstream)<<"Track"<<
1476         //
1477         "run="<<fRun<<              //  run number
1478         "event="<<fEvent<<          //  event number
1479         "time="<<fTime<<            //  time stamp of event
1480         "trigger="<<fTrigger<<      //  trigger
1481         "mag="<<fMagF<<             //  magnetic field
1482
1483         "id="<<id<<
1484         "accept="<<accept<<
1485         "driftA.="<<fFitAside<<
1486         "driftC.="<<fFitCside<<
1487         "time="<<time<<
1488         "dist3D="<<dist3D<<
1489         "dist0="<<dist0<<
1490         "distphi="<<distphi<<
1491         //
1492         //
1493         "counter="<<fCounter[id]<<
1494         "clcounter="<<fClusterCounter[id]<<
1495         "satur="<<fClusterSatur[id]<<
1496         "fz="<<fFitZ[id]<<
1497         //
1498         "LTr.="<<ltrp<<
1499         "Esd.="<<track<<
1500         "Tr.="<<param<<
1501         "x0="<<xyz[0]<<
1502         "x1="<<xyz[1]<<
1503         "x2="<<xyz[2]<<
1504         "px0="<<pxyz[0]<<
1505         "px1="<<pxyz[1]<<
1506         "px2="<<pxyz[2]<<
1507         //
1508         "lx0="<<lxyz[0]<<
1509         "lx1="<<lxyz[1]<<
1510         "lx2="<<lxyz[2]<<
1511         "lpx0="<<lpxyz[0]<<
1512         "lpx1="<<lpxyz[1]<<
1513         "lpx2="<<lpxyz[2]<<
1514         "\n";
1515     }
1516   }
1517 }
1518
1519 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1520   //
1521   // Refit the track with different tracklet models:
1522   // 1. Per ROC using the kalman filter, different edge cuts
1523   // 2. Per ROC linear in y and z
1524   // 3. Per ROC quadratic in y and z
1525   // 4. Per track offset for each sector, linear for each sector, common quadratic
1526   // store x, y, z information for all models and the cluster to calculate the residuals
1527   //
1528   
1529   // The clusters which do not fulfill given  criteria are skipped
1530   //
1531   // Cluters removed from fit
1532   //  1. Extended shape     > kShapeCut
1533   //  2. In saturation  Max > kMax
1534   //  3. Distance to edge   < cutEdge
1535   //
1536   // Clusters not used for the calibration:
1537   //
1538   //  1. Extended shape     > kShapeCut
1539   //  2. In saturation  Max > kMax
1540
1541
1542   AliTPCseed *track      = (AliTPCseed*)fTracksTPC.At(id);
1543   AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1544   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1545
1546   AliTPCclusterMI dummyCl;
1547
1548   //two tracklets
1549   Int_t kMaxTracklets=2;  
1550   Float_t kShapeCut = 1.3;
1551   Float_t kRatioCut = 2.;
1552
1553   Float_t kMax      = 900.;
1554
1555
1556   //=============================================//
1557   // Linear Fitters for the Different Approaches //
1558   //=============================================//
1559   //linear fit model in y and z; inner - outer sector, combined with offset
1560   static TLinearFitter fy1I(2,"hyp1");
1561   static TLinearFitter fy1O(2,"hyp1");
1562   static TLinearFitter fz1I(2,"hyp1");
1563   static TLinearFitter fz1O(2,"hyp1");
1564   static TLinearFitter fy1IO(3,"hyp2");
1565   static TLinearFitter fz1IO(3,"hyp2");
1566   //quadratic fit model in y and z; inner - sector
1567   static TLinearFitter fy2I(3,"hyp2");
1568   static TLinearFitter fy2O(3,"hyp2");
1569   static TLinearFitter fz2I(3,"hyp2");
1570   static TLinearFitter fz2O(3,"hyp2");
1571   //common quadratic fit for IROC and OROC in y and z
1572   static TLinearFitter fy4(5,"hyp4");
1573   static TLinearFitter fz4(5,"hyp4");
1574
1575
1576   //set standard cuts
1577   if ( fNcuts==0 ){
1578       fNcuts=1;
1579       fEdgeXcuts[0]=4;
1580       fEdgeYcuts[0]=3;
1581       fNClCuts[0]=20;
1582   }
1583   //=============================//
1584   // Loop over all Tracklet Cuts //
1585   //=============================//
1586   for (Int_t icut=0; icut<fNcuts; icut++){
1587     Float_t xinMin = 2500, xinMax=-90;
1588     Float_t xoutMin=2500, xoutMax=-90;
1589     Float_t msigmaYIn=0;
1590     Float_t msigmaYOut=0;
1591     Float_t mqratioIn=0;
1592     Float_t mqratioOut=0;
1593
1594
1595       AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1596       //cut parameters
1597       Double_t edgeCutX = fEdgeXcuts[icut];
1598       Double_t edgeCutY = fEdgeYcuts[icut];
1599       Int_t    nclCut   = (Int_t)fNClCuts[icut];
1600       //=========================//
1601       // Parameters to calculate //
1602       //=========================//
1603       //fit parameter inner, outer tracklet and combined fit
1604       TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1605       TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1606       //
1607       TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1608       TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1609       TVectorD vecy4res(5),vecz4res(5);
1610       TVectorD vecy1resIO(3),vecz1resIO(3);
1611       // cluster and track positions for each row - used for residuals
1612       TVectorD vecgX(159);        // global X
1613       TVectorD vecgY(159);        // global Y
1614       TVectorD vecgZ(159);        // global Z
1615
1616       TVectorD vecX(159);        // x is the same for all (row center)
1617       TVectorD vecYkalman(159);  // y from kalman fit
1618       TVectorD vecZkalman(159);  // z from kalman fit
1619       TVectorD vecY1(159);       // y from pol1 fit per ROC
1620       TVectorD vecZ1(159);       // z from pol1 fit per ROC
1621       TVectorD vecY1IO(159);     // y from pol1 fit per ROC
1622       TVectorD vecZ1IO(159);     // z from pol1 fit per ROC
1623       TVectorD vecY2(159);       // y from pol2 fit per ROC
1624       TVectorD vecZ2(159);       // z from pol2 fit per ROC
1625       TVectorD vecY4(159);       // y from sector fit
1626       TVectorD vecZ4(159);       // z from sector fit
1627       TVectorD vecClY(159);      // y cluster position
1628       TVectorD vecClZ(159);      // z cluster position
1629       TVectorD vecSec(159);      // sector for each row
1630       TVectorD  isReject(159);    // flag - cluster to be rejected 
1631       //chi2 of fits
1632       Double_t chi2I1z=0;       // chi2 of pol1 fit in z (inner)
1633       Double_t chi2I1y=0;       // chi2 of pol1 fit in y (inner)
1634       Double_t chi2O1z=0;       // chi2 of pol1 fit in z (outer)
1635       Double_t chi2O1y=0;       // chi2 of pol1 fit in y (outer)
1636       Double_t chi2IO1z=0;       // chi2 of pol1 fit in z (outer)
1637       Double_t chi2IO1y=0;       // chi2 of pol1 fit in y (outer)
1638       Double_t chi2I2z=0;       // chi2 of pol2 fit in z (inner)
1639       Double_t chi2I2y=0;       // chi2 of pol2 fit in y (inner)
1640       Double_t chi2O2z=0;       // chi2 of pol2 fit in z (outer)
1641       Double_t chi2O2y=0;       // chi2 of pol2 fit in y (outer)
1642       Double_t chi2IOz=0;       // chi2 of hyp4 fit in z (inner+outer)
1643       Double_t chi2IOy=0;       // chi2 of hyp4 fit in y (inner+outer)
1644       //more
1645       Int_t innerSector = -1;    // number of inner sector
1646       Int_t outerSector = -1;    // number of outer sector
1647       Int_t nclI=0;              // number of clusters (inner)
1648       Int_t nclO=0;              // number of clusters (outer)
1649       //=================================================//
1650       // Perform the Kalman Fit using the Tracklet Class //
1651       //=================================================//
1652       AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1653       TObjArray tracklets=
1654           AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1655                                           kFALSE,nclCut,kMaxTracklets);
1656       // tracklet pointers
1657       AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1658       AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1659       AliTPCTracklet *tr=0x0;
1660       AliTPCTracklet dummy;
1661       //continue if we didn't find a tracklet
1662       if ( !trInner && !trOuter ) continue;
1663       //================================================================================//
1664       // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1665       //================================================================================//
1666       if ( trInner && trOuter ){
1667           if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1668           if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1669               tr = trInner;
1670               trInner=trOuter;
1671               trOuter=tr;
1672               AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1673           }
1674       } else {
1675           if ( trInner ){
1676               if ( !trInner->GetInner() ) continue;
1677               trOuter=&dummy;
1678               if ( trInner->GetSector()>35 ){
1679                   trOuter=trInner;
1680                   trInner=&dummy;
1681               }
1682           } else { //trOuter
1683               if ( !trOuter->GetInner() ) continue;
1684               trInner=&dummy;
1685               if ( trOuter->GetSector()<36 ){
1686                   trInner=trOuter;
1687                   trOuter=&dummy;
1688               }
1689           }
1690       }
1691       innerSector = trInner->GetSector();
1692       if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1693       outerSector = trOuter->GetSector();
1694       if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1695
1696       // array of clusters
1697       TClonesArray arrCl("AliTPCclusterMI",159);
1698       arrCl.ExpandCreateFast(159);
1699       //=======================================//
1700       // fill fitters with cluster information //
1701       //=======================================//
1702       AliDebug(3,"Filing Cluster Information");
1703
1704       //
1705       //
1706
1707       for (Int_t irow=158;irow>-1;--irow) {
1708           AliTPCclusterMI *c=track->GetClusterPointer(irow);
1709           AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1710           cl=dummyCl;
1711           vecX[irow]   = 0;
1712           vecClY[irow] = 0;
1713           vecClZ[irow] = 0;
1714           Float_t meanY=0, sumY=0;
1715           for (Int_t drow=-1;drow<=1;drow++) {
1716             if (irow+drow<0) continue;
1717             if (irow+drow>158) continue;
1718             AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1719             if (!ccurrent) continue;
1720             Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1721             if ( roc!=innerSector && roc!=outerSector ) continue;
1722             if (ccurrent->GetX()<10) continue;
1723             if (ccurrent->GetY()==0) continue;
1724             meanY+=ccurrent->GetY();
1725             sumY++;
1726           }
1727           if (sumY>0)  meanY/=sumY; 
1728          
1729           //
1730           vecSec[irow]=-1;
1731           if (!c) continue;
1732           Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1733           Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1734           
1735           //
1736           Int_t roc = static_cast<Int_t>(c->GetDetector());
1737           if ( roc!=innerSector && roc!=outerSector ) continue;
1738           vecSec[irow]=roc;
1739           //store clusters in clones array
1740           cl=*c;          
1741           //
1742           if (c->GetMax()<4)   continue;  // noise cluster?
1743           if (TMath::Abs(c->GetY())<0.0001)   continue;  // noise cluster?
1744           //cluster position
1745           vecX[irow]   = c->GetX();
1746           vecClY[irow] = c->GetY();
1747           vecClZ[irow] = c->GetZ();
1748           //
1749 //        Float_t gxyz[3];
1750 //        c->GetGlobalXYZ(gxyz);
1751 //        vecgX[irow]   = gxyz[0];
1752 //        vecgY[irow]   = gxyz[1];
1753 //        vecgZ[irow]   = gxyz[2];
1754           //
1755           Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1756           Double_t y = vecClY[irow];
1757           Double_t z = vecClZ[irow];
1758           //
1759           Double_t x2[2]={x,x*x};   //linear and parabolic parameters
1760           Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1761           Double_t xIO[2]={0,x};    //common linear + offset IROC-OROC
1762           if ( roc == innerSector ) {
1763               x4[0]=1; //offset inner - outer sector
1764               x4[1]=x; //slope parameter inner sector
1765               xIO[0]=1;
1766           } else {
1767               x4[2]=x; //slope parameter outer sector
1768           }
1769           x4[3]=x*x;   //common parabolic shape
1770           if (pedgeX < fEdgeXcuts[icut]) continue;
1771           if (pedgeY < fEdgeYcuts[icut]) continue;
1772           if (c->GetMax()>900) continue;  // cluster in saturation        
1773           //
1774           if ( roc==innerSector ){
1775               fy1I.AddPoint(x2,y);
1776               fz1I.AddPoint(x2,z);
1777               fy2I.AddPoint(x2,y);
1778               fz2I.AddPoint(x2,z);
1779               ++nclI;
1780               if (c->GetX()<xinMin) xinMin=c->GetX();
1781               if (c->GetX()>xinMax) xinMax=c->GetX();
1782
1783               msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1784               mqratioIn +=c->GetMax()/c->GetQ();
1785
1786           }
1787           if ( roc==outerSector ){
1788               fy1O.AddPoint(x2,y);
1789               fz1O.AddPoint(x2,z);
1790               fy2O.AddPoint(x2,y);
1791               fz2O.AddPoint(x2,z);
1792               ++nclO;
1793               if (c->GetX()<xoutMin) xoutMin=c->GetX();
1794               if (c->GetX()>xoutMax) xoutMax=c->GetX();
1795               msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1796               mqratioOut +=c->GetMax()/c->GetQ();
1797
1798           }
1799           fy4.AddPoint(x4,y);
1800           fz4.AddPoint(x4,z);
1801           fy1IO.AddPoint(xIO,y);
1802           fz1IO.AddPoint(xIO,z);
1803       }
1804       if (nclI>0)  {
1805         msigmaYIn/=nclI;
1806         mqratioIn/=nclI;
1807       }
1808       if (nclO>0)  {
1809         msigmaYOut/=nclO;
1810         mqratioOut/=nclO;
1811       }
1812       //======================================//
1813       // Evaluate and retrieve fit parameters //
1814       //======================================//
1815       AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1816       //inner sector
1817       if (  (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1818           fy1I.Eval();
1819           fz1I.Eval();
1820           fy2I.Eval();
1821           fz2I.Eval();
1822           fy1I.GetParameters(vecy1resInner);
1823           fz1I.GetParameters(vecz1resInner);
1824           fy2I.GetParameters(vecy2resInner);
1825           fz2I.GetParameters(vecz2resInner);
1826           chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1827           chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1828           chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1829           chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1830       }
1831       //outer sector
1832       if (  (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1833           fy1O.Eval();
1834           fz1O.Eval();
1835           fy2O.Eval();
1836           fz2O.Eval();
1837           fy1O.GetParameters(vecy1resOuter);
1838           fz1O.GetParameters(vecz1resOuter);
1839           fy2O.GetParameters(vecy2resOuter);
1840           fz2O.GetParameters(vecz2resOuter);
1841           chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1842           chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1843           chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1844           chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1845       }
1846       //combined hyp4 fit
1847       if ( innerSector>0 && outerSector>0 ){
1848           if (fy4.GetNpoints()>0) {
1849               fy4.Eval();
1850               fy4.GetParameters(vecy4res);
1851               chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1852           }
1853           if (fz4.GetNpoints()>0) {
1854               fz4.Eval();
1855               fz4.GetParameters(vecz4res);
1856               chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1857           }
1858           if (fy1IO.GetNpoints()>0) {
1859             fy1IO.Eval();
1860             fy1IO.GetParameters(vecy1resIO);
1861             chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1862           }
1863           if (fz1IO.GetNpoints()>0) {
1864             fz1IO.Eval();
1865             fz1IO.GetParameters(vecz1resIO);
1866             chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1867           }
1868       }
1869       //clear points
1870       fy4.ClearPoints();  fz4.ClearPoints();
1871       fy1I.ClearPoints(); fy1O.ClearPoints();
1872       fz1I.ClearPoints(); fz1O.ClearPoints();
1873       fy2I.ClearPoints(); fy2O.ClearPoints();
1874       fz2I.ClearPoints(); fz2O.ClearPoints();
1875       fy1IO.ClearPoints(); fz1IO.ClearPoints();
1876       //==============================//
1877       // calculate tracklet positions //
1878       //==============================//
1879       AliDebug(4,"Calculate tracklet positions");
1880       for (Int_t irow=158;irow>-1;--irow) {
1881         isReject[irow]=0;
1882          AliTPCclusterMI *c=track->GetClusterPointer(irow);
1883          if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1884            isReject[irow]+=1;
1885          }
1886
1887          if (!c) { //no cluster    
1888            isReject[irow]+=2;
1889          }else{
1890            if (c->GetMax()>kMax){   //saturation
1891              isReject[irow]+=4;
1892            }
1893            if ( vecSec[irow] == outerSector ) {   //extended shape
1894              if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut)                isReject[irow]+=8;
1895              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut)    isReject[irow]+=16;
1896            }else{
1897              if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut)                isReject[irow]+=8;
1898              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut)    isReject[irow]+=16;
1899            }       
1900          }
1901
1902
1903
1904           if ( vecSec[irow]==-1 ) continue;  //no cluster info
1905           if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1906           tr=&dummy;
1907           Double_t x=vecX[irow];
1908           Double_t xref=x-133.4;
1909           //
1910           Double_t yoffInner=0;
1911           Double_t zoffInner=0;
1912           Double_t yoffInner1=0;
1913           Double_t zoffInner1=0;
1914           Double_t yslopeInner=0;
1915           Double_t yslopeOuter=0;
1916           Double_t zslopeInner=0;
1917           Double_t zslopeOuter=0;
1918           //positions of hyperplane fits
1919           if ( vecSec[irow] == outerSector ) {
1920               tr=trOuter;
1921               vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1922               vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1923               vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1924               vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1925               yslopeOuter=vecy4res[3];
1926               zslopeOuter=vecz4res[3];
1927           } else {
1928               tr=trInner;
1929               vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1930               vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1931               vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1932               vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1933               yoffInner=vecy4res[1];
1934               zoffInner=vecz4res[1];
1935               yoffInner1=vecy1resIO[1];
1936               zoffInner1=vecz1resIO[1];
1937               yslopeInner=vecy4res[2];
1938               zslopeInner=vecz4res[2];
1939           }
1940           vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1941           vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1942           vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1943           vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1944           //positions of kalman fits
1945           Double_t gxyz[3],xyz[3];
1946           AliExternalTrackParam *param = 0x0;
1947           //
1948           param=tr->GetInner();
1949           if (param){
1950               param->GetXYZ(gxyz);
1951               Float_t bz = AliTracker::GetBz(gxyz);
1952               param->GetYAt(x, bz, xyz[1]);
1953               param->GetZAt(x, bz, xyz[2]);
1954               vecYkalman[irow]=xyz[1];
1955               vecZkalman[irow]=xyz[2];
1956           }
1957           //
1958           //
1959           //
1960          
1961       }
1962       //=====================================================================//
1963       // write results from the different tracklet fits with debug streamers //
1964       //=====================================================================//
1965       if (fStreamLevel>4){
1966           TTreeSRedirector *cstream = GetDebugStreamer();
1967           if (cstream){
1968             Float_t dedx = track->GetdEdx();
1969             (*cstream)<<"FitModels"<<
1970               "run="<<fRun<<              //  run number
1971               "event="<<fEvent<<          //  event number
1972               "time="<<fTime<<            //  time stamp of event
1973               "trigger="<<fTrigger<<      //  trigger
1974               "mag="<<fMagF<<             //  magnetic field          
1975               //
1976               "cutNr="      << icut <<
1977               "edgeCutX="   << edgeCutX <<
1978               "edgeCutY="   << edgeCutY <<
1979               "nclCut="     << nclCut <<
1980               "innerSector="<< innerSector <<
1981               "outerSector="<< outerSector <<
1982               "dEdx="       << dedx <<
1983               "LTr.="       << ltrp <<
1984               "Tr.="        << extparam <<
1985               "yPol1In.="   << &vecy1resInner <<
1986               "zPol1In.="   << &vecz1resInner <<
1987               "yPol1InOut.="<< &vecy1resIO <<
1988               "zPol1InOut.="<< &vecz1resIO <<
1989               "yPol2In.="   << &vecy2resInner <<
1990               "zPol2In.="   << &vecz2resInner <<
1991               "yPol1Out.="  << &vecy1resOuter <<
1992               "zPol1Out.="  << &vecz1resOuter <<
1993               "yPol2Out.="  << &vecy2resOuter <<
1994               "zPol2Out.="  << &vecz2resOuter <<
1995               "yInOut.="    << &vecy4res <<
1996               "zInOut.="    << &vecz4res <<
1997               "chi2y1In="   << chi2I1y <<
1998               "chi2z1In="   << chi2I1z <<
1999               "chi2y1InOut="<< chi2IO1y <<
2000               "chi2z1InOut="<< chi2IO1z <<
2001               "chi2y1Out="  << chi2O1y <<
2002               "chi2z1Out="  << chi2O1z <<
2003               "chi2y2In="   << chi2I2y <<
2004               "chi2z2In="   << chi2I2z <<
2005               "chi2y2Out="  << chi2O2y <<
2006               "chi2z2Out="  << chi2O2z <<
2007               "chi2yInOut=" << chi2IOy <<
2008               "chi2zInOut=" << chi2IOz <<
2009               "trletIn.="   << trInner <<
2010               "trletOut.="  << trOuter <<
2011               "nclI="       << nclI <<
2012               "nclO="       << nclO <<
2013               "xinMin="     << xinMin<<
2014               "xinMax="     << xinMax<<
2015               "xoutMin="    << xoutMin<<
2016               "xoutMax="     << xoutMax<<
2017               "msigmaYIn=" <<msigmaYIn<<
2018               "msigmaYOut=" <<msigmaYOut<<
2019               "mqratioIn=" <<mqratioIn<<
2020               "mqratioOut=" <<  mqratioOut      <<
2021               "\n";
2022           }
2023       }
2024
2025       // wirte residuals information
2026       if (fStreamLevel>5){
2027           TTreeSRedirector *cstream = GetDebugStreamer();
2028           if (cstream){
2029             Float_t dedx = track->GetdEdx();
2030             (*cstream)<<"Residuals"<<
2031               "run="<<fRun<<              //  run number
2032               "event="<<fEvent<<          //  event number
2033               "time="<<fTime<<            //  time stamp of event
2034               "trigger="<<fTrigger<<      //  trigger
2035               "mag="<<fMagF<<             //  magnetic field          
2036               //
2037               "cutNr="      << icut <<
2038               "edgeCutX="   << edgeCutX <<
2039               "edgeCutY="   << edgeCutY <<
2040               "nclCut="     << nclCut   <<
2041               "LTr.="       << ltrp <<
2042               "Tr.="        << extparam<<
2043               "dEdx="       << dedx <<
2044               "Cl.="        << &arrCl <<
2045               "vX.="        << &vecgX<<   // global x
2046               "vY.="        << &vecgY<<   // global y
2047               "vZ.="        << &vecgZ<<   // global z
2048               "TrX.="       << &vecX <<
2049               "TrYpol1.="   << &vecY1 <<
2050               "TrZpol1.="   << &vecZ1 <<
2051               "TrYpol2.="   << &vecY2 <<
2052               "TrZpol2.="   << &vecZ2 <<
2053               "TrYpol1InOut.="<< &vecY1IO <<
2054               "TrZpol1InOut.="<< &vecZ1IO <<
2055               "TrYInOut.="  << &vecY4 <<
2056               "TrZInOut.="  << &vecZ4 <<
2057               "ClY.="       << &vecClY <<
2058               "ClZ.="       << &vecClZ <<
2059               "isReject.="  << &isReject<<
2060               "sec.="       << &vecSec <<
2061               "nclI="       << nclI <<
2062               "nclO="       << nclO <<
2063               "xinMin="     << xinMin<<
2064               "xinMax="     << xinMax<<
2065               "xoutMin="    << xoutMin<<
2066               "xoutMax="     << xoutMax<<
2067               "msigmaYIn=" <<msigmaYIn<<
2068               "msigmaYOut=" <<msigmaYOut<<
2069               "mqratioIn=" <<mqratioIn<<
2070               "mqratioOut=" <<  mqratioOut      <<
2071               "yInOut.="    << &vecy4res <<
2072               "zInOut.="    << &vecz4res <<
2073               //chi2s
2074               "chi2y1In="   << chi2I1y <<    //
2075               "chi2z1In="   << chi2I1z <<
2076               "chi2y1Out="  << chi2O1y <<
2077               "chi2z1Out="  << chi2O1z <<
2078               "chi2y1InOut="<< chi2IO1y <<
2079               "chi2z1InOut="<< chi2IO1z <<
2080               "chi2y2In="   << chi2I2y <<
2081               "chi2z2In="   << chi2I2z <<
2082               "chi2y2Out="  << chi2O2y <<
2083               "chi2z2Out="  << chi2O2z <<
2084               "chi2yInOut=" << chi2IOy <<
2085               "chi2zInOut=" << chi2IOz <<
2086               // fit parameters
2087               "yPol1In.="   << &vecy1resInner <<
2088               "zPol1In.="   << &vecz1resInner <<
2089               "yPol2In.="   << &vecy2resInner <<
2090               "zPol2In.="   << &vecz2resInner <<
2091               "yPol1Out.="  << &vecy1resOuter <<
2092               "zPol1Out.="  << &vecz1resOuter <<
2093               "yPol1InOut.="<< &vecy1resIO <<
2094               "zPol1InOut.="<< &vecz1resIO <<
2095               "yPol2Out.="  << &vecy2resOuter <<
2096               "zPol2Out.="  << &vecz2resOuter <<
2097
2098               "\n";
2099
2100           }
2101       }
2102       //==========================//
2103       // Fill Residual Histograms //
2104       //==========================//
2105       if (!fHisNclIn) MakeFitHistos(); 
2106
2107       TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2108       TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2109       TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2110       TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2111       TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2112       TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2113       //      TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2114       //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2115       //
2116       for (Int_t irow=158;irow>-1;--irow) {
2117         if (vecSec[irow]==-1)continue; // no cluster info
2118         if (isReject[irow]>0.5) continue;  // 
2119         Double_t ycl  = vecClY[irow];
2120         Double_t yfit = vecY1[irow];
2121         Double_t yfit2 = vecY2[irow];
2122         Double_t x    = vecX[irow];
2123         Double_t yabsbeam = -1000;
2124         if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2125           yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2126         else if(innerSector==fBeamSectorInner[id])
2127           yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2128
2129         //      Double_t yfit3 = vecY2[irow];
2130         Double_t zcl  = vecClZ[irow];
2131         Double_t zfit = vecZ1[irow];
2132         Double_t zfit2 = vecZ2[irow];
2133         //Double_t zfit3 = vecZ2[irow];
2134
2135         // dz abs 
2136         // The expressions for zcorrected has been obtained by
2137         // inverting the fits in the FitDriftV() method (ignoring the
2138         // global y dependence for now):
2139         // A side:
2140         // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2141         // =>
2142         // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2143         // 
2144         // C side:
2145         // 250 + zmeasured  = [0] + [1]*(250+zreal) + .... (yglobal)
2146         // =>
2147         // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2148
2149         Double_t dzabs = -1000;
2150         Double_t zcorrected = -1000;
2151         if (ltrp->GetSide()==0){
2152           if ((*fFitAside)[1]>0. || fUseFixedDriftV) { 
2153             // ignore global y dependence for now
2154             zcorrected = 0;         
2155             if(!fUseFixedDriftV) 
2156               zcorrected = (zcl + (*fFitAside)[0] - 
2157                             (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2158             else
2159               zcorrected = (zcl + fFixedFitAside0 - 
2160                             (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2161             //      zcorrected = zcl;
2162             if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2163               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2164             else if(innerSector==fBeamSectorInner[id])
2165               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2166           }
2167         } else {
2168           if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2169             
2170             if(!fUseFixedDriftV) 
2171               zcorrected = (zcl - (*fFitCside)[0] + 
2172                             (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2173             else
2174               zcorrected = (zcl - fFixedFitCside0 + 
2175                             (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2176             
2177             //      zcorrected = zcl;
2178             if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2179               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2180             else if(innerSector==fBeamSectorInner[id])
2181               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2182           }
2183         }
2184         
2185         if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2186           if (profy){       
2187               profy->Fill(irow,ycl-yfit);
2188               profy2->Fill(irow,ycl-yfit2);
2189               if(yabsbeam<-100) {
2190                 fHisYAbsErrors->Fill(id);
2191                 //              profyabs->Fill(irow,-0.99);
2192               } else
2193                 profyabs->Fill(irow,ycl-yabsbeam);
2194
2195               //              profy3->Fill(irow,ycl-yfit3);
2196           }
2197           if (profz) {
2198               profz->Fill(irow,zcl-zfit);
2199               profz2->Fill(irow,zcl-zfit2);
2200               //profz3->Fill(irow,zcl-zfit3);
2201               if(dzabs<-100) {
2202
2203                 fHisZAbsErrors->Fill(id);
2204               }else
2205                 profzabs->Fill(irow,dzabs);
2206           }
2207         }
2208       }
2209       //
2210       //
2211       // Fill laser fit histograms
2212       //
2213       Float_t dedx = track->GetdEdx();
2214       if (nclI>20&&nclO>20){
2215         fHisNclIn->Fill(id,nclI);      //->Number of clusters inner
2216         fHisNclOut->Fill(id,nclO);     //->Number of clusters outer
2217         fHisNclIO->Fill(id,nclI+nclO);      //->Number of cluster inner outer
2218         //
2219         fHisLclIn->Fill(id,xinMax-xinMin);        //->Level arm inner
2220         fHisLclOut->Fill(id,xoutMax-xoutMin);     //->Level arm outer
2221         fHisLclIO->Fill(id,xoutMax-xinMin);       //->Number of cluster inner outer
2222         //
2223         fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2224         fHisdZfit->Fill(id,fFitZ[id]);
2225         // 
2226         //
2227         fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y));      //->chi2 y inner - line
2228         fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y));     //->chi2 y inner - line
2229         fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y));      //->chi2 y inner - parabola
2230         fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y));     //->chi2 y inner - parabola
2231         fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy));      //->chi2 y IO    - common
2232         
2233         
2234         fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z));      //->chi2 z inner - line
2235         fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z));     //->chi2 z inner - line
2236         fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y));      //->chi2 z inner - parabola
2237         fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z));     //->chi2 z inner - parabola
2238         fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz));      //->chi2 z IO    - common
2239         //
2240         //
2241         fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]);     //-> delta y   P0outer-P0inner - line
2242         fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]);     //-> delta y   P0outer-P0inner - parabola
2243         fHisPy3vP0->Fill(id,vecy4res[1]);                           //-> delta y   P0outer-P0inner - common parabola
2244         //
2245         fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]);     //-> delta ky  P1outer-P1inner - line
2246         fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]);     //-> delta ky  P1outer-P1inner - parabola
2247         fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]);               //-> delta ky  P1outer-P1inner - common parabola
2248         //
2249         fHisPy3vP2IO->Fill(id,vecy4res[4]);        //-> Curv  P2outerinner - common parabola
2250         fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]);     //-> delta z   P0outer-P0inner - line
2251         fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]);     //-> delta z   P0outer-P0inner - parabola
2252         fHisPz3vP0->Fill(id,vecz4res[1]);                           //-> delta z   P0outer-P0inner - common parabola
2253         //
2254         fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]);     //-> delta kz  P1outer-P1inner - line
2255         fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]);     //-> delta kz  P1outer-P1inner - parabola
2256         fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]);               //-> delta kz  P1outer-P1inner - common parabola
2257         fHisPz3vP2IO->Fill(id,vecz4res[4]);        //-> Curv  P2outerinner - common parabola
2258       }
2259       if(nclI>20){
2260         fHisPy2vP2In->Fill(id,vecy2resInner[2]);   //-> Curv  P2inner - parabola
2261         fHisPz2vP2In->Fill(id,vecz2resInner[2]);   //-> Curv  P2inner - parabola
2262       }
2263       //
2264       if (nclO>20){
2265         fHisPz2vP2Out->Fill(id,vecz2resOuter[2]);  //-> Curv  P2outer - parabola
2266         fHisPy2vP2Out->Fill(id,vecy2resOuter[2]);  //-> Curv  P2outer - parabola
2267       }
2268       
2269   }
2270   //
2271   // Fill raw THnSparses
2272   //
2273   for (Int_t irow=0;irow<159;irow++) {
2274     AliTPCclusterMI *c=track->GetClusterPointer(irow);
2275     if (!c) continue;
2276     if (c->GetMax()>800) continue; // saturation cut
2277     //if (TMath::Sqrt(TMath::Abs(c->GetSigmaY2()))>1) continue; // saturation cut
2278     //
2279     Double_t deltaY=c->GetY()-(*ltrp->GetVecLY())[irow];
2280     Double_t deltaZ=c->GetZ()-(*ltrp->GetVecLZ())[irow];
2281     //TString axisName[6]={"Delta","bin", "rms shape", "Q", "row","trackID"}
2282     Double_t xyz[6]={0, 0, 0,TMath::Sqrt(c->GetMax()),irow,id};
2283     xyz[0]=deltaY;
2284     xyz[1]=c->GetPad();
2285     xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaY2()));
2286     fHisLaserPad->Fill(xyz);
2287     xyz[0]=deltaZ;
2288     xyz[1]=c->GetTimeBin();
2289     xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaZ2()));
2290     fHisLaserTime->Fill(xyz);
2291   }
2292 }
2293
2294
2295
2296 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2297   //
2298   //  Dump information about laser beams
2299   //  isOK variable indicates usability of the beam  
2300   //  Beam is not usable if:
2301   //     a.  No entries in range (krmsCut0)
2302   //     b.  Big sperad          (krmscut1)
2303   //     c.  RMSto fit sigma bigger then (kmultiCut)
2304   //     d.  Too big angular spread 
2305   //  
2306
2307   const Float_t krmsCut0=0.001;
2308   const Float_t krmsCut1=0.16;
2309   const Float_t kmultiCut=2;
2310   const Float_t kcutP0=0.002;
2311   AliMagF* magF=  dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
2312   Double_t xyz[3]={90,0,10};         // tmp. global position 
2313   Double_t bxyz[3]={90,0,10};        // tmp. mag field  integral - cylindrical
2314   Double_t bgxyz[3]={90,0,10};       // tmp. mag field  integral - local
2315   //
2316   AliTPCcalibLaser *laser = this;
2317   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2318   TF1 fg("fg","gaus");
2319   AliTPCParam  * tpcparam    = 0;   
2320   // start set up for absolute residuals analysis
2321   //
2322   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
2323   tpcparam    = calib->GetParameters(); 
2324   if (!tpcparam) tpcparam    = new AliTPCParamSR;
2325   tpcparam->Update();
2326   AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2327   Float_t current=0;
2328   Float_t bfield      = 0, bz=0;
2329
2330   if (grp){
2331     Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
2332     current = grp->GetL3Current((AliGRPObject::Stats)0);
2333     bfield = polarity*5*current/30000.;
2334     bz = polarity*5*current/30000.;
2335     printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2336   }
2337
2338   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2339   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2340   TLinearFitter lfabsyInner(2);
2341   lfabsyInner.SetFormula("1 ++ x");
2342   TLinearFitter lfabszInner(2);
2343   lfabszInner.SetFormula("1 ++ x");
2344
2345   TLinearFitter lfabsyOuter(2);
2346   lfabsyOuter.SetFormula("1 ++ x");
2347   TLinearFitter lfabszOuter(2);
2348   lfabszOuter.SetFormula("1 ++ x");
2349   // end set up for absolute residuals analysis
2350
2351   //
2352   //
2353   for (Int_t id=0; id<336; id++){
2354     Bool_t isOK=kTRUE;
2355     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
2356     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2357     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
2358     TH1F * hisP3    = (TH1F*)laser->fDeltaP3.At(id);
2359     TH1F * hisP4    = (TH1F*)laser->fDeltaP4.At(id);
2360     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
2361     //if (!hisphi) continue;
2362     Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2363     //if (entries<minEntries) continue;
2364     //
2365     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2366     if (!ltrp) {
2367      AliTPCLaserTrack::LoadTracks();
2368       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2369     }
2370     ltrp->UpdatePoints();
2371     pcstream->GetFile()->cd();
2372     if (hisphi)  hisphi->Write();
2373     if (hisphiP) hisphiP->Write();
2374     if (hisZ)    hisZ->Write();
2375     if (hisP3)    hisP3->Write();
2376     if (hisP4)    hisP4->Write();
2377     
2378     Float_t meanphi = hisphi->GetMean();
2379     Float_t rmsphi = hisphi->GetRMS();
2380     //
2381     Float_t meanphiP = hisphiP->GetMean();
2382     Float_t rmsphiP = hisphiP->GetRMS();
2383     Float_t meanZ = hisZ->GetMean();
2384     Float_t rmsZ = hisZ->GetRMS();
2385     if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2386       hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2387     Double_t gphi1 = fg.GetParameter(1);
2388     Double_t gphi2 = fg.GetParameter(2);
2389     if (hisphiP->GetRMS()>0)
2390       hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2391     Double_t gphiP1 = fg.GetParameter(1);
2392     Double_t gphiP2 = fg.GetParameter(2);
2393     //
2394     if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2395       hisZ->Fit(&fg,"","");
2396     Double_t gz1 = fg.GetParameter(1);
2397     Double_t gz2 = fg.GetParameter(2);
2398     //
2399     if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2400       hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2401     Double_t gp31 = fg.GetParameter(1);
2402     Double_t gp32 = fg.GetParameter(2);
2403     Double_t meanp3 = hisP3->GetMean();
2404     Double_t rmsp3  = hisP3->GetRMS();
2405     //
2406     if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2407       hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2408     Double_t gp41 = fg.GetParameter(1);
2409     Double_t gp42 = fg.GetParameter(2);
2410     Double_t meanp4 = hisP4->GetMean();
2411     Double_t rmsp4  = hisP4->GetRMS();
2412     //
2413     Float_t meanS=hisS->GetMean();
2414     //
2415     Double_t lxyz[3];
2416     Double_t lpxyz[3];
2417     ltrp->GetXYZ(lxyz);
2418     ltrp->GetPxPyPz(lpxyz);
2419
2420     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2421     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2422     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
2423     if (gphiP2>kcutP0) isOK=kFALSE;
2424     //
2425     //
2426     //
2427     TH1 *his =0;
2428     //
2429     his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2430     Float_t mnclIn =  his->GetMean();
2431     delete his;
2432     his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2433     Float_t mnclOut =  his->GetMean();
2434     delete his;
2435     his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2436     Float_t mnclIO =  his->GetMean();
2437     delete his;
2438     his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2439     Float_t mLclIn =  his->GetMean();
2440     delete his;
2441     his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2442     Float_t mLclOut =  his->GetMean();
2443     delete his;
2444     his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2445     Float_t mLclIO =  his->GetMean();
2446     delete his;
2447     //
2448     his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2449     Float_t mdEdx =  his->GetMean();
2450     delete his;
2451     //
2452     //
2453     //
2454     //
2455     his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - line
2456     Float_t  mChi2YIn1=  his->GetMean();
2457     delete his;
2458     his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - line
2459     Float_t mChi2YOut1 =  his->GetMean();
2460     delete his;
2461     his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - parabola
2462     Float_t mChi2YIn2 =  his->GetMean();
2463     delete his;
2464     his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - parabola
2465     Float_t mChi2YOut2 =  his->GetMean();
2466     delete his;
2467     his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 y IO    - common
2468     Float_t mChi2YIO1 =  his->GetMean();
2469     delete his;
2470     his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - line
2471     Float_t mChi2ZIn1 =  his->GetMean();
2472     delete his;
2473     his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - line
2474     Float_t mChi2ZOut1 =  his->GetMean();
2475     delete his;
2476     his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - parabola
2477     Float_t mChi2ZIn2 =  his->GetMean();
2478     delete his;
2479     his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - parabola
2480     Float_t mChi2ZOut2 =  his->GetMean();
2481     delete his;
2482     his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 z IO    - common
2483     Float_t mChi2ZIO1 =  his->GetMean();
2484     delete his;
2485     //
2486     // fit res. histos
2487     // 
2488     his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2489     Float_t  edZfit  = his->GetEntries();
2490     Float_t  mdZfit =  his->GetMean();
2491     Float_t  rdZfit  = his->GetRMS();
2492     delete his;
2493
2494     his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - line
2495     Float_t  ePy1vP0  = his->GetEntries();
2496     Float_t  mPy1vP0  = his->GetMean();
2497     Float_t  rPy1vP0  = his->GetRMS();
2498     delete his;
2499     
2500     his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - parabola
2501     Float_t  ePy2vP0  = his->GetEntries();
2502     Float_t  mPy2vP0  = his->GetMean();
2503     Float_t  rPy2vP0  = his->GetRMS();
2504     delete his;
2505     
2506     his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - common parabola
2507     Float_t  ePy3vP0  = his->GetEntries();
2508     Float_t  mPy3vP0  = his->GetMean();
2509     Float_t  rPy3vP0  = his->GetRMS();
2510     delete his;
2511     
2512     his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - line
2513     Float_t  ePy1vP1  = his->GetEntries();
2514     Float_t  mPy1vP1  = his->GetMean();
2515     Float_t  rPy1vP1  = his->GetRMS();
2516     delete his;
2517     
2518     his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - parabola
2519     Float_t  ePy2vP1  = his->GetEntries();
2520     Float_t  mPy2vP1  = his->GetMean();
2521     Float_t  rPy2vP1  = his->GetRMS();
2522     delete his;
2523     
2524     his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - common parabola
2525     Float_t  ePy3vP1  = his->GetEntries();
2526     Float_t  mPy3vP1  = his->GetMean();
2527     Float_t  rPy3vP1  = his->GetRMS();
2528     delete his;
2529     
2530     his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2531     Float_t  ePy2vP2In  = his->GetEntries();
2532     Float_t  mPy2vP2In  = his->GetMean();
2533     Float_t  rPy2vP2In  = his->GetRMS();
2534     delete his;
2535     
2536     his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2537     Float_t  ePy2vP2Out  = his->GetEntries();
2538     Float_t  mPy2vP2Out  = his->GetMean();
2539     Float_t  rPy2vP2Out  = his->GetRMS();
2540     delete his;
2541     
2542     his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2543     Float_t  ePy3vP2IO  = his->GetEntries();
2544     Float_t  mPy3vP2IO  = his->GetMean();
2545     Float_t  rPy3vP2IO  = his->GetRMS();
2546     delete his;
2547     
2548     //
2549     //
2550     his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - line
2551     Float_t  ePz1vP0  = his->GetEntries();
2552     Float_t  mPz1vP0  = his->GetMean();
2553     Float_t  rPz1vP0  = his->GetRMS();
2554     delete his;
2555     
2556     his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - parabola
2557     Float_t  ePz2vP0  = his->GetEntries();
2558     Float_t  mPz2vP0  = his->GetMean();
2559     Float_t  rPz2vP0  = his->GetRMS();
2560     delete his;
2561     
2562     his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - common parabola
2563     Float_t  ePz3vP0  = his->GetEntries();
2564     Float_t  mPz3vP0  = his->GetMean();
2565     Float_t  rPz3vP0  = his->GetRMS();
2566     delete his;
2567     
2568     his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - line
2569     Float_t  ePz1vP1  = his->GetEntries();
2570     Float_t  mPz1vP1  = his->GetMean();
2571     Float_t  rPz1vP1  = his->GetRMS();
2572     delete his;
2573     
2574     his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - parabola
2575     Float_t  ePz2vP1  = his->GetEntries();
2576     Float_t  mPz2vP1  = his->GetMean();
2577     Float_t  rPz2vP1  = his->GetRMS();
2578     delete his;
2579     
2580     his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - common parabola
2581     Float_t  ePz3vP1  = his->GetEntries();
2582     Float_t  mPz3vP1  = his->GetMean();
2583     Float_t  rPz3vP1  = his->GetRMS();
2584     delete his;
2585     
2586     his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2587     Float_t  ePz2vP2In  = his->GetEntries();
2588     Float_t  mPz2vP2In  = his->GetMean();
2589     Float_t  rPz2vP2In  = his->GetRMS();
2590     delete his;
2591     
2592     his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2593     Float_t  ePz2vP2Out  = his->GetEntries();
2594     Float_t  mPz2vP2Out  = his->GetMean();
2595     Float_t  rPz2vP2Out  = his->GetRMS();
2596     delete his;
2597     
2598     his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2599     Float_t  ePz3vP2IO  = his->GetEntries();
2600     Float_t  mPz3vP2IO  = his->GetMean();
2601     Float_t  rPz3vP2IO  = his->GetRMS();
2602     delete his;
2603
2604     // Fit absolute laser residuals
2605     TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2606     TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2607
2608     Int_t secInner = TMath::Nint(fBeamSectorInner[id]); 
2609     Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]); 
2610
2611     TVectorD vecX(159);        // X
2612     TVectorD vecY(159);        // Y
2613     TVectorD vecR(159);        // R
2614     TVectorD vecDY(159);       // absolute residuals in Y
2615     TVectorD vecDZ(159);       // absolute residuals in Z
2616     TVectorD vecN(159);        // number of clusters
2617     TVectorD vecEy(159);       //error y
2618     TVectorD vecEz(159);       //error z
2619     TVectorD vecPhi(159);      // local tangent
2620     TVectorD vecPhiR(159);     // local tangent
2621     // magnetic field integrals
2622     TVectorD vecIBR(159);        // radial
2623     TVectorD vecIBRPhi(159);     // r-phi
2624     TVectorD vecIBLX(159);       // local x
2625     TVectorD vecIBLY(159);       // local y
2626     TVectorD vecIBGX(159);       // local x
2627     TVectorD vecIBGY(159);       // local y
2628     TVectorD vecIBZ(159);        // z
2629     //
2630     for (Int_t irow=0;irow<159;irow++){
2631       vecIBR[irow]=0;
2632       vecIBRPhi[irow]=0;
2633       vecIBLX[irow]=0;
2634       vecIBLY[irow]=0;
2635       vecIBGX[irow]=0;
2636       vecIBGY[irow]=0;
2637       vecIBZ[irow]=0;
2638       Double_t gx    =(*(ltrp->fVecGX))[irow];
2639       Double_t gy    =(*(ltrp->fVecGY))[irow];
2640       Int_t    lsec  =TMath::Nint((*(ltrp->fVecSec))[irow]);
2641       Double_t   ca  =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2642       Double_t   sa  =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2643       xyz[2]=(*(ltrp->fVecGZ))[irow];
2644       xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2645       xyz[1]=TMath::ATan2(gy,gx);
2646       Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2647       if (magF){
2648         magF->GetTPCIntCyl(xyz,bxyz);
2649         magF->GetTPCInt(gxyz,bgxyz);
2650         vecIBR[irow]=bxyz[0];
2651         vecIBRPhi[irow]=bxyz[1];
2652         //
2653         vecIBGX[irow]=bgxyz[0];
2654         vecIBGY[irow]=bgxyz[1];
2655         //
2656         vecIBLX[irow]=  bgxyz[0]*ca+bgxyz[1]*sa;
2657         vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2658         //
2659
2660         vecIBZ[irow]=bxyz[2];
2661       }
2662     }
2663
2664
2665     lfabsyInner.ClearPoints();    
2666     lfabszInner.ClearPoints();    
2667     lfabsyOuter.ClearPoints();    
2668     lfabszOuter.ClearPoints();    
2669     // dummy fit values
2670     Int_t    nClInY         = 0;
2671     Float_t  yAbsInOffset  = -100;
2672     Float_t  yAbsInSlope   = -100;
2673     Float_t  yAbsInDeltay  = -100;
2674     Int_t    nClInZ         = 0;
2675     Float_t  zAbsInOffset  = -100;
2676     Float_t  zAbsInSlope   = -100;
2677     Float_t  zAbsInDeltaz  = -100;
2678     Int_t    nClOutY         = 0;
2679     Float_t  yAbsOutOffset  = -100;
2680     Float_t  yAbsOutSlope   = -100;
2681     Float_t  yAbsOutDeltay  = -100;
2682     Int_t    nClOutZ         = 0;
2683     Float_t  zAbsOutOffset  = -100;
2684     Float_t  zAbsOutSlope   = -100;
2685     Float_t  zAbsOutDeltaz  = -100;
2686
2687     Float_t  lasTanPhiLocIn = -100;
2688     Float_t  lasTanPhiLocOut = -100;
2689
2690     if(histAbsY && histAbsY->GetEntries()>0) {
2691       
2692       Double_t rotAngOut = 10;
2693       Double_t rotAngIn = 10;
2694       if((secInner%36)!=(secOuter%36))
2695         rotAngIn += 20; // 30 degrees
2696       
2697       // Calculate laser mirror X position in local frame
2698       Double_t laserposOut = 
2699         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2700       Double_t laserposIn = 
2701         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2702       
2703       // Calculate laser tan phi in local frame
2704       lasTanPhiLocOut  = TMath::ASin(ltrp->GetSnp());
2705       if(lasTanPhiLocOut<0) {
2706         lasTanPhiLocIn  = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2707         lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2708       } else {
2709         
2710         lasTanPhiLocIn  = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2711         lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2712       }
2713       
2714       lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2715       lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2716       
2717       TProfile* yprof  = histAbsY->ProfileX("yprof");
2718       TProfile* zprof  = histAbsZ->ProfileX("zprof");
2719       
2720       for(Int_t bin = 1; bin<=159; bin++) {
2721         
2722         if(yprof->GetBinEntries(bin)<5&&
2723            zprof->GetBinEntries(bin)<5) {
2724           continue;     
2725         }
2726         
2727         // There is a problem in defining inner and outer sectors for
2728         // the outer beams (0 and 6) where both sectors are OROCs. To
2729         // make sure there is no overlap row 94 to 99 are cutted.
2730         if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2731           continue;
2732         
2733         Int_t row = (bin-1);
2734         if(row>62)
2735           row -= 63;
2736
2737         Bool_t isOuter = kTRUE;
2738         Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2739         
2740         if(bin<=62 ||                                        
2741            (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2742
2743           isOuter = kFALSE;
2744           sector = TMath::Nint(fBeamSectorInner[id]);
2745         }
2746
2747
2748         Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2749         vecN[bin-1] =yprof->GetBinEntries(bin);
2750         vecEy[bin-1]=yprof->GetBinError(bin);
2751         vecEz[bin-1]=zprof->GetBinError(bin);
2752         vecX[bin-1] = x;
2753         vecDY[bin-1] = yprof->GetBinContent(bin);
2754         vecDZ[bin-1] = zprof->GetBinContent(bin);
2755         if (bin>0&&bin<159){
2756           //
2757           //truncated mean - skip first and the last  pad row
2758           //
2759           Int_t firstBin=TMath::Max(bin-5,0);
2760           Int_t lastBin =TMath::Min(bin+5,158);
2761           histAbsY->GetXaxis()->SetRangeUser(firstBin,lastBin);
2762           histAbsY->GetYaxis()->SetRangeUser(-2,2);
2763           vecEy[bin-1]=histAbsY->GetRMS(2);
2764           vecDY[bin-1]=histAbsY->GetMean(2);
2765           histAbsY->GetXaxis()->SetRangeUser(firstBin+2,lastBin-2);//use+-2 bins
2766           histAbsY->GetYaxis()->SetRangeUser(vecDY[bin-1]-4*vecEy[bin-1],vecDY[bin-1]+4*vecEy[bin-1]);
2767           if (yprof->GetBinEntries(bin-1)>0) vecEy[bin-1]=histAbsY->GetRMS(2)/TMath::Sqrt(yprof->GetBinEntries(bin-1));
2768           vecDY[bin-1]=histAbsY->GetMean(2);
2769         }
2770
2771         if(!isOuter) { // inner   
2772           vecPhi[bin-1]=lasTanPhiLocIn;
2773           // Calculate local y from residual and database
2774           Double_t y  = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2775             + vecDY[bin-1];
2776           vecY[bin-1] = y;
2777           Double_t r  = TMath::Sqrt(x*x + y*y);
2778           vecR[bin-1] = r;
2779           // Find angle between laser vector and R vector
2780           // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2781           Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2782           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2783           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2784           if(lasTanPhiLocIn<0)
2785             vecPhiR[bin-1]*=-1; // to have the same sign
2786           
2787           if(yprof->GetBinEntries(bin)>=10) {
2788             lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin), 
2789                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2790           } 
2791           if(zprof->GetBinEntries(bin)>=10) {
2792             lfabszInner.AddPoint(&x, zprof->GetBinContent(bin), 
2793                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2794           }
2795         } else { // outer
2796           vecPhi[bin-1]=lasTanPhiLocOut;          
2797           // Calculate local y from residual and database
2798           Double_t y  = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2799             + vecDY[bin-1];
2800           vecY[bin-1] = y;
2801           Double_t r  = TMath::Sqrt(x*x + y*y);
2802           vecR[bin-1] = r;
2803
2804           Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2805           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2806           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2807           if(lasTanPhiLocOut<0)
2808             vecPhiR[bin-1]*=-1; // to have the same sign
2809
2810           if(yprof->GetBinEntries(bin)>=10) {
2811             lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin), 
2812                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2813           }
2814           if(zprof->GetBinEntries(bin)>=10) {
2815             lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin), 
2816                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2817           }
2818         }
2819         // global position
2820         
2821       }
2822         
2823       delete yprof; delete zprof;
2824
2825       
2826       // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2827       nClInY = lfabsyInner.GetNpoints();
2828       if(lfabsyInner.GetNpoints()>10) {
2829         lfabsyInner.EvalRobust(0.95);
2830         
2831         TVectorD result(2);
2832         lfabsyInner.GetParameters(result);
2833         yAbsInOffset = result[0];
2834         yAbsInSlope  = result[1];
2835         yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2836       }
2837       nClInZ = lfabszInner.GetNpoints();
2838       if(lfabszInner.GetNpoints()>10) {
2839         lfabszInner.EvalRobust(0.95);
2840         
2841         TVectorD result(2);
2842         lfabszInner.GetParameters(result);
2843         zAbsInOffset = result[0];
2844         zAbsInSlope  = result[1];
2845         zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2846       } 
2847       nClOutY = lfabsyOuter.GetNpoints();
2848       if(lfabsyOuter.GetNpoints()>10) {
2849         lfabsyOuter.EvalRobust(0.95);
2850         
2851         TVectorD result(2);
2852         lfabsyOuter.GetParameters(result);
2853         yAbsOutOffset = result[0];
2854         yAbsOutSlope  = result[1];
2855         yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2856       }
2857       nClOutZ = lfabszOuter.GetNpoints();
2858       if(lfabszOuter.GetNpoints()>10) {
2859         lfabszOuter.EvalRobust(0.95);
2860         
2861         TVectorD result(2);
2862         lfabszOuter.GetParameters(result);
2863         zAbsOutOffset = result[0];
2864         zAbsOutSlope  = result[1];
2865         zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2866       }
2867     }
2868
2869     
2870     Int_t itime=-1;
2871     Float_t  coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2872     Float_t  coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2873     Float_t  coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2874     Float_t  coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2875     //
2876     Float_t  skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2877     Float_t  skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2878     //
2879     Float_t  ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2880     Float_t  ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2881
2882     //
2883     (*pcstream)<<"Mean"<<
2884       "run="<<run<<               //
2885       //voltages
2886       "UcIA="<<coverIA<<
2887       "UcIC="<<coverIC<<
2888       "UcOA="<<coverOA<<
2889       "UcOC="<<coverOC<<
2890       "UsA="<<skirtA<<
2891       "UsC="<<skirtC<<
2892       "UggA="<<ggOffA<<
2893       "UggC="<<ggOffC<<
2894       //
2895       "isOK="<<isOK<<             //
2896       "id="<<id<<                 // track id
2897       "entries="<<entries<<       // number of entries
2898       "bz="<<bfield<<             // bfield
2899       "LTr.="<<ltrp<<             // refernece track
2900       //
2901       "mnclIn="<<mnclIn<<         // mean number of clusters in inner
2902       "mnclOut="<<mnclOut<<       // mean number of clusters in outer
2903       "mnclIO="<<mnclIO<<         // mean number of clusters in inner+outer
2904       "mLclIn="<<mLclIn<<         // mean number of clusters in inner
2905       "mLclOut="<<mLclOut<<       // mean number of clusters in outer
2906       "mLclIO="<<mLclIO<<         // mean number of clusters in inner+outer
2907       "mdEdx="<<mdEdx<<           // mean dEdx
2908       "edZfit="<<edZfit<<           // entries z fit
2909       "mdZfit="<<mdZfit<<           // mean z fit
2910       "rdZfit="<<rdZfit<<           // RMS z fit
2911       //
2912       //
2913       "mChi2YIn1="<<mChi2YIn1<<       //->chi2 y inner - line
2914       "mChi2YOut1="<<mChi2YOut1<<     //->chi2 y inner - line
2915       "mChi2YIn2="<<mChi2YIn2<<       //->chi2 y inner - parabola
2916       "mChi2YOut2="<<mChi2YOut2<<     //->chi2 y inner - parabola
2917       "mChi2YIO1="<<mChi2YIO1<<       //->chi2 y IO    - common
2918       "mChi2ZIn1="<<mChi2ZIn1<<       //->chi2 z inner - line
2919       "mChi2ZOut1="<<mChi2ZOut1<<     //->chi2 z inner - line
2920       "mChi2ZIn2="<<mChi2ZIn2<<       //->chi2 z inner - parabola
2921       "mChi2ZOut2="<<mChi2ZOut2<<     //->chi2 z inner - parabola
2922       "mChi2ZIO1="<<mChi2ZIO1<<       //->chi2 z IO    - common
2923       //
2924       //
2925       "ePy1vP0="<<ePy1vP0<<
2926       "mPy1vP0="<<mPy1vP0<<
2927       "rPy1vP0="<<rPy1vP0<<
2928       "ePy2vP0="<<ePy2vP0<<
2929       "mPy2vP0="<<mPy2vP0<<
2930       "rPy2vP0="<<rPy2vP0<<      
2931       "ePy3vP0="<<ePy3vP0<<
2932       "mPy3vP0="<<mPy3vP0<<
2933       "rPy3vP0="<<rPy3vP0<<
2934       "ePy1vP1="<<ePy1vP1<<
2935       "mPy1vP1="<<mPy1vP1<<
2936       "rPy1vP1="<<rPy1vP1<<
2937       "ePy2vP1="<<ePy2vP1<<
2938       "mPy2vP1="<<mPy2vP1<<
2939       "rPy2vP1="<<rPy2vP1<<
2940       "ePy3vP1="<<ePy3vP1<<
2941       "mPy3vP1="<<mPy3vP1<<
2942       "rPy3vP1="<<rPy3vP1<<
2943       "ePy2vP2In="<<ePy2vP2In<<
2944       "mPy2vP2In="<<mPy2vP2In<<
2945       "rPy2vP2In="<<rPy2vP2In<<
2946       "ePy2vP2Out="<<ePy2vP2Out<<
2947       "mPy2vP2Out="<<mPy2vP2Out<<
2948       "rPy2vP2Out="<<rPy2vP2Out<<
2949       "ePy3vP2IO="<<ePy3vP2IO<<
2950       "mPy3vP2IO="<<mPy3vP2IO<<
2951       "rPy3vP2IO="<<rPy3vP2IO<<
2952       //
2953       //
2954       "ePz1vP0="<<ePz1vP0<<
2955       "mPz1vP0="<<mPz1vP0<<
2956       "rPz1vP0="<<rPz1vP0<<
2957       "ePz2vP0="<<ePz2vP0<<
2958       "mPz2vP0="<<mPz2vP0<<
2959       "rPz2vP0="<<rPz2vP0<<
2960       "ePz3vP0="<<ePz3vP0<<
2961       "mPz3vP0="<<mPz3vP0<<
2962       "rPz3vP0="<<rPz3vP0<<
2963       "ePz1vP1="<<ePz1vP1<<
2964       "mPz1vP1="<<mPz1vP1<<
2965       "rPz1vP1="<<rPz1vP1<<
2966       "ePz2vP1="<<ePz2vP1<<
2967       "mPz2vP1="<<mPz2vP1<<
2968       "rPz2vP1="<<rPz2vP1<<
2969       "ePz3vP1="<<ePz3vP1<<
2970       "mPz3vP1="<<mPz3vP1<<
2971       "rPz3vP1="<<rPz3vP1<<
2972       "ePz2vP2In="<<ePz2vP2In<<
2973       "mPz2vP2In="<<mPz2vP2In<<
2974       "rPz2vP2In="<<rPz2vP2In<<
2975       "ePz2vP2Out="<<ePz2vP2Out<<
2976       "mPz2vP2Out="<<mPz2vP2Out<<
2977       "rPz2vP2Out="<<rPz2vP2Out<<
2978       "ePz3vP2IO="<<ePz3vP2IO<<  
2979       "mPz3vP2IO="<<mPz3vP2IO<<
2980       "rPz3vP2IO="<<rPz3vP2IO<<     
2981       //
2982       //
2983       //
2984       "lx0="<<lxyz[0]<<            // reference x
2985       "lx1="<<lxyz[1]<<            // reference y
2986       "lx2="<<lxyz[2]<<            // refernece z
2987       "lpx0="<<lpxyz[0]<<           // reference x
2988       "lpx1="<<lpxyz[1]<<          // reference y
2989       "lpx2="<<lpxyz[2]<<          // refernece z
2990       //
2991       "msig="<<meanS<<
2992       //
2993       "mphi="<<meanphi<<         //
2994       "rmsphi="<<rmsphi<<        //
2995       "gphi1="<<gphi1<<
2996       "gphi2="<<gphi2<<
2997       //
2998       "mphiP="<<meanphiP<<       //
2999       "rmsphiP="<<rmsphiP<<      //
3000       "gphiP1="<<gphiP1<<
3001       "gphiP2="<<gphiP2<<
3002       //
3003       "meanZ="<<meanZ<<
3004       "rmsZ="<<rmsZ<<
3005       "gz1="<<gz1<<
3006       "gz2="<<gz2<<
3007       //
3008       "gp31="<<gp31<<            //gaus mean - tgl
3009       "gp32="<<gp32<<            //gaus rms  -tgl
3010       "meanp3="<<meanp3<<            //mean - tgl
3011       "rmsp3="<<rmsp3<<            //rms  -tgl
3012       "gp41="<<gp41<<            //gaus mean - P4
3013       "gp42="<<gp42<<            //gaus rms  - P4
3014       "meanp4="<<meanp4<<            //mean - P4
3015       "rmsp4="<<rmsp4<<            //rms  - P4
3016       // Parameters from abs res analysis
3017       "SecIn="<<secInner<<              // inner sector
3018       "SecOut="<<secOuter<<             // outer sector
3019       "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
3020       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
3021       "ibr.="<<&vecIBR<<   // radial filed integral
3022       "ibrphi.="<<&vecIBRPhi<<   // r=phifiled integral
3023       "ibr.="<<&vecIBR<<   // radial filed integral
3024       "ibz.="<<&vecIBZ<<   // z filed integral
3025       //
3026       "iblx.="<<&vecIBLX<<   // local bx  integral
3027       "ibly.="<<&vecIBLY<<   // local by integral
3028       "ibgx.="<<&vecIBGX<<   // global bx  integral
3029       "ibgy.="<<&vecIBGY<<   // global by integral
3030       //
3031       "X.="<<&vecX<<       // local x 
3032       "Y.="<<&vecY<<       // local y 
3033       "R.="<<&vecR<<       // radius 
3034       "dY.="<<&vecDY<<     // abs y residuals
3035       "dZ.="<<&vecDZ<<     // abs z residuals
3036       "eY.="<<&vecEy<<     // error of y residuals
3037       "eZ.="<<&vecEz<<     // error of z residuals
3038       "kY.="<<&vecPhi<<    // local tangent y (fixed for sector)
3039       "kYR.="<<&vecPhiR<<  // tangent between laser and R vector (varies inside sector) 
3040       "nCl.="<<&vecN<<     // number of clusters
3041       //
3042       "nClInY="<<nClInY<<               // Number of clusters for inner y
3043       "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
3044       "yAbsInSlope="<<yAbsInSlope <<  // fitted slope absres (inner y)
3045       "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
3046       "nClInZ="<<nClInZ<<               // Number of clusters for inner z
3047       "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
3048       "zAbsInSlope="<<zAbsInSlope <<  // fitted slope absres (inner z)
3049       "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
3050       //
3051       "nClOutY="<<nClOutY<<               // Number of clusters for outer y
3052       "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
3053       "yAbsOutSlope="<<yAbsOutSlope <<  // fitted slope absres (outer y)
3054       "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
3055       "nClOutZ="<<nClOutZ<<               // Number of clusters for outer z
3056       "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
3057       "zAbsOutSlope="<<zAbsOutSlope <<  // fitted slope absres (outer z)
3058       "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
3059       //
3060       "\n";
3061   }
3062   delete pcstream;
3063   /*
3064     Browse the content
3065     TFile fmean("laserMean.root");
3066     
3067
3068    */
3069
3070
3071 }
3072
3073
3074
3075 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3076   //
3077   //
3078   //
3079   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3080   TFile * f = pcstream->GetFile();
3081   f->mkdir("dirphi");
3082   f->mkdir("dirphiP");
3083   f->mkdir("dirZ");
3084   TF1 fp("p1","pol1");
3085   //
3086   //
3087   char    cut[1000];
3088   char grname[1000];
3089   char grnamefull[1000];
3090
3091   Double_t mphi[100];
3092   Double_t mphiP[100];
3093   Double_t smphi[100];
3094   Double_t smphiP[100];
3095   Double_t mZ[100];
3096   Double_t smZ[100];
3097   Double_t bz[100];
3098   Double_t sbz[100];
3099   // fit parameters
3100   Double_t pphi[3];
3101   Double_t pphiP[3];
3102   Double_t pmZ[3];
3103   
3104   //
3105   for (Int_t id=0; id<336; id++){
3106     // id =205;
3107     snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
3108     Int_t entries = chain->Draw("bz",cut,"goff");
3109     if (entries<3) continue;
3110     AliTPCLaserTrack *ltrp = 0;
3111     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3112     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3113     Double_t lxyz[3];
3114     Double_t lpxyz[3];
3115     ltrp->GetXYZ(lxyz);
3116     ltrp->GetPxPyPz(lpxyz);
3117
3118     chain->Draw("bz",cut,"goff");
3119     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3120     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3121     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3122     //
3123     chain->Draw("gphi1",cut,"goff");
3124     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3125     chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3126     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3127     //
3128     chain->Draw("gphiP1",cut,"goff");
3129     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3130     chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3131     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3132     //
3133     chain->Draw("gz1",cut,"goff");
3134     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3135     chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3136     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3137     //
3138     //
3139     snprintf(grnamefull,1000,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3140             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3141     // store data  
3142     // phi
3143     f->cd("dirphi");
3144     Float_t phiB0 =0;
3145     Float_t phiPB0=0;
3146     Float_t zB0=0;
3147     for (Int_t ientry=0; ientry<entries; ientry++){
3148       if (TMath::Abs(bz[ientry])<0.05){
3149         phiB0  = mphi[ientry];
3150         phiPB0 = mphiP[ientry];
3151         zB0    = mZ[ientry];
3152       }
3153     }
3154     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3155     grphi->Draw("a*");
3156     grphi->Fit(&fp);
3157     pphi[0] = fp.GetParameter(0);                          // offset
3158     pphi[1] = fp.GetParameter(1);                          // slope
3159     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3160     snprintf(grname,1000,"phi_id%d",id);
3161     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
3162     grphi->GetXaxis()->SetTitle("b_{z} (T)");
3163     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3164     grphi->SetMaximum(1.2);
3165     grphi->SetMinimum(-1.2);
3166     grphi->Draw("a*");
3167
3168     grphi->Write();
3169     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3170     // phiP
3171     f->cd("dirphiP)");
3172     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3173     grphiP->Draw("a*");
3174     grphiP->Fit(&fp);
3175     pphiP[0] = fp.GetParameter(0);                          // offset
3176     pphiP[1] = fp.GetParameter(1);                          // slope
3177     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3178     snprintf(grname,1000,"phiP_id%d",id);
3179     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
3180     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3181     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3182     grphiP->SetMaximum(pphiP[0]+0.005);
3183     grphiP->SetMinimum(pphiP[0]-0.005);
3184
3185     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3186     grphiP->Write();
3187     //
3188     //Z
3189     f->cd("dirZ");
3190     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3191     grmZ->Draw("a*");
3192     grmZ->Fit(&fp);
3193     pmZ[0] = fp.GetParameter(0);                          // offset
3194     pmZ[1] = fp.GetParameter(1);                          // slope
3195     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3196     snprintf(grname,1000,"mZ_id%d",id);
3197     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
3198     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3199     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3200
3201     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3202     grmZ->Write();
3203     //
3204     // P4
3205     //
3206
3207     for (Int_t ientry=0; ientry<entries; ientry++){
3208       (*pcstream)<<"Mean"<<
3209         "id="<<id<<
3210         "LTr.="<<ltrp<<
3211         "entries="<<entries<<
3212         "bz="<<bz[ientry]<<
3213         "lx0="<<lxyz[0]<<          // reference x
3214         "lx1="<<lxyz[1]<<          // reference y
3215         "lx2="<<lxyz[2]<<          // refernece z      
3216         "lpx0="<<lpxyz[0]<<          // reference x
3217         "lpx1="<<lpxyz[1]<<          // reference y
3218         "lpx2="<<lpxyz[2]<<          // refernece z            
3219         //values
3220         "phiB0="<<phiB0<<          // position shift at 0 field
3221         "phiPB0="<<phiPB0<<        // angular  shift at 0 field
3222         "zB0="<<zB0<<              // z shift for 0 field
3223         //
3224         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3225         "pphi0="<<pphi[0]<<   // offset
3226         "pphi1="<<pphi[1]<<   // slope
3227         "pphi2="<<pphi[2]<<   // norm chi2
3228         //
3229         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3230         "pphiP0="<<pphiP[0]<< // offset
3231         "pphiP1="<<pphiP[1]<< // slope
3232         "pphiP2="<<pphiP[2]<< // norm chi2
3233         //
3234         "gz1="<<mZ[ientry]<<
3235         "pmZ0="<<pmZ[0]<<     // offset
3236         "pmZ1="<<pmZ[1]<<     // slope
3237         "pmZ2="<<pmZ[2]<<     // norm chi2
3238         "\n";
3239     }
3240   }
3241   
3242   delete pcstream;
3243   
3244 }
3245
3246
3247 void AliTPCcalibLaser::Analyze(){
3248   //
3249   //
3250   //
3251 }
3252
3253
3254 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3255
3256   TIterator* iter = li->MakeIterator();
3257   AliTPCcalibLaser* cal = 0;
3258   static Int_t counter0=0;
3259   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3260     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3261       //      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3262       return -1;
3263     }
3264     printf("Marging number %d\n", counter0);
3265     counter0++;
3266     //
3267     MergeFitHistos(cal);
3268     TH1F *h=0x0;
3269     TH1F *hm=0x0;
3270     TH2F *h2=0x0;
3271     TH2F *h2m=0x0;
3272     //    TProfile *hp=0x0;
3273     //TProfile *hpm=0x0;
3274
3275     for (Int_t id=0; id<336; id++){
3276       // merge fDeltaZ histograms
3277       hm = (TH1F*)cal->fDeltaZ.At(id);
3278       h  = (TH1F*)fDeltaZ.At(id);
3279       if (!h &&hm &&hm->GetEntries()>0) {
3280         h=(TH1F*)hm->Clone();
3281         h->SetDirectory(0);
3282         fDeltaZ.AddAt(h,id);      
3283       }
3284       if (h && hm) h->Add(hm);
3285
3286       // merge fP3 histograms
3287       hm = (TH1F*)cal->fDeltaP3.At(id);
3288       h  = (TH1F*)fDeltaP3.At(id);
3289       if (!h&&hm &&hm->GetEntries()>0) {
3290         h=(TH1F*)hm->Clone();
3291         h->SetDirectory(0);
3292         fDeltaP3.AddAt(h,id);
3293       }
3294       if (h && hm) h->Add(hm);
3295       // merge fP4 histograms
3296       hm = (TH1F*)cal->fDeltaP4.At(id);
3297       h  = (TH1F*)fDeltaP4.At(id);
3298       if (!h &&hm &&hm->GetEntries()>0) {
3299         h=(TH1F*)hm->Clone();
3300         h->SetDirectory(0);
3301         fDeltaP4.AddAt(h,id);
3302       }
3303       if (h&&hm) h->Add(hm);
3304
3305       //
3306       // merge fDeltaPhi histograms
3307       hm = (TH1F*)cal->fDeltaPhi.At(id);
3308       h  = (TH1F*)fDeltaPhi.At(id);
3309       if (!h &&hm &&hm->GetEntries()>0) {
3310         h= (TH1F*)hm->Clone();
3311         h->SetDirectory(0);
3312         fDeltaPhi.AddAt(h,id);
3313       }
3314       if (h&&hm) h->Add(hm);
3315       // merge fDeltaPhiP histograms
3316       hm = (TH1F*)cal->fDeltaPhiP.At(id);
3317       h  = (TH1F*)fDeltaPhiP.At(id);
3318       if (!h&&hm &&hm->GetEntries()>0) {
3319         h=(TH1F*)hm->Clone();
3320         h->SetDirectory(0);
3321         fDeltaPhiP.AddAt(h,id);
3322       }
3323       if (h&&hm) h->Add(hm);
3324       // merge fSignals histograms
3325       hm = (TH1F*)cal->fSignals.At(id);
3326       h  = (TH1F*)fSignals.At(id);
3327       if (!h&&hm &&hm->GetEntries()>0) {
3328         h=(TH1F*)hm->Clone();
3329         h->SetDirectory(0);
3330         fSignals.AddAt(h,id);
3331       }
3332       if (h&&hm) h->Add(hm);
3333       //
3334       //
3335       // merge ProfileY histograms -0
3336       h2m = (TH2F*)cal->fDeltaYres.At(id);
3337       h2  = (TH2F*)fDeltaYres.At(id);
3338       if (h2m&&h2) h2->Add(h2m);
3339       //
3340       h2m = (TH2F*)cal->fDeltaZres.At(id);
3341       h2  = (TH2F*)fDeltaZres.At(id);
3342       if (h2m&&h2) h2->Add(h2m);
3343       // merge ProfileY histograms - 2
3344       h2m = (TH2F*)cal->fDeltaYres2.At(id);
3345       h2  = (TH2F*)fDeltaYres2.At(id);
3346       if (h2m&&h2) h2->Add(h2m);
3347       //
3348       h2m = (TH2F*)cal->fDeltaZres2.At(id);
3349       h2  = (TH2F*)fDeltaZres2.At(id);
3350       if (h2m&&h2) h2->Add(h2m);
3351
3352       // merge ProfileY histograms - abs
3353       h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3354       h2  = (TH2F*)fDeltaYresAbs.At(id);
3355       if (h2m&&h2) h2->Add(h2m);
3356       if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaYresAbs.AddAt(h2,id);}