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