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