Adding fitting of B scanned values (Marian)
[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 "TClonesArray.h"
85 #include "TPad.h"
86
87
88 #include "TTreeStream.h"
89 #include <iostream>
90 #include <sstream>
91 #include "AliTPCLaserTrack.h"
92
93 using namespace std;
94
95 ClassImp(AliTPCcalibLaser)
96
97 AliTPCcalibLaser::AliTPCcalibLaser():
98   AliTPCcalibBase(),
99   fESD(0),
100   fESDfriend(),
101   fTracksMirror(336),
102   fTracksEsd(336),
103   fTracksEsdParam(336),
104   fTracksTPC(336),
105   fDeltaZ(336),          // array of histograms of delta z for each track
106   fDeltaPhi(336),          // array of histograms of delta z for each track
107   fDeltaPhiP(336),          // array of histograms of delta z for each track
108   fSignals(336),           // array of dedx signals
109   fFitAside(new TVectorD(3)),        // drift fit - A side
110   fFitCside(new TVectorD(3)),        // drift fit - C- side
111   fRun(0)
112 {
113   //
114   // Constructor
115   //
116   fTracksEsdParam.SetOwner(kTRUE);
117 }
118
119 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
120   AliTPCcalibBase(),
121   fESD(0),
122   fESDfriend(0),
123   fTracksMirror(336),
124   fTracksEsd(336),
125   fTracksEsdParam(336),
126   fTracksTPC(336),
127   fDeltaZ(336),          // array of histograms of delta z for each track
128   fDeltaPhi(336),          // array of histograms of delta z for each track
129   fDeltaPhiP(336),          // array of histograms of delta z for each track
130   fSignals(336),           // array of dedx signals
131   fFitAside(new TVectorD(3)),        // drift fit - A side
132   fFitCside(new TVectorD(3)),        // drift fit - C- side
133   fRun(0)
134 {
135   SetName(name);
136   SetTitle(title);
137   //
138   // Constructor
139   //
140   fTracksEsdParam.SetOwner(kTRUE);  
141 }
142
143 AliTPCcalibLaser::~AliTPCcalibLaser() {
144   //
145   // destructor
146   //
147 }
148
149
150
151 void AliTPCcalibLaser::Process(AliESDEvent * event) {
152   //
153   // 
154   // Loop over tracks and call  Process function
155   //
156   fESD = event;
157   if (!fESD) {
158     return;
159   }
160   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
161   if (!fESDfriend) {
162     return;
163   }
164   fTracksTPC.Clear();
165   fTracksEsd.Clear();
166   fTracksEsdParam.Delete();
167   //
168   Int_t n=fESD->GetNumberOfTracks();
169   Int_t run = fESD->GetRunNumber();
170   fRun = run;
171   for (Int_t i=0;i<n;++i) {
172     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
173     AliESDtrack *track=fESD->GetTrack(i);
174     TObject *calibObject=0;
175     AliTPCseed *seed=0;
176     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
177       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
178         break;
179     if (track&&seed) FindMirror(track,seed);
180     //
181   }
182   
183   FitDriftV();
184   MakeDistHisto();
185   //
186   for (Int_t id=0; id<336; id++){
187     //
188     //
189     if (!fTracksEsdParam.At(id)) continue;
190     DumpLaser(id);
191     RefitLaser(id);    
192     
193   }
194 }
195
196 void AliTPCcalibLaser::MakeDistHisto(){
197   //
198   //
199   //
200   for (Int_t id=0; id<336; id++){
201     //
202     //
203     if (!fTracksEsdParam.At(id)) continue;  
204     if (!AcceptLaser(id)) continue;
205     //
206     //
207     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
208     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
209     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
210     TH1F * hisSignal = (TH1F*)fSignals.At(id);
211
212     if (!hisdz){      
213       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
214       hisdz->SetDirectory(0);
215       fDeltaZ.AddAt(hisdz,id);
216       //
217       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
218       hisdphi->SetDirectory(0);
219       fDeltaPhi.AddAt(hisdphi,id);
220       //
221       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
222       hisdphiP->SetDirectory(0);
223       fDeltaPhiP.AddAt(hisdphiP,id);
224       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
225       hisSignal->SetDirectory(0);
226       fSignals.AddAt(hisSignal,id);
227     }
228
229     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
230     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
231     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
232     if (!param) return;
233     if (!ltrp) return;
234     if (!track) return;
235     Double_t xyz[3];
236     Double_t pxyz[3];
237     Double_t lxyz[3];
238     Double_t lpxyz[3];
239     param->GetXYZ(xyz);
240     param->GetPxPyPz(pxyz);
241     ltrp->GetXYZ(lxyz);
242     ltrp->GetPxPyPz(lpxyz);
243     //
244     Float_t dz   = param->GetZ()-ltrp->GetZ();
245     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
246     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
247     if (hisdz) hisdz->Fill(dz);
248     if (hisdphi) hisdphi->Fill(dphi);
249     if (hisdphiP) hisdphiP->Fill(dphiP); 
250     if (hisSignal) hisSignal->Fill(track->GetTPCsignal());
251   }
252 }
253
254 void AliTPCcalibLaser::FitDriftV(){
255   //
256   // Fit drift velocity - linear approximation in the z and global y
257   //
258   static TLinearFitter fdriftA(3,"hyp2");
259   static TLinearFitter fdriftC(3,"hyp2");
260   fdriftA.ClearPoints();
261   fdriftC.ClearPoints();
262   //
263   for (Int_t id=0; id<336; id++){
264     if (!fTracksEsdParam.At(id)) continue;  
265     if (!AcceptLaser(id)) continue;
266     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
267     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
268     Double_t xyz[3];
269     Double_t pxyz[3];
270     Double_t lxyz[3];
271     Double_t lpxyz[3];
272     param->GetXYZ(xyz);
273     param->GetPxPyPz(pxyz);
274     ltrp->GetXYZ(lxyz);
275     ltrp->GetPxPyPz(lpxyz);
276     Double_t xxx[2] = {lxyz[2],lxyz[1]};
277     if (ltrp->GetSide()==0){
278       fdriftA.AddPoint(xxx,xyz[2],1);
279     }else{
280       fdriftC.AddPoint(xxx,xyz[2],1);
281     }
282   }
283   Float_t chi2A = 0;
284   Float_t chi2C = 0;
285   Int_t npointsA=0;
286   Int_t npointsC=0;
287   //
288   if (fdriftA.GetNpoints()>10){
289     fdriftA.Eval();
290     fdriftA.EvalRobust(0.8);
291     fdriftA.GetParameters(*fFitAside);
292     npointsA= fdriftA.GetNpoints();
293     chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
294   }
295   if (fdriftC.GetNpoints()>10){
296     fdriftC.Eval();
297     fdriftC.EvalRobust(0.8);
298     fdriftC.GetParameters(*fFitCside);
299     npointsC= fdriftC.GetNpoints();
300     chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
301   }
302   
303   if (fStreamLevel>0){
304     TTreeSRedirector *cstream = GetDebugStreamer();
305     Int_t time = fESD->GetTimeStamp();
306     if (cstream){
307       (*cstream)<<"driftv"<<
308         "driftA.="<<fFitAside<<
309         "driftC.="<<fFitCside<<
310         "chi2A="<<chi2A<<
311         "chi2C="<<chi2C<<
312         "nA="<<npointsA<<
313         "nC="<<npointsC<<
314         "time="<<time<<
315         "\n";
316     }
317   }
318   //
319 }
320
321
322 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
323   //
324   //
325   //
326   /*
327   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
328   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
329   TCut cutN("cutN","fTPCncls>70");
330   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
331   TCut cutA = cutT+cutPt+cutP;
332   */
333   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
334   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
335   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
336
337   if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE;
338   if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE;  
339   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE;
340   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE;
341   //
342   // dedx cut
343   //
344   if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE;
345   if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE;
346   //
347   return kTRUE;  
348 }
349
350 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
351   //
352   // Find corresponding mirror
353   // add the corresponding tracks
354   //
355   Float_t kRadius0 = 252;
356   Float_t kRadius  = 253.4;
357   if (!track->GetOuterParam()) return -1;
358   AliExternalTrackParam param(*(track->GetOuterParam()));
359   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
360   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
361   AliTPCLaserTrack ltr;
362   AliTPCLaserTrack *ltrp=0x0;
363   //
364   Int_t id = AliTPCLaserTrack::IdentifyTrack(&param);
365   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id))) 
366     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
367   else 
368     ltrp=&ltr;
369   
370   if (id>=0){
371     //
372     //
373     Float_t radius=TMath::Abs(ltrp->GetX());
374     AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
375     //
376     if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
377     fTracksEsdParam.AddAt(param.Clone(),id);
378     fTracksEsd.AddAt(track,id);
379     fTracksTPC.AddAt(seed,id);
380     //
381   }
382
383   return id;
384 }
385
386
387
388 void AliTPCcalibLaser::DumpLaser(Int_t id) {
389   //
390   //  Dump Laser info to the tree
391   //
392   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
393   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
394   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
395   //
396   // Fast laser ID
397   //
398   Double_t xyz[3];
399   Double_t pxyz[3];
400   Double_t lxyz[3];
401   Double_t lpxyz[3];
402   param->GetXYZ(xyz);
403   param->GetPxPyPz(pxyz);
404   ltrp->GetXYZ(lxyz);
405   ltrp->GetPxPyPz(lpxyz);
406
407   if (fStreamLevel>0){
408     TTreeSRedirector *cstream = GetDebugStreamer();
409     Int_t time = fESD->GetTimeStamp();
410     Bool_t accept = AcceptLaser(id);
411     if (cstream){
412       (*cstream)<<"Track"<<
413         "run="<<fRun<<
414         "id="<<id<<
415         "accept="<<accept<<
416         "driftA.="<<fFitAside<<
417         "driftC.="<<fFitCside<<
418         "time="<<time<<
419         //
420         "LTr.="<<ltrp<<
421         "Esd.="<<track<<
422         "Tr.="<<param<<
423         "x0="<<xyz[0]<<
424         "x1="<<xyz[1]<<
425         "x2="<<xyz[2]<<
426         "px0="<<pxyz[0]<<
427         "px1="<<pxyz[1]<<
428         "px2="<<pxyz[2]<<
429         //
430         "lx0="<<lxyz[0]<<
431         "lx1="<<lxyz[1]<<
432         "lx2="<<lxyz[2]<<
433         "lpx0="<<lpxyz[0]<<
434         "lpx1="<<lpxyz[1]<<
435         "lpx2="<<lpxyz[2]<<
436         "\n";
437     }
438   }
439 }
440
441
442 void AliTPCcalibLaser::RefitLaser(Int_t id){
443   //
444   // Refit the track store residuals
445   //
446
447   AliTPCseed *track    = (AliTPCseed*)fTracksTPC.At(id);
448   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
449   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
450                              
451   //
452   static TLinearFitter fy2(3,"hyp2");
453   static TLinearFitter fz2(3,"hyp2");
454   static TLinearFitter fy1(2,"hyp1");
455   static TLinearFitter fz1(2,"hyp1");
456   static TVectorD vecy2,vecz2,vecy1,vecz1;
457
458   const Int_t kMinClusters=20;
459   Int_t nclusters[72]; 
460   //
461   for (Int_t i=0;i<72;++i) nclusters[i]=0;
462
463   for (Int_t i=0;i<160;++i) {    
464     AliTPCclusterMI *c=track->GetClusterPointer(i);
465     if (c) nclusters[c->GetDetector()]++;
466   }
467    
468   for (Int_t isec=0; isec<72;isec++){
469     if (nclusters[isec]<kMinClusters) continue;
470     fy2.ClearPoints();
471     fz2.ClearPoints();
472     fy1.ClearPoints();
473     fz1.ClearPoints();
474     //
475     for (Int_t irow=0;irow<160;++irow) {      
476       AliTPCclusterMI *c=track->GetClusterPointer(irow);
477       //if (c && RejectCluster(c)) continue;
478       if (c&&c->GetDetector()==isec) {
479         Double_t xd = c->GetX()-120;;
480         Double_t x[2]={xd,xd*xd};
481         fy2.AddPoint(x,c->GetY());
482         fz2.AddPoint(x,c->GetZ());
483         //
484         fy1.AddPoint(x,c->GetY());
485         fz1.AddPoint(x,c->GetZ());
486       }
487     }
488     fy2.Eval();
489     fz2.Eval();
490     fy1.Eval();
491     fz1.Eval();
492     fy1.GetParameters(vecy1);
493     fy2.GetParameters(vecy2);
494     fz1.GetParameters(vecz1);
495     fz2.GetParameters(vecz2);
496     
497     if (fStreamLevel>0){
498       TTreeSRedirector *cstream = GetDebugStreamer();
499       if (cstream){
500         Float_t dedx = track->GetdEdx();
501         (*cstream)<<"Tracklet"<<
502           "LTr.="<<ltrp<<
503           "Tr.="<<param<<
504           "sec="<<isec<<
505           "ncl="<<nclusters[isec]<<
506           "dedx="<<dedx<<
507           "dedx="<<dedx<<
508           "vy1.="<<&vecy1<<
509           "vy2.="<<&vecy2<<
510           "vz1.="<<&vecz1<<
511           "vz2.="<<&vecz2<<
512           "\n";
513       }
514     }
515   }
516   //
517   //
518   //
519   //   for (Int_t irow=0;irow<160;++irow) {      
520   //       AliTPCclusterMI *c=track->GetClusterPointer(irow);
521   //       if (c && RejectCluster(c)) continue;
522   //       if (c&&c->GetDetector()==isec) {
523   //    Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
524   //    fy2.AddPoint(&x,c->GetY());
525   //    fz2.AddPoint(&x,c->GetZ());
526   //    //
527   //    fy1.AddPoint(&x,c->GetY());
528   //    fz1.AddPoint(&x,c->GetZ());
529   //       }
530   //     }    
531   
532 }
533
534
535 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
536   //
537   //  Dump information about laser beams
538   //  isOK variable indicates usability of the beam  
539   //  Beam is not usable if:
540   //     a.  No entries in range (krmsCut0)
541   //     b.  Big sperad          (krmscut1)
542   //     c.  RMSto fit sigma bigger then (kmultiCut)
543   //     d.  Too big angular spread 
544   //  
545
546   const Float_t krmsCut0=0.001;
547   const Float_t krmsCut1=0.16;
548   const Float_t kmultiCut=2;
549   const Float_t kcutP0=0.002;
550   //
551   AliTPCcalibLaser *laser = this;
552   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
553   TF1 fg("fg","gaus");
554   //
555   //
556   for (Int_t id=0; id<336; id++){
557     Bool_t isOK=kTRUE;
558     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
559     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
560     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
561     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
562     if (!hisphi) continue;;
563     Double_t entries = hisphi->GetEntries();
564     if (entries<minEntries) continue;
565     //
566     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
567     if (!ltrp) {
568      AliTPCLaserTrack::LoadTracks();
569       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
570     }
571     Float_t meanphi = hisphi->GetMean();
572     Float_t rmsphi = hisphi->GetRMS();
573     //
574     Float_t meanphiP = hisphiP->GetMean();
575     Float_t rmsphiP = hisphiP->GetRMS();
576     Float_t meanZ = hisZ->GetMean();
577     Float_t rmsZ = hisZ->GetRMS();
578     hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
579     Double_t gphi1 = fg.GetParameter(1); 
580     Double_t gphi2 = fg.GetParameter(2); 
581     hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
582     Double_t gphiP1 = fg.GetParameter(1); 
583     Double_t gphiP2 = fg.GetParameter(2); 
584     hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
585     Double_t gz1 = fg.GetParameter(1); 
586     Double_t gz2 = fg.GetParameter(2); 
587     //
588     Float_t meanS=hisS->GetMean();
589     //
590     Double_t lxyz[3];
591     Double_t lpxyz[3];
592     ltrp->GetXYZ(lxyz);
593     ltrp->GetPxPyPz(lpxyz);
594
595     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
596     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
597     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
598     if (gphiP2>kcutP0) isOK=kFALSE;
599     //
600     (*pcstream)<<"Mean"<<
601       "isOK="<<isOK<<
602       "entries="<<entries<<      // number of entries
603       "bz="<<bfield<<            // bfield
604       "LTr.="<<ltrp<<             // refernece track
605       //
606       "lx0="<<lxyz[0]<<          // reference x
607       "lx1="<<lxyz[1]<<          // reference y
608       "lx2="<<lxyz[2]<<          // refernece z      
609       "lpx0="<<lpxyz[0]<<          // reference x
610       "lpx1="<<lpxyz[1]<<          // reference y
611       "lpx2="<<lpxyz[2]<<          // refernece z            
612       //
613       "msig="<<meanS<<
614       //
615       "mphi="<<meanphi<<         //
616       "rmsphi="<<rmsphi<<        //
617       "gphi1="<<gphi1<<
618       "gphi2="<<gphi2<<
619       //
620       "mphiP="<<meanphiP<<       //
621       "rmsphiP="<<rmsphiP<<      //
622       "gphiP1="<<gphiP1<<
623       "gphiP2="<<gphiP2<<
624       //
625       "meanZ="<<meanZ<<
626       "rmsZ="<<rmsZ<<
627       "gz1="<<gz1<<
628       "gz2="<<gz2<<
629
630       "\n";
631   }
632   delete pcstream;
633 }
634
635
636
637 void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
638   //
639   //
640   //
641   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
642   TFile * f = pcstream->GetFile();
643   f->mkdir("dirphi");
644   f->mkdir("dirphiP");
645   f->mkdir("dirZ");
646   TF1 fp("p1","pol1");
647   //
648   //
649   char    cut[1000];
650   char grname[1000];
651   char grnamefull[1000];
652
653   Double_t mphi[100];
654   Double_t mphiP[100];
655   Double_t smphi[100];
656   Double_t smphiP[100];
657   Double_t mZ[100];
658   Double_t smZ[100];
659   Double_t bz[100];
660   Double_t sbz[100];
661   // fit parameters
662   Double_t pphi[3];
663   Double_t pphiP[3];
664   Double_t pmZ[3];
665   //
666   for (Int_t id=0; id<336; id++){
667     // id =205;
668     sprintf(cut,"isOK&&fId==%d",id);
669     Int_t entries = chain->Draw("bz",cut,"goff");
670     if (entries<3) continue;
671     AliTPCLaserTrack *ltrp = 0;;
672     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
673     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
674     Double_t lxyz[3];
675     Double_t lpxyz[3];
676     ltrp->GetXYZ(lxyz);
677     ltrp->GetPxPyPz(lpxyz);
678
679     chain->Draw("bz",cut,"goff");
680     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
681     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
682     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
683     //
684     chain->Draw("gphi1",cut,"goff");
685     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
686     chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
687     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
688     //
689     chain->Draw("gphiP1",cut,"goff");
690     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
691     chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
692     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
693     //
694     chain->Draw("gz1",cut,"goff");
695     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
696     chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
697     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
698     //
699     //
700     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
701             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
702     // store data  
703     // phi
704     f->cd("dirphi");
705     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
706     grphi->Draw("a*");
707     grphi->Fit(&fp);
708     pphi[0] = fp.GetParameter(0);                          // offset
709     pphi[1] = fp.GetParameter(1);                          // slope
710     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
711     sprintf(grname,"phi_id%d",id);
712     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
713     grphi->GetXaxis()->SetTitle("b_{z} (T)");
714     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
715     grphi->Draw("a*");
716
717     grphi->Write();
718     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
719     // phiP
720     f->cd("dirphiP)");
721     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
722     grphiP->Draw("a*");
723     grphiP->Fit(&fp);
724     pphiP[0] = fp.GetParameter(0);                          // offset
725     pphiP[1] = fp.GetParameter(1);                          // slope
726     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
727     sprintf(grname,"phiP_id%d",id);
728     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
729     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
730     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
731
732     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
733     grphiP->Write();
734     //
735     //Z
736     f->cd("dirZ");
737     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
738     grmZ->Draw("a*");
739     grmZ->Fit(&fp);
740     pmZ[0] = fp.GetParameter(0);                          // offset
741     pmZ[1] = fp.GetParameter(1);                          // slope
742     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
743     sprintf(grname,"mZ_id%d",id);
744     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
745     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
746     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
747
748     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
749     grmZ->Write();
750     
751
752     for (Int_t ientry=0; ientry<entries; ientry++){
753       (*pcstream)<<"Mean"<<
754         "id="<<id<<
755         "LTr.="<<ltrp<<
756         "entries="<<entries<<
757         "bz="<<bz[ientry]<<
758         "lx0="<<lxyz[0]<<          // reference x
759         "lx1="<<lxyz[1]<<          // reference y
760         "lx2="<<lxyz[2]<<          // refernece z      
761         "lpx0="<<lpxyz[0]<<          // reference x
762         "lpx1="<<lpxyz[1]<<          // reference y
763         "lpx2="<<lpxyz[2]<<          // refernece z            
764         //values
765         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
766         "pphi0="<<pphi[0]<<   // offset
767         "pphi1="<<pphi[1]<<   // mean
768         "pphi2="<<pphi[2]<<   // norm chi2
769         //
770         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
771         "pphiP0="<<pphiP[0]<< // offset
772         "pphiP1="<<pphiP[1]<< // mean
773         "pphiP2="<<pphiP[2]<< // norm chi2
774         //
775         "gz1="<<mZ[ientry]<<
776         "pmZ0="<<pmZ[0]<<     // offset
777         "pmZ1="<<pmZ[1]<<     // mean
778         "pmZ2="<<pmZ[2]<<     // norm chi2
779         "\n";
780     }
781   }
782   
783   delete pcstream;
784   
785 }
786
787
788 void AliTPCcalibLaser::Analyze(){
789   //
790   //
791   //
792 }
793
794
795 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
796
797   TIterator* iter = li->MakeIterator();
798   AliTPCcalibLaser* cal = 0;
799
800   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
801     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
802       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
803       return -1;
804     }
805     //
806    //  fHistNTracks->Add(cal->fHistNTracks);
807 //     fClusters->Add(cal->fClusters);
808 //     fModules->Add(cal->fModules);
809 //     fHistPt->Add(cal->fHistPt);
810 //     fPtResolution->Add(cal->fPtResolution);
811 //     fDeDx->Add(cal->fDeDx);
812     
813
814     TH1F *h=0x0;
815     TH1F *hm=0x0;
816
817     for (Int_t id=0; id<336; id++){
818       // merge fDeltaZ histograms
819       hm = (TH1F*)cal->fDeltaZ.At(id); 
820       h  = (TH1F*)fDeltaZ.At(id);
821       if (!h) {
822         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
823         fDeltaZ.AddAt(h,id);
824       }
825       if (hm) h->Add(hm);
826       // merge fDeltaPhi histograms
827       hm = (TH1F*)cal->fDeltaPhi.At(id); 
828       h  = (TH1F*)fDeltaPhi.At(id);
829       if (!h) {
830         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
831         fDeltaPhi.AddAt(h,id);
832       }
833       if (hm) h->Add(hm);
834       // merge fDeltaPhiP histograms
835       hm = (TH1F*)cal->fDeltaPhiP.At(id); 
836       h  = (TH1F*)fDeltaPhiP.At(id);
837       if (!h) {
838         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
839         fDeltaPhiP.AddAt(h,id);
840       }
841       if (hm) h->Add(hm);
842       // merge fSignals histograms
843       hm = (TH1F*)cal->fSignals.At(id); 
844       h  = (TH1F*)fSignals.At(id);
845       if (!h) {
846         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
847         fSignals.AddAt(h,id);
848       }
849       if (hm) h->Add(hm);      
850     }
851   }
852   return 0;
853 }
854
855
856
857 /*
858  gSystem->Load("libSTAT.so")
859  TStatToolkit toolkit;
860  Double_t chi2;
861  TVectorD fitParam;
862  TMatrixD covMatrix;
863  Int_t npoints;
864  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003");
865
866
867 TString fstring="";
868 //
869 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
870 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
871 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
872 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4 
873 //
874 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
875 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
876 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
877 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
878 //   
879 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
880 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
881 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
882 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
883
884
885
886
887  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
888
889  treeT->SetAlias("fit",strq0->Data());
890  
891
892  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
893
894  treeT->SetAlias("fitP",strqP->Data());
895
896
897
898
899
900
901
902  
903  */