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