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