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  */