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