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