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