Minor bug fixes
[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   //
30   //
31   // To make laser scan the user interaction neccessary
32   //
33   .x ~/UliStyle.C
34   gSystem->Load("libANALYSIS");
35   gSystem->Load("libTPCcalib");
36   TFile fcalib("CalibObjects.root");
37   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
38   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
39   laser->DumpMeanInfo(-0.4)
40   TFile fmean("laserMean.root")
41   //
42   //  laser track clasification;
43   //
44   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
45   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
46   TCut cutN("cutN","fTPCncls>70");
47   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
48   TCut cutA = cutT+cutPt+cutP;
49   TFile f("laserTPCDebug.root");
50   TTree * treeT = (TTree*)f.Get("Track");
51   //
52   //
53   // Analyze  LASER scan 
54   //
55   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
56   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
57   AliXRDPROOFtoolkit tool;
58   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
59   chain->Lookup();
60   AliTPCcalibLaser::DumpScanInfo(chain)
61   TFile fscan("laserScan.root")
62   TTree * treeT = (TTree*)fscan.Get("Mean")
63
64 */
65
66
67
68 #include "TLinearFitter.h"
69 #include "AliTPCcalibLaser.h"
70 #include "AliExternalTrackParam.h"
71 #include "AliESDEvent.h"
72 #include "AliESDfriend.h"
73 #include "AliESDtrack.h"
74 #include "AliTPCTracklet.h"
75 #include "TH1D.h"
76 #include "TVectorD.h"
77 #include "TTreeStream.h"
78 #include "TFile.h"
79 #include "TF1.h"
80 #include "TGraphErrors.h"
81 #include "AliTPCclusterMI.h"
82 #include "AliTPCseed.h"
83 #include "AliTracker.h"
84 #include "AliLog.h"
85 #include "TClonesArray.h"
86 #include "TPad.h"
87
88
89 #include "TTreeStream.h"
90 #include <iostream>
91 #include <sstream>
92 #include "AliTPCLaserTrack.h"
93
94 using namespace std;
95
96 ClassImp(AliTPCcalibLaser)
97
98 AliTPCcalibLaser::AliTPCcalibLaser():
99   AliTPCcalibBase(),
100   fESD(0),
101   fESDfriend(),
102   fTracksMirror(336),
103   fTracksEsd(336),
104   fTracksEsdParam(336),
105   fTracksTPC(336),
106   fDeltaZ(336),
107   fDeltaPhi(336),
108   fDeltaPhiP(336),
109   fSignals(336),  
110   fFitAside(new TVectorD(3)),      
111   fFitCside(new TVectorD(3)),      
112   fEdgeXcuts(5),    
113   fEdgeYcuts(5),    
114   fNClCuts(5),      
115   fNcuts(0),        
116   fRun(0),
117   fEvent(0)
118 {
119   //
120   // Constructor
121   //
122   fTracksEsdParam.SetOwner(kTRUE);
123 }
124
125 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
126   AliTPCcalibBase(),
127   fESD(0),
128   fESDfriend(0),
129   fTracksMirror(336),
130   fTracksEsd(336),
131   fTracksEsdParam(336),
132   fTracksTPC(336),
133   fDeltaZ(336),          // array of histograms of delta z for each track
134   fDeltaPhi(336),          // array of histograms of delta z for each track
135   fDeltaPhiP(336),          // array of histograms of delta z for each track
136   fSignals(336),           // array of dedx signals
137   fFitAside(new TVectorD(3)),        // drift fit - A side
138   fFitCside(new TVectorD(3)),        // drift fit - C- side
139   fEdgeXcuts(5),       // cuts in local x direction; used in the refit of the laser tracks
140   fEdgeYcuts(5),       // cuts in local y direction; used in the refit of the laser tracks
141   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
142   fNcuts(0),           // number of cuts
143   fRun(0),
144   fEvent(0)
145 {
146   SetName(name);
147   SetTitle(title);
148   //
149   // Constructor
150   //
151   fTracksEsdParam.SetOwner(kTRUE);
152 }
153
154 AliTPCcalibLaser::~AliTPCcalibLaser() {
155   //
156   // destructor
157   //
158 }
159
160
161
162 void AliTPCcalibLaser::Process(AliESDEvent * event) {
163   //
164   //
165   // Loop over tracks and call  Process function
166   //
167   fESD = event;
168   if (!fESD) {
169     return;
170   }
171   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
172   if (!fESDfriend) {
173     return;
174   }
175   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
176   fTracksTPC.Clear();
177   fTracksEsd.Clear();
178   fTracksEsdParam.Delete();
179   //
180   Int_t n=fESD->GetNumberOfTracks();
181   Int_t run = fESD->GetRunNumber();
182   fRun = run;
183   for (Int_t i=0;i<n;++i) {
184     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
185     AliESDtrack *track=fESD->GetTrack(i);
186     TObject *calibObject=0;
187     AliTPCseed *seed=0;
188     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
189       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
190         break;
191     if (track&&seed) FindMirror(track,seed);
192     //
193   }
194
195   FitDriftV();
196   MakeDistHisto();
197   //
198   for (Int_t id=0; id<336; id++){
199     //
200     //
201     if (!fTracksEsdParam.At(id)) continue;
202     DumpLaser(id);
203 //    RefitLaser(id);
204     RefitLaserJW(id);
205
206   }
207 //  fEvent++;
208 }
209
210 void AliTPCcalibLaser::MakeDistHisto(){
211   //
212   //
213   //
214   for (Int_t id=0; id<336; id++){
215     //
216     //
217     if (!fTracksEsdParam.At(id)) continue;
218     if (!AcceptLaser(id)) continue;
219     //
220     //
221     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
222     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
223     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
224     TH1F * hisSignal = (TH1F*)fSignals.At(id);
225
226     if (!hisdz){
227       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
228       hisdz->SetDirectory(0);
229       fDeltaZ.AddAt(hisdz,id);
230       //
231       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
232       hisdphi->SetDirectory(0);
233       fDeltaPhi.AddAt(hisdphi,id);
234       //
235       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
236       hisdphiP->SetDirectory(0);
237       fDeltaPhiP.AddAt(hisdphiP,id);
238       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
239       hisSignal->SetDirectory(0);
240       fSignals.AddAt(hisSignal,id);
241     }
242
243     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
244     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
245     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
246     if (!param) return;
247     if (!ltrp) return;
248     if (!track) return;
249     Double_t xyz[3];
250     Double_t pxyz[3];
251     Double_t lxyz[3];
252     Double_t lpxyz[3];
253     param->GetXYZ(xyz);
254     param->GetPxPyPz(pxyz);
255     ltrp->GetXYZ(lxyz);
256     ltrp->GetPxPyPz(lpxyz);
257     //
258     Float_t dz   = param->GetZ()-ltrp->GetZ();
259     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
260     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
261     if (hisdz) hisdz->Fill(dz);
262     if (hisdphi) hisdphi->Fill(dphi);
263     if (hisdphiP) hisdphiP->Fill(dphiP);
264     if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
265   }
266 }
267
268 void AliTPCcalibLaser::FitDriftV(){
269   //
270   // Fit drift velocity - linear approximation in the z and global y
271   //
272   static TLinearFitter fdriftA(3,"hyp2");
273   static TLinearFitter fdriftC(3,"hyp2");
274   fdriftA.ClearPoints();
275   fdriftC.ClearPoints();
276   //
277   for (Int_t id=0; id<336; id++){
278     if (!fTracksEsdParam.At(id)) continue;
279     if (!AcceptLaser(id)) continue;
280     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
281     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
282     Double_t xyz[3];
283     Double_t pxyz[3];
284     Double_t lxyz[3];
285     Double_t lpxyz[3];
286     param->GetXYZ(xyz);
287     param->GetPxPyPz(pxyz);
288     ltrp->GetXYZ(lxyz);
289     ltrp->GetPxPyPz(lpxyz);
290     Double_t xxx[2] = {lxyz[2],lxyz[1]};
291     if (ltrp->GetSide()==0){
292       fdriftA.AddPoint(xxx,xyz[2],1);
293     }else{
294       fdriftC.AddPoint(xxx,xyz[2],1);
295     }
296   }
297   Float_t chi2A = 0;
298   Float_t chi2C = 0;
299   Int_t npointsA=0;
300   Int_t npointsC=0;
301   //
302   if (fdriftA.GetNpoints()>10){
303     fdriftA.Eval();
304     fdriftA.EvalRobust(0.8);
305     fdriftA.GetParameters(*fFitAside);
306     npointsA= fdriftA.GetNpoints();
307     chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
308   }
309   if (fdriftC.GetNpoints()>10){
310     fdriftC.Eval();
311     fdriftC.EvalRobust(0.8);
312     fdriftC.GetParameters(*fFitCside);
313     npointsC= fdriftC.GetNpoints();
314     chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
315   }
316
317   if (fStreamLevel>0){
318     TTreeSRedirector *cstream = GetDebugStreamer();
319     Int_t time = fESD->GetTimeStamp();
320     if (cstream){
321       (*cstream)<<"driftv"<<
322         "driftA.="<<fFitAside<<
323         "driftC.="<<fFitCside<<
324         "chi2A="<<chi2A<<
325         "chi2C="<<chi2C<<
326         "nA="<<npointsA<<
327         "nC="<<npointsC<<
328         "time="<<time<<
329         "\n";
330     }
331   }
332   //
333 }
334
335
336 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
337   //
338   //
339   //
340   /*
341   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
342   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
343   TCut cutN("cutN","fTPCncls>70");
344   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
345   TCut cutA = cutT+cutPt+cutP;
346   */
347   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
348   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
349   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
350
351   if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
352   if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;
353   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
354   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
355   //
356   // dedx cut
357   //
358   if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
359   if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
360   //
361   return kTRUE;
362 }
363
364 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
365   //
366   // Find corresponding mirror
367   // add the corresponding tracks
368   //
369   Float_t kRadius0 = 252;
370   Float_t kRadius  = 253.4;
371   if (!track->GetOuterParam()) return -1;
372   AliExternalTrackParam param(*(track->GetOuterParam()));
373   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
374   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
375   AliTPCLaserTrack ltr;
376   AliTPCLaserTrack *ltrp=0x0;
377   AliTPCLaserTrack *ltrpjw=0x0;
378   //
379   Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param);
380  // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
381   //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
382
383   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
384     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
385   else
386     ltrp=&ltr;
387   /*
388     if (idjw!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw)))
389     ltrpjw=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(idjw);
390   else
391     ltrpjw=&ltr;
392
393
394     if (fStreamLevel>0){
395     TTreeSRedirector *cstream = GetDebugStreamer();
396     if (cstream){
397         (*cstream)<<"idcmp"<<
398             "id=" << id <<
399             "idjw=" << idjw <<
400             "tr.="  << ltrp <<
401             "trjw.="<< ltrpjw <<
402             "seed.="<<seed<<
403             "event="<<fEvent <<
404             "\n";
405     }
406   }
407
408       */
409
410
411   if (id>=0){
412     //
413     //
414     Float_t radius=TMath::Abs(ltrp->GetX());
415     AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
416     //
417     if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
418     fTracksEsdParam.AddAt(param.Clone(),id);
419     fTracksEsd.AddAt(track,id);
420     fTracksTPC.AddAt(seed,id);
421     //
422   }
423
424   return id;
425 }
426
427
428
429 void AliTPCcalibLaser::DumpLaser(Int_t id) {
430   //
431   //  Dump Laser info to the tree
432   //
433   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
434   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
435   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
436   //
437   // Fast laser ID
438   //
439   Double_t xyz[3];
440   Double_t pxyz[3];
441   Double_t lxyz[3];
442   Double_t lpxyz[3];
443   param->GetXYZ(xyz);
444   param->GetPxPyPz(pxyz);
445   ltrp->GetXYZ(lxyz);
446   ltrp->GetPxPyPz(lpxyz);
447
448   if (fStreamLevel>0){
449     TTreeSRedirector *cstream = GetDebugStreamer();
450     Int_t time = fESD->GetTimeStamp();
451     Bool_t accept = AcceptLaser(id);
452     if (cstream){
453       (*cstream)<<"Track"<<
454         "run="<<fRun<<
455         "id="<<id<<
456         "accept="<<accept<<
457         "driftA.="<<fFitAside<<
458         "driftC.="<<fFitCside<<
459         "time="<<time<<
460         //
461         "LTr.="<<ltrp<<
462         "Esd.="<<track<<
463         "Tr.="<<param<<
464         "x0="<<xyz[0]<<
465         "x1="<<xyz[1]<<
466         "x2="<<xyz[2]<<
467         "px0="<<pxyz[0]<<
468         "px1="<<pxyz[1]<<
469         "px2="<<pxyz[2]<<
470         //
471         "lx0="<<lxyz[0]<<
472         "lx1="<<lxyz[1]<<
473         "lx2="<<lxyz[2]<<
474         "lpx0="<<lpxyz[0]<<
475         "lpx1="<<lpxyz[1]<<
476         "lpx2="<<lpxyz[2]<<
477         "\n";
478     }
479   }
480 }
481
482 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
483   //
484   // Refit the track with different tracklet models:
485   // 1. Per ROC using the kalman filter, different edge cuts
486   // 2. Per ROC linear in y and z
487   // 3. Per ROC quadratic in y and z
488   // 4. Per track offset for each sector, linear for each sector, common quadratic
489   // store x, y, z information for all models and the cluster to calculate the residuals
490   //
491   AliTPCseed *track      = (AliTPCseed*)fTracksTPC.At(id);
492   AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
493   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
494
495   AliTPCclusterMI dummyCl;
496
497   //two tracklets
498   Int_t kMaxTracklets=2;
499   //=============================================//
500   // Linear Fitters for the Different Approaches //
501   //=============================================//
502   //linear fit model in y and z; inner - outer sector
503   static TLinearFitter fy1I(2,"hyp1");
504   static TLinearFitter fy1O(2,"hyp1");
505   static TLinearFitter fz1I(2,"hyp1");
506   static TLinearFitter fz1O(2,"hyp1");
507   //quadratic fit model in y and z; inner - sector
508   static TLinearFitter fy2I(3,"hyp2");
509   static TLinearFitter fy2O(3,"hyp2");
510   static TLinearFitter fz2I(3,"hyp2");
511   static TLinearFitter fz2O(3,"hyp2");
512   //common quadratic fit for IROC and OROC in y and z
513   static TLinearFitter fy4(5,"hyp4");
514   static TLinearFitter fz4(5,"hyp4");
515
516
517   //set standard cuts
518   if ( fNcuts==0 ){
519       fNcuts=1;
520       fEdgeXcuts[0]=4;
521       fEdgeYcuts[0]=3;
522       fNClCuts[0]=20;
523   }
524   //=============================//
525   // Loop over all Tracklet Cuts //
526   //=============================//
527   for (Int_t icut=0; icut<fNcuts; icut++){
528       AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
529       //cut parameters
530       Double_t edgeCutX = fEdgeXcuts[icut];
531       Double_t edgeCutY = fEdgeYcuts[icut];
532       Int_t    nclCut   = fNClCuts[icut];
533       //=========================//
534       // Parameters to calculate //
535       //=========================//
536       //fit parameter inner, outer tracklet and combined fit
537       TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
538       TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
539       //
540       TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
541       TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
542       TVectorD vecy4res(5),vecz4res(5);
543       // cluster and track positions for each row - used for residuals
544       TVectorD vecX(159);        // x is the same for all (row center)
545       TVectorD vecYkalman(159);  // y from kalman fit
546       TVectorD vecZkalman(159);  // z from kalman fit
547       TVectorD vecY1(159);       // y from pol1 fit per ROC
548       TVectorD vecZ1(159);       // z from pol1 fit per ROC
549       TVectorD vecY2(159);       // y from pol2 fit per ROC
550       TVectorD vecZ2(159);       // z from pol2 fit per ROC
551       TVectorD vecY4(159);       // y from sector fit
552       TVectorD vecZ4(159);       // z from sector fit
553       TVectorD vecClY(159);      // y cluster position
554       TVectorD vecClZ(159);      // z cluster position
555       TVectorD vecSec(159);      // sector for each row
556       //chi2 of fits
557       Double_t chi2I1z=-1;       // chi2 of pol1 fit in z (inner)
558       Double_t chi2I1y=-1;       // chi2 of pol1 fit in y (inner)
559       Double_t chi2O1z=-1;       // chi2 of pol1 fit in z (outer)
560       Double_t chi2O1y=-1;       // chi2 of pol1 fit in y (outer)
561       Double_t chi2I2z=-1;       // chi2 of pol2 fit in z (inner)
562       Double_t chi2I2y=-1;       // chi2 of pol2 fit in y (inner)
563       Double_t chi2O2z=-1;       // chi2 of pol2 fit in z (outer)
564       Double_t chi2O2y=-1;       // chi2 of pol2 fit in y (outer)
565       Double_t chi2IOz=-1;       // chi2 of hyp4 fit in z (inner+outer)
566       Double_t chi2IOy=-1;       // chi2 of hyp4 fit in y (inner+outer)
567       //more
568       Int_t innerSector = -1;    // number of inner sector
569       Int_t outerSector = -1;    // number of outer sector
570       Int_t nclI=0;              // number of clusters (inner)
571       Int_t nclO=0;              // number of clusters (outer)
572       //=================================================//
573       // Perform the Kalman Fit using the Tracklet Class //
574       //=================================================//
575       AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
576       TObjArray tracklets=
577           AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
578                                           kFALSE,nclCut,kMaxTracklets);
579       // tracklet pointers
580       AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
581       AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
582       AliTPCTracklet *tr=0x0;
583       AliTPCTracklet dummy;
584       //continue if we didn't find a tracklet
585       if ( !trInner && !trOuter ) continue;
586       //================================================================================//
587       // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
588       //================================================================================//
589       if ( trInner && trOuter ){
590           if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
591           if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
592               tr = trInner;
593               trInner=trOuter;
594               trOuter=tr;
595               AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
596           }
597       } else {
598           if ( trInner ){
599               if ( !trInner->GetInner() ) continue;
600               trOuter=&dummy;
601               if ( trInner->GetSector()>35 ){
602                   trOuter=trInner;
603                   trInner=&dummy;
604               }
605           } else { //trOuter
606               if ( !trOuter->GetInner() ) continue;
607               trInner=&dummy;
608               if ( trOuter->GetSector()<36 ){
609                   trInner=trOuter;
610                   trOuter=&dummy;
611               }
612           }
613       }
614       innerSector = trInner->GetSector();
615       if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
616       outerSector = trOuter->GetSector();
617       if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
618
619       // array of clusters
620       TClonesArray arrCl("AliTPCclusterMI",159);
621       arrCl.ExpandCreateFast(159);
622       //=======================================//
623       // fill fitters with cluster information //
624       //=======================================//
625       AliDebug(3,"Filing Cluster Information");
626       for (Int_t irow=158;irow>-1;--irow) {
627           AliTPCclusterMI *c=track->GetClusterPointer(irow);
628           AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
629           cl=dummyCl;
630           //
631           vecSec[irow]=-1;
632           if (!c) continue;
633           //
634           Int_t roc = static_cast<Int_t>(c->GetDetector());
635           if ( roc!=innerSector && roc!=outerSector ) continue;
636           vecSec[irow]=roc;
637           //store clusters in clones array
638           cl=*c;
639           //cluster position
640           vecX[irow]   = c->GetX();
641           vecClY[irow] = c->GetY();
642           vecClZ[irow] = c->GetZ();
643           //
644           Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
645           Double_t y = vecClY[irow];
646           Double_t z = vecClZ[irow];
647           //
648           Double_t x2[2]={x,x*x};   //linear and parabolic parameters
649           Double_t x4[4]={0,0,0,0}; //hyp4 parameters
650           if ( roc == innerSector ) {
651               x4[0]=1; //offset inner - outer sector
652               x4[1]=x; //slope parameter inner sector
653           } else {
654               x4[2]=x; //slope parameter outer sector
655           }
656           x4[3]=x*x;   //common parabolic shape
657           //
658           if ( roc==innerSector ){
659               fy1I.AddPoint(x2,y);
660               fz1I.AddPoint(x2,z);
661               fy2I.AddPoint(x2,y);
662               fz2I.AddPoint(x2,z);
663               ++nclI;
664           }
665           if ( roc==outerSector ){
666               fy1O.AddPoint(x2,y);
667               fz1O.AddPoint(x2,z);
668               fy2O.AddPoint(x2,y);
669               fz2O.AddPoint(x2,z);
670               ++nclO;
671           }
672           fy4.AddPoint(x4,y);
673           fz4.AddPoint(x4,z);
674       }
675       //======================================//
676       // Evaluate and retrieve fit parameters //
677       //======================================//
678       AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
679       //inner sector
680       if (  (innerSector>-1) && (fy1I.GetNpoints()>0) ){
681           fy1I.Eval();
682           fz1I.Eval();
683           fy2I.Eval();
684           fz2I.Eval();
685           fy1I.GetParameters(vecy1resInner);
686           fz1I.GetParameters(vecz1resInner);
687           fy2I.GetParameters(vecy2resInner);
688           fz2I.GetParameters(vecz2resInner);
689           chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
690           chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
691           chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
692           chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
693       }
694       //outer sector
695       if (  (outerSector>-1) && (fy1O.GetNpoints()>0) ){
696           fy1O.Eval();
697           fz1O.Eval();
698           fy2O.Eval();
699           fz2O.Eval();
700           fy1O.GetParameters(vecy1resOuter);
701           fz1O.GetParameters(vecz1resOuter);
702           fy2O.GetParameters(vecy2resOuter);
703           fz2O.GetParameters(vecz2resOuter);
704           chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
705           chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
706           chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
707           chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
708       }
709       //combined hyp4 fit
710       if ( innerSector>0 && outerSector>0 ){
711           if (fy4.GetNpoints()>0) {
712               fy4.Eval();
713               fy4.GetParameters(vecy4res);
714               chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
715           }
716           if (fz4.GetNpoints()>0) {
717               fz4.Eval();
718               fz4.GetParameters(vecz4res);
719               chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
720           }
721       }
722       //clear points
723       fy4.ClearPoints();  fz4.ClearPoints();
724       fy1I.ClearPoints(); fy1O.ClearPoints();
725       fz1I.ClearPoints(); fz1O.ClearPoints();
726       fy2I.ClearPoints(); fy2O.ClearPoints();
727       fz2I.ClearPoints(); fz2O.ClearPoints();
728       //==============================//
729       // calculate tracklet positions //
730       //==============================//
731       AliDebug(4,"Calculate tracklet positions");
732       for (Int_t irow=158;irow>-1;--irow) {
733           if ( vecSec[irow]==-1 ) continue;  //no cluster info
734           if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
735           tr=&dummy;
736           Double_t x=vecX[irow];
737           Double_t xref=x-133.4;
738           //
739           Double_t yoffInner=0;
740           Double_t zoffInner=0;
741           Double_t yslopeInner=0;
742           Double_t yslopeOuter=0;
743           Double_t zslopeInner=0;
744           Double_t zslopeOuter=0;
745           //positions of hyperplane fits
746           if ( vecSec[irow] == outerSector ) {
747               tr=trOuter;
748               vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
749               vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
750               vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
751               vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
752               yslopeOuter=vecy4res[3];
753               zslopeOuter=vecz4res[3];
754           } else {
755               tr=trInner;
756               vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
757               vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
758               vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
759               vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
760               yoffInner=vecy4res[1];
761               zoffInner=vecz4res[1];
762               yslopeInner=vecy4res[2];
763               zslopeInner=vecz4res[2];
764           }
765           vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
766           vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
767           //positions of kalman fits
768           Double_t gxyz[3],xyz[3];
769           AliExternalTrackParam *param = 0x0;
770           //
771           param=tr->GetInner();
772           if (param){
773               param->GetXYZ(gxyz);
774               Float_t bz = AliTracker::GetBz(gxyz);
775               param->GetYAt(x, bz, xyz[1]);
776               param->GetZAt(x, bz, xyz[2]);
777               vecYkalman[irow]=xyz[1];
778               vecZkalman[irow]=xyz[2];
779           }
780       }
781       //=====================================================================//
782       // write results from the different tracklet fits with debug streamers //
783       //=====================================================================//
784       if (fStreamLevel>4){
785           TTreeSRedirector *cstream = GetDebugStreamer();
786           if (cstream){
787               Float_t dedx = track->GetdEdx();
788               (*cstream)<<"FitModels"<<
789                   "cutNr="      << icut <<
790                   "edgeCutX="   << edgeCutX <<
791                   "edgeCutY="   << edgeCutY <<
792                   "nclCut="     << nclCut <<
793                   "innerSector="<< innerSector <<
794                   "outerSector="<< outerSector <<
795                   "dEdx="       << dedx <<
796                   "LTr.="       << ltrp <<
797                   "Tr.="        << extparam <<
798                   "yPol1In.="   << &vecy1resInner <<
799                   "zPol1In.="   << &vecz1resInner <<
800                   "yPol2In.="   << &vecy2resInner <<
801                   "zPol2In.="   << &vecz2resInner <<
802                   "yPol1Out.="  << &vecy1resOuter <<
803                   "zPol1Out.="  << &vecz1resOuter <<
804                   "yPol2Out.="  << &vecy2resOuter <<
805                   "zPol2Out.="  << &vecz2resOuter <<
806                   "yInOut.="    << &vecy4res <<
807                   "zInOut.="    << &vecz4res <<
808                   "chi2y1In="   << chi2I1y <<
809                   "chi2z1In="   << chi2I1z <<
810                   "chi2y1Out="  << chi2O1y <<
811                   "chi2z1Out="  << chi2O1z <<
812                   "chi2y2In="   << chi2I2y <<
813                   "chi2z2In="   << chi2I2z <<
814                   "chi2y2Out="  << chi2O2y <<
815                   "chi2z2Out="  << chi2O2z <<
816                   "chi2yInOut=" << chi2IOy <<
817                   "chi2zInOut=" << chi2IOz <<
818                   "trletIn.="   << trInner <<
819                   "trletOut.="  << trOuter <<
820                   "nclI="       << nclI <<
821                   "nclO="       << nclO <<
822                   "\n";
823           }
824       }
825
826       // wirte residuals information
827       if (fStreamLevel>5){
828           TTreeSRedirector *cstream = GetDebugStreamer();
829           if (cstream){
830               Float_t dedx = track->GetdEdx();
831               (*cstream)<<"Residuals"<<
832                   "cutNr="      << icut <<
833                   "edgeCutX="   << edgeCutX <<
834                   "edgeCutY="   << edgeCutY <<
835                   "nclCut="     << nclCut   <<
836                   "LTr.="       << ltrp <<
837                   "Tr.="        << extparam<<
838                   "dEdx="       << dedx <<
839                   "Cl.="        << &arrCl <<
840                   "TrX.="       << &vecX <<
841                   "TrYpol1.="   << &vecY1 <<
842                   "TrZpol1.="   << &vecZ1 <<
843                   "TrYpol2.="   << &vecY2 <<
844                   "TrZpol2.="   << &vecZ2 <<
845                   "TrYInOut.="  << &vecY4 <<
846                   "TrZInOut.="  << &vecZ4 <<
847                   "ClY.="       << &vecClY <<
848                   "ClZ.="       << &vecClZ <<
849                   "sec.="       << &vecSec <<
850                   "nclI="       << nclI <<
851                   "nclO="       << nclO <<
852                   "yInOut.="    << &vecy4res <<
853                   "zInOut.="    << &vecz4res <<
854                   "chi2y1In="   << chi2I1y <<
855                   "chi2z1In="   << chi2I1z <<
856                   "chi2y1Out="  << chi2O1y <<
857                   "chi2z1Out="  << chi2O1z <<
858                   "chi2y2In="   << chi2I2y <<
859                   "chi2z2In="   << chi2I2z <<
860                   "chi2y2Out="  << chi2O2y <<
861                   "chi2z2Out="  << chi2O2z <<
862                   "chi2yInOut=" << chi2IOy <<
863                   "chi2zInOut=" << chi2IOz <<
864                   "\n";
865
866           }
867       }
868   }
869  /*
870
871   Int_t indexMaxCut[kMaxTracklets];
872
873   Float_t xMinCut[kMaxTracklets];
874   Float_t xMedCut[kMaxTracklets];
875   Float_t xMaxCut[kMaxTracklets];
876
877   for (Int_t i=0; i<kMaxTracklets; i++){
878       trMinCut = (AliTPCTracklet*)trackletsMinCuts->At(i);
879       trMedCut = (AliTPCTracklet*)trackletsMedCuts->At(i);
880       trMaxCut = (AliTPCTracklet*)trackletsMaxCuts->At(i);
881       if (!trMinCut ) trMinCut=&dummy;
882       if (!trMedCut ) trMedCut=&dummy;
883       if (!trMaxCut ) trMaxCut=&dummy;
884       xMinCut[i]=trMinCut->GetInner()->GetX();
885       xMedCut[i]=trMedCut->GetInner()->GetX();
886       xMaxCut[i]=trMaxCut->GetInner()->GetX();
887   }
888   TMath::Sort(kMaxTracklets, xMinCut, indexMinCut);
889   TMath::Sort(kMaxTracklets, xMedCut, indexMedCut);
890   TMath::Sort(kMaxTracklets, xMaxCut, indexMaxCut);
891    */
892
893 }
894
895
896 void AliTPCcalibLaser::RefitLaser(Int_t id){
897   //
898   // Refit the track store residuals
899   //
900
901   AliTPCseed *track    = (AliTPCseed*)fTracksTPC.At(id);
902   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
903   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
904
905   //linear fit model in y and z per sector
906   static TLinearFitter fy1(2,"hyp1");
907   static TLinearFitter fz1(2,"hyp1");
908   //quadratic fit model in y and z per sector
909   static TLinearFitter fy2(3,"hyp2");
910   static TLinearFitter fz2(3,"hyp2");
911   static TVectorD vecy2,vecz2,vecy1,vecz1;
912
913   const Int_t kMinClusters=20;
914   Int_t nclusters[72];
915   //
916   for (Int_t i=0;i<72;++i) nclusters[i]=0;
917
918   for (Int_t i=0;i<160;++i) {
919     AliTPCclusterMI *c=track->GetClusterPointer(i);
920     if (c) nclusters[c->GetDetector()]++;
921   }
922
923   for (Int_t isec=0; isec<72;isec++){
924     if (nclusters[isec]<kMinClusters) continue;
925     fy2.ClearPoints();
926     fz2.ClearPoints();
927     fy1.ClearPoints();
928     fz1.ClearPoints();
929     //
930     for (Int_t irow=0;irow<160;++irow) {
931       AliTPCclusterMI *c=track->GetClusterPointer(irow);
932       //if (c && RejectCluster(c)) continue;
933       Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
934       if (c&&c->GetDetector()==isec) {
935         Double_t x[2]={xd,xd*xd};
936         fy2.AddPoint(x,c->GetY());
937         fz2.AddPoint(x,c->GetZ());
938         //
939         fy1.AddPoint(x,c->GetY());
940         fz1.AddPoint(x,c->GetZ());
941       }
942     }
943     fy2.Eval();
944     fz2.Eval();
945     fy1.Eval();
946     fz1.Eval();
947     fy1.GetParameters(vecy1);
948     fy2.GetParameters(vecy2);
949     fz1.GetParameters(vecz1);
950     fz2.GetParameters(vecz2);
951
952     if (fStreamLevel>0){
953       TTreeSRedirector *cstream = GetDebugStreamer();
954       if (cstream){
955         Float_t dedx = track->GetdEdx();
956         (*cstream)<<"Tracklet"<<
957           "LTr.="<<ltrp<<
958           "Tr.="<<param<<
959           "sec="<<isec<<
960           "ncl="<<nclusters[isec]<<
961           "dedx="<<dedx<<
962           "dedx="<<dedx<<
963           "vy1.="<<&vecy1<<
964           "vy2.="<<&vecy2<<
965           "vz1.="<<&vecz1<<
966           "vz2.="<<&vecz2<<
967           "\n";
968       }
969     }
970   }
971   //
972   //
973   //
974   //   for (Int_t irow=0;irow<160;++irow) {
975   //       AliTPCclusterMI *c=track->GetClusterPointer(irow);
976   //       if (c && RejectCluster(c)) continue;
977   //       if (c&&c->GetDetector()==isec) {
978   //    Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
979   //    fy2.AddPoint(&x,c->GetY());
980   //    fz2.AddPoint(&x,c->GetZ());
981   //    //
982   //    fy1.AddPoint(&x,c->GetY());
983   //    fz1.AddPoint(&x,c->GetZ());
984   //       }
985   //     }
986
987 }
988
989
990 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
991   //
992   //  Dump information about laser beams
993   //  isOK variable indicates usability of the beam  
994   //  Beam is not usable if:
995   //     a.  No entries in range (krmsCut0)
996   //     b.  Big sperad          (krmscut1)
997   //     c.  RMSto fit sigma bigger then (kmultiCut)
998   //     d.  Too big angular spread 
999   //  
1000
1001   const Float_t krmsCut0=0.001;
1002   const Float_t krmsCut1=0.16;
1003   const Float_t kmultiCut=2;
1004   const Float_t kcutP0=0.002;
1005   //
1006   AliTPCcalibLaser *laser = this;
1007   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1008   TF1 fg("fg","gaus");
1009   //
1010   //
1011   for (Int_t id=0; id<336; id++){
1012     Bool_t isOK=kTRUE;
1013     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
1014     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1015     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
1016     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
1017     if (!hisphi) continue;;
1018     Double_t entries = hisphi->GetEntries();
1019     if (entries<minEntries) continue;
1020     //
1021     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1022     if (!ltrp) {
1023      AliTPCLaserTrack::LoadTracks();
1024       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1025     }
1026     Float_t meanphi = hisphi->GetMean();
1027     Float_t rmsphi = hisphi->GetRMS();
1028     //
1029     Float_t meanphiP = hisphiP->GetMean();
1030     Float_t rmsphiP = hisphiP->GetRMS();
1031     Float_t meanZ = hisZ->GetMean();
1032     Float_t rmsZ = hisZ->GetRMS();
1033     hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1034     Double_t gphi1 = fg.GetParameter(1);
1035     Double_t gphi2 = fg.GetParameter(2);
1036     hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1037     Double_t gphiP1 = fg.GetParameter(1);
1038     Double_t gphiP2 = fg.GetParameter(2);
1039     hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
1040     Double_t gz1 = fg.GetParameter(1);
1041     Double_t gz2 = fg.GetParameter(2);
1042     //
1043     Float_t meanS=hisS->GetMean();
1044     //
1045     Double_t lxyz[3];
1046     Double_t lpxyz[3];
1047     ltrp->GetXYZ(lxyz);
1048     ltrp->GetPxPyPz(lpxyz);
1049
1050     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1051     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1052     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
1053     if (gphiP2>kcutP0) isOK=kFALSE;
1054     //
1055     (*pcstream)<<"Mean"<<
1056       "isOK="<<isOK<<
1057       "entries="<<entries<<      // number of entries
1058       "bz="<<bfield<<            // bfield
1059       "LTr.="<<ltrp<<             // refernece track
1060       //
1061       "lx0="<<lxyz[0]<<          // reference x
1062       "lx1="<<lxyz[1]<<          // reference y
1063       "lx2="<<lxyz[2]<<          // refernece z
1064       "lpx0="<<lpxyz[0]<<          // reference x
1065       "lpx1="<<lpxyz[1]<<          // reference y
1066       "lpx2="<<lpxyz[2]<<          // refernece z
1067       //
1068       "msig="<<meanS<<
1069       //
1070       "mphi="<<meanphi<<         //
1071       "rmsphi="<<rmsphi<<        //
1072       "gphi1="<<gphi1<<
1073       "gphi2="<<gphi2<<
1074       //
1075       "mphiP="<<meanphiP<<       //
1076       "rmsphiP="<<rmsphiP<<      //
1077       "gphiP1="<<gphiP1<<
1078       "gphiP2="<<gphiP2<<
1079       //
1080       "meanZ="<<meanZ<<
1081       "rmsZ="<<rmsZ<<
1082       "gz1="<<gz1<<
1083       "gz2="<<gz2<<
1084
1085       "\n";
1086   }
1087   delete pcstream;
1088 }
1089
1090
1091
1092 void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1093   //
1094   //
1095   //
1096   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1097   TFile * f = pcstream->GetFile();
1098   f->mkdir("dirphi");
1099   f->mkdir("dirphiP");
1100   f->mkdir("dirZ");
1101   TF1 fp("p1","pol1");
1102   //
1103   //
1104   char    cut[1000];
1105   char grname[1000];
1106   char grnamefull[1000];
1107
1108   Double_t mphi[100];
1109   Double_t mphiP[100];
1110   Double_t smphi[100];
1111   Double_t smphiP[100];
1112   Double_t mZ[100];
1113   Double_t smZ[100];
1114   Double_t bz[100];
1115   Double_t sbz[100];
1116   // fit parameters
1117   Double_t pphi[3];
1118   Double_t pphiP[3];
1119   Double_t pmZ[3];
1120   //
1121   for (Int_t id=0; id<336; id++){
1122     // id =205;
1123     sprintf(cut,"isOK&&fId==%d",id);
1124     Int_t entries = chain->Draw("bz",cut,"goff");
1125     if (entries<3) continue;
1126     AliTPCLaserTrack *ltrp = 0;;
1127     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1128     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1129     Double_t lxyz[3];
1130     Double_t lpxyz[3];
1131     ltrp->GetXYZ(lxyz);
1132     ltrp->GetPxPyPz(lpxyz);
1133
1134     chain->Draw("bz",cut,"goff");
1135     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1136     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1137     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1138     //
1139     chain->Draw("gphi1",cut,"goff");
1140     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1141     chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1142     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1143     //
1144     chain->Draw("gphiP1",cut,"goff");
1145     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1146     chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1147     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1148     //
1149     chain->Draw("gz1",cut,"goff");
1150     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1151     chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1152     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1153     //
1154     //
1155     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1156             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1157     // store data  
1158     // phi
1159     f->cd("dirphi");
1160     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1161     grphi->Draw("a*");
1162     grphi->Fit(&fp);
1163     pphi[0] = fp.GetParameter(0);                          // offset
1164     pphi[1] = fp.GetParameter(1);                          // slope
1165     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1166     sprintf(grname,"phi_id%d",id);
1167     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
1168     grphi->GetXaxis()->SetTitle("b_{z} (T)");
1169     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
1170     grphi->SetMaximum(1.2);
1171     grphi->SetMinimum(-1.2);
1172     grphi->Draw("a*");
1173
1174     grphi->Write();
1175     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1176     // phiP
1177     f->cd("dirphiP)");
1178     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1179     grphiP->Draw("a*");
1180     grphiP->Fit(&fp);
1181     pphiP[0] = fp.GetParameter(0);                          // offset
1182     pphiP[1] = fp.GetParameter(1);                          // slope
1183     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1184     sprintf(grname,"phiP_id%d",id);
1185     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
1186     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1187     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
1188     grphiP->SetMaximum(pphiP[0]+0.005);
1189     grphiP->SetMinimum(pphiP[0]-0.005);
1190
1191     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1192     grphiP->Write();
1193     //
1194     //Z
1195     f->cd("dirZ");
1196     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1197     grmZ->Draw("a*");
1198     grmZ->Fit(&fp);
1199     pmZ[0] = fp.GetParameter(0);                          // offset
1200     pmZ[1] = fp.GetParameter(1);                          // slope
1201     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1202     sprintf(grname,"mZ_id%d",id);
1203     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
1204     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1205     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1206
1207     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1208     grmZ->Write();
1209     
1210
1211     for (Int_t ientry=0; ientry<entries; ientry++){
1212       (*pcstream)<<"Mean"<<
1213         "id="<<id<<
1214         "LTr.="<<ltrp<<
1215         "entries="<<entries<<
1216         "bz="<<bz[ientry]<<
1217         "lx0="<<lxyz[0]<<          // reference x
1218         "lx1="<<lxyz[1]<<          // reference y
1219         "lx2="<<lxyz[2]<<          // refernece z      
1220         "lpx0="<<lpxyz[0]<<          // reference x
1221         "lpx1="<<lpxyz[1]<<          // reference y
1222         "lpx2="<<lpxyz[2]<<          // refernece z            
1223         //values
1224         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1225         "pphi0="<<pphi[0]<<   // offset
1226         "pphi1="<<pphi[1]<<   // mean
1227         "pphi2="<<pphi[2]<<   // norm chi2
1228         //
1229         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1230         "pphiP0="<<pphiP[0]<< // offset
1231         "pphiP1="<<pphiP[1]<< // mean
1232         "pphiP2="<<pphiP[2]<< // norm chi2
1233         //
1234         "gz1="<<mZ[ientry]<<
1235         "pmZ0="<<pmZ[0]<<     // offset
1236         "pmZ1="<<pmZ[1]<<     // mean
1237         "pmZ2="<<pmZ[2]<<     // norm chi2
1238         "\n";
1239     }
1240   }
1241   
1242   delete pcstream;
1243   
1244 }
1245
1246
1247 void AliTPCcalibLaser::Analyze(){
1248   //
1249   //
1250   //
1251 }
1252
1253
1254 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
1255
1256   TIterator* iter = li->MakeIterator();
1257   AliTPCcalibLaser* cal = 0;
1258
1259   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1260     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1261       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1262       return -1;
1263     }
1264     //
1265    //  fHistNTracks->Add(cal->fHistNTracks);
1266 //     fClusters->Add(cal->fClusters);
1267 //     fModules->Add(cal->fModules);
1268 //     fHistPt->Add(cal->fHistPt);
1269 //     fPtResolution->Add(cal->fPtResolution);
1270 //     fDeDx->Add(cal->fDeDx);
1271
1272
1273     TH1F *h=0x0;
1274     TH1F *hm=0x0;
1275
1276     for (Int_t id=0; id<336; id++){
1277       // merge fDeltaZ histograms
1278       hm = (TH1F*)cal->fDeltaZ.At(id);
1279       h  = (TH1F*)fDeltaZ.At(id);
1280       if (!h) {
1281         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1282         fDeltaZ.AddAt(h,id);
1283       }
1284       if (hm) h->Add(hm);
1285       // merge fDeltaPhi histograms
1286       hm = (TH1F*)cal->fDeltaPhi.At(id);
1287       h  = (TH1F*)fDeltaPhi.At(id);
1288       if (!h) {
1289         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
1290         fDeltaPhi.AddAt(h,id);
1291       }
1292       if (hm) h->Add(hm);
1293       // merge fDeltaPhiP histograms
1294       hm = (TH1F*)cal->fDeltaPhiP.At(id);
1295       h  = (TH1F*)fDeltaPhiP.At(id);
1296       if (!h) {
1297         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
1298         fDeltaPhiP.AddAt(h,id);
1299       }
1300       if (hm) h->Add(hm);
1301       // merge fSignals histograms
1302       hm = (TH1F*)cal->fSignals.At(id);
1303       h  = (TH1F*)fSignals.At(id);
1304       if (!h) {
1305         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
1306         fSignals.AddAt(h,id);
1307       }
1308       if (hm) h->Add(hm);
1309     }
1310   }
1311   return 0;
1312 }
1313
1314
1315
1316
1317 /*
1318  gSystem->Load("libSTAT.so")
1319  TStatToolkit toolkit;
1320  Double_t chi2;
1321  TVectorD fitParam;
1322  TMatrixD covMatrix;
1323  Int_t npoints;
1324  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
1325
1326
1327 TString fstring="";
1328 //
1329 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
1330 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
1331 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
1332 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4 
1333 //
1334 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
1335 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
1336 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
1337 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
1338 //   
1339 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
1340 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
1341 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
1342 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
1343
1344
1345
1346
1347  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1348
1349  treeT->SetAlias("fit",strq0->Data());
1350  
1351
1352  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1353
1354  treeT->SetAlias("fitP",strqP->Data());
1355
1356
1357  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1358  treeT->SetAlias("fitD",strqDrift->Data());
1359
1360
1361 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,""); 
1362 {
1363 for (Int_t i=0; i<6;i++){
1364 treeT->SetLineColor(i+2);
1365 treeT->SetMarkerSize(1);
1366 treeT->SetMarkerStyle(22+i);
1367 treeT->SetMarkerColor(i+2);
1368
1369 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same"); 
1370 }
1371
1372  */
1373  
1374
1375
1376 /*
1377   TTree * tree = (TTree*)f.Get("FitModels");
1378
1379   TEventList listLFit0("listLFit0","listLFit0");
1380   TEventList listLFit1("listLFit1","listLFit1");
1381   
1382   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1383   tree->SetEventList(&listLFit0);
1384   
1385
1386 */