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