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         profy->SetDirectory(0);
 906         fDeltaYres.AddAt(profy,id);
 907       }
 908       if (!profz){
 909         profz=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
 910         profz->SetDirectory(0);
 911         fDeltaZres.AddAt(profz,id);
 912       }
 913       for (Int_t irow=158;irow>-1;--irow) {
 914         if (vecSec[irow]==-1)continue; //no cluster info
 915         Double_t x    = vecX[irow];
 916         Double_t ycl  = vecClY[irow];
 917         Double_t yfit = vecY1[irow];
 918         Double_t zcl  = vecClZ[irow];
 919         Double_t zfit = vecZ1[irow];
 920         if (profy)
 921           if (profy->GetEntries()<1000000)
 922             profy->Fill(x,yfit-ycl);
 923         if (profz)
 924           if (profz->GetEntries()<1000000)
 925             profz->Fill(x,zfit-zcl);
 926       }
 927       //===============================//
 928       // Fill Fit Parameter Histograms //
 929       //===============================//
 930       TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
 931       TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
 932       TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
 933       TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
 934       TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
 935       TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
 936       TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
 937       TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
 938       //create histograms if the do not already exist
 939       if (!pol2InnerY){
 940           pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
 941                                Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
 942                                500,-.005,.005);
 943           pol2InnerY->SetDirectory(0);
 944           pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
 945                                Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
 946                                500,0.01,.01);
 947           pol2OuterY->SetDirectory(0);
 948
 949           diff1InnerY=new TH1F(Form("diff1inY%03d",id),
 950                                Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
 951                                500,-.5,.5);
 952           diff1InnerY->SetDirectory(0);
 953
 954           diff1OuterY=new TH1F(Form("diff1outY%03d",id),
 955                                Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
 956                                500,-1,1);
 957           diff1OuterY->SetDirectory(0);
 958
 959           pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
 960                                Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
 961                                500,-.002,.002);
 962           pol2InnerZ->SetDirectory(0);
 963
 964           pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
 965                                Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
 966                                500,-.005,.005);
 967           pol2OuterZ->SetDirectory(0);
 968           diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
 969                                Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
 970                                500,-.02,.02);
 971           diff1InnerZ->SetDirectory(0);
 972           diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
 973                                Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
 974                                500,-.03,.03);
 975           diff1OuterZ->SetDirectory(0);
 976           //add
 977           fPol2Par2InY.AddAt(pol2InnerY,id);
 978           fDiffPar1InY.AddAt(diff1InnerY,id);
 979           fPol2Par2OutY.AddAt(pol2OuterY,id);
 980           fDiffPar1OutY.AddAt(diff1OuterY,id);
 981           fPol2Par2InZ.AddAt(pol2InnerZ,id);
 982           fDiffPar1InZ.AddAt(diff1InnerZ,id);
 983           fPol2Par2OutZ.AddAt(pol2OuterZ,id);
 984           fDiffPar1OutZ.AddAt(diff1OuterZ,id);
 985       }
 986       //fill histograms
 987       pol2InnerY ->Fill(vecy2resInner[2]);
 988       pol2OuterY ->Fill(vecy2resOuter[2]);
 989       diff1InnerY->Fill(vecy2resInner[1]-vecy1resInner[1]);
 990       diff1OuterY->Fill(vecy2resOuter[1]-vecy1resOuter[1]);
 991       pol2InnerZ ->Fill(vecz2resInner[2]);
 992       pol2OuterZ ->Fill(vecz2resOuter[2]);
 993       diff1InnerZ->Fill(vecz2resInner[1]-vecz1resInner[1]);
 994       diff1OuterZ->Fill(vecz2resOuter[1]-vecz1resOuter[1]);
 995
 996
 997
 998   }
 999 }
1000
1001
1002 void AliTPCcalibLaser::RefitLaser(Int_t id){
1003   //
1004   // Refit the track store residuals
1005   //
1006
1007   AliTPCseed *track    = (AliTPCseed*)fTracksTPC.At(id);
1008   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1009   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1010
1011   //linear fit model in y and z per sector
1012   static TLinearFitter fy1(2,"hyp1");
1013   static TLinearFitter fz1(2,"hyp1");
1014   //quadratic fit model in y and z per sector
1015   static TLinearFitter fy2(3,"hyp2");
1016   static TLinearFitter fz2(3,"hyp2");
1017   static TVectorD vecy2,vecz2,vecy1,vecz1;
1018
1019   const Int_t kMinClusters=20;
1020   Int_t nclusters[72];
1021   //
1022   for (Int_t i=0;i<72;++i) nclusters[i]=0;
1023
1024   for (Int_t i=0;i<160;++i) {
1025     AliTPCclusterMI *c=track->GetClusterPointer(i);
1026     if (c) nclusters[c->GetDetector()]++;
1027   }
1028
1029   for (Int_t isec=0; isec<72;isec++){
1030     if (nclusters[isec]<kMinClusters) continue;
1031     fy2.ClearPoints();
1032     fz2.ClearPoints();
1033     fy1.ClearPoints();
1034     fz1.ClearPoints();
1035     //
1036     for (Int_t irow=0;irow<160;++irow) {
1037       AliTPCclusterMI *c=track->GetClusterPointer(irow);
1038       //if (c && RejectCluster(c)) continue;
1039       Double_t xd = c->GetX()-133.4; // reference x is beteen iroc and oroc
1040       if (c&&c->GetDetector()==isec) {
1041         Double_t x[2]={xd,xd*xd};
1042         fy2.AddPoint(x,c->GetY());
1043         fz2.AddPoint(x,c->GetZ());
1044         //
1045         fy1.AddPoint(x,c->GetY());
1046         fz1.AddPoint(x,c->GetZ());
1047       }
1048     }
1049     fy2.Eval();
1050     fz2.Eval();
1051     fy1.Eval();
1052     fz1.Eval();
1053     fy1.GetParameters(vecy1);
1054     fy2.GetParameters(vecy2);
1055     fz1.GetParameters(vecz1);
1056     fz2.GetParameters(vecz2);
1057
1058     if (fStreamLevel>0){
1059       TTreeSRedirector *cstream = GetDebugStreamer();
1060       if (cstream){
1061         Float_t dedx = track->GetdEdx();
1062         (*cstream)<<"Tracklet"<<
1063           "LTr.="<<ltrp<<
1064           "Tr.="<<param<<
1065           "sec="<<isec<<
1066           "ncl="<<nclusters[isec]<<
1067           "dedx="<<dedx<<
1068           "dedx="<<dedx<<
1069           "vy1.="<<&vecy1<<
1070           "vy2.="<<&vecy2<<
1071           "vz1.="<<&vecz1<<
1072           "vz2.="<<&vecz2<<
1073           "\n";
1074       }
1075     }
1076   }
1077   //
1078   //
1079   //
1080   //   for (Int_t irow=0;irow<160;++irow) {
1081   //       AliTPCclusterMI *c=track->GetClusterPointer(irow);
1082   //       if (c && RejectCluster(c)) continue;
1083   //       if (c&&c->GetDetector()==isec) {
1084   //    Double_t x[2]={c->GetX(),c->GetX()*c->GetX()};
1085   //    fy2.AddPoint(&x,c->GetY());
1086   //    fz2.AddPoint(&x,c->GetZ());
1087   //    //
1088   //    fy1.AddPoint(&x,c->GetY());
1089   //    fz1.AddPoint(&x,c->GetZ());
1090   //       }
1091   //     }
1092
1093 }
1094
1095
1096 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield,Int_t minEntries){
1097   //
1098   //  Dump information about laser beams
1099   //  isOK variable indicates usability of the beam
1100   //  Beam is not usable if:
1101   //     a.  No entries in range (krmsCut0)
1102   //     b.  Big sperad          (krmscut1)
1103   //     c.  RMSto fit sigma bigger then (kmultiCut)
1104   //     d.  Too big angular spread
1105   //
1106
1107   const Float_t krmsCut0=0.001;
1108   const Float_t krmsCut1=0.16;
1109   const Float_t kmultiCut=2;
1110   const Float_t kcutP0=0.002;
1111   //
1112   AliTPCcalibLaser *laser = this;
1113   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1114   TF1 fg("fg","gaus");
1115   //
1116   //
1117   for (Int_t id=0; id<336; id++){
1118     Bool_t isOK=kTRUE;
1119     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
1120     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1121     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
1122     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
1123     if (!hisphi) continue;;
1124     Double_t entries = hisphi->GetEntries();
1125     if (entries<minEntries) continue;
1126     //
1127     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1128     if (!ltrp) {
1129      AliTPCLaserTrack::LoadTracks();
1130       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1131     }
1132     Float_t meanphi = hisphi->GetMean();
1133     Float_t rmsphi = hisphi->GetRMS();
1134     //
1135     Float_t meanphiP = hisphiP->GetMean();
1136     Float_t rmsphiP = hisphiP->GetRMS();
1137     Float_t meanZ = hisZ->GetMean();
1138     Float_t rmsZ = hisZ->GetRMS();
1139     hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1140     Double_t gphi1 = fg.GetParameter(1);
1141     Double_t gphi2 = fg.GetParameter(2);
1142     hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1143     Double_t gphiP1 = fg.GetParameter(1);
1144     Double_t gphiP2 = fg.GetParameter(2);
1145     hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
1146     Double_t gz1 = fg.GetParameter(1);
1147     Double_t gz2 = fg.GetParameter(2);
1148     //
1149     Float_t meanS=hisS->GetMean();
1150     //
1151     Double_t lxyz[3];
1152     Double_t lpxyz[3];
1153     ltrp->GetXYZ(lxyz);
1154     ltrp->GetPxPyPz(lpxyz);
1155
1156     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1157     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1158     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
1159     if (gphiP2>kcutP0) isOK=kFALSE;
1160     //
1161     (*pcstream)<<"Mean"<<
1162       "isOK="<<isOK<<
1163       "entries="<<entries<<      // number of entries
1164       "bz="<<bfield<<            // bfield
1165       "LTr.="<<ltrp<<             // refernece track
1166       //
1167       "lx0="<<lxyz[0]<<          // reference x
1168       "lx1="<<lxyz[1]<<          // reference y
1169       "lx2="<<lxyz[2]<<          // refernece z
1170       "lpx0="<<lpxyz[0]<<          // reference x
1171       "lpx1="<<lpxyz[1]<<          // reference y
1172       "lpx2="<<lpxyz[2]<<          // refernece z
1173       //
1174       "msig="<<meanS<<
1175       //
1176       "mphi="<<meanphi<<         //
1177       "rmsphi="<<rmsphi<<        //
1178       "gphi1="<<gphi1<<
1179       "gphi2="<<gphi2<<
1180       //
1181       "mphiP="<<meanphiP<<       //
1182       "rmsphiP="<<rmsphiP<<      //
1183       "gphiP1="<<gphiP1<<
1184       "gphiP2="<<gphiP2<<
1185       //
1186       "meanZ="<<meanZ<<
1187       "rmsZ="<<rmsZ<<
1188       "gz1="<<gz1<<
1189       "gz2="<<gz2<<
1190
1191       "\n";
1192   }
1193   delete pcstream;
1194 }
1195
1196
1197
1198 void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1199   //
1200   //
1201   //
1202   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1203   TFile * f = pcstream->GetFile();
1204   f->mkdir("dirphi");
1205   f->mkdir("dirphiP");
1206   f->mkdir("dirZ");
1207   TF1 fp("p1","pol1");
1208   //
1209   //
1210   char    cut[1000];
1211   char grname[1000];
1212   char grnamefull[1000];
1213
1214   Double_t mphi[100];
1215   Double_t mphiP[100];
1216   Double_t smphi[100];
1217   Double_t smphiP[100];
1218   Double_t mZ[100];
1219   Double_t smZ[100];
1220   Double_t bz[100];
1221   Double_t sbz[100];
1222   // fit parameters
1223   Double_t pphi[3];
1224   Double_t pphiP[3];
1225   Double_t pmZ[3];
1226   //
1227   for (Int_t id=0; id<336; id++){
1228     // id =205;
1229     sprintf(cut,"isOK&&fId==%d",id);
1230     Int_t entries = chain->Draw("bz",cut,"goff");
1231     if (entries<3) continue;
1232     AliTPCLaserTrack *ltrp = 0;;
1233     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1234     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1235     Double_t lxyz[3];
1236     Double_t lpxyz[3];
1237     ltrp->GetXYZ(lxyz);
1238     ltrp->GetPxPyPz(lpxyz);
1239
1240     chain->Draw("bz",cut,"goff");
1241     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1242     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1243     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1244     //
1245     chain->Draw("gphi1",cut,"goff");
1246     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1247     chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1248     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1249     //
1250     chain->Draw("gphiP1",cut,"goff");
1251     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1252     chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1253     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1254     //
1255     chain->Draw("gz1",cut,"goff");
1256     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1257     chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1258     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1259     //
1260     //
1261     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1262             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1263     // store data
1264     // phi
1265     f->cd("dirphi");
1266     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1267     grphi->Draw("a*");
1268     grphi->Fit(&fp);
1269     pphi[0] = fp.GetParameter(0);                          // offset
1270     pphi[1] = fp.GetParameter(1);                          // slope
1271     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1272     sprintf(grname,"phi_id%d",id);
1273     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
1274     grphi->GetXaxis()->SetTitle("b_{z} (T)");
1275     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
1276     grphi->SetMaximum(1.2);
1277     grphi->SetMinimum(-1.2);
1278     grphi->Draw("a*");
1279
1280     grphi->Write();
1281     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1282     // phiP
1283     f->cd("dirphiP)");
1284     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1285     grphiP->Draw("a*");
1286     grphiP->Fit(&fp);
1287     pphiP[0] = fp.GetParameter(0);                          // offset
1288     pphiP[1] = fp.GetParameter(1);                          // slope
1289     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1290     sprintf(grname,"phiP_id%d",id);
1291     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
1292     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1293     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
1294     grphiP->SetMaximum(pphiP[0]+0.005);
1295     grphiP->SetMinimum(pphiP[0]-0.005);
1296
1297     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1298     grphiP->Write();
1299     //
1300     //Z
1301     f->cd("dirZ");
1302     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1303     grmZ->Draw("a*");
1304     grmZ->Fit(&fp);
1305     pmZ[0] = fp.GetParameter(0);                          // offset
1306     pmZ[1] = fp.GetParameter(1);                          // slope
1307     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1308     sprintf(grname,"mZ_id%d",id);
1309     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
1310     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1311     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1312
1313     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1314     grmZ->Write();
1315
1316
1317     for (Int_t ientry=0; ientry<entries; ientry++){
1318       (*pcstream)<<"Mean"<<
1319         "id="<<id<<
1320         "LTr.="<<ltrp<<
1321         "entries="<<entries<<
1322         "bz="<<bz[ientry]<<
1323         "lx0="<<lxyz[0]<<          // reference x
1324         "lx1="<<lxyz[1]<<          // reference y
1325         "lx2="<<lxyz[2]<<          // refernece z
1326         "lpx0="<<lpxyz[0]<<          // reference x
1327         "lpx1="<<lpxyz[1]<<          // reference y
1328         "lpx2="<<lpxyz[2]<<          // refernece z
1329         //values
1330         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1331         "pphi0="<<pphi[0]<<   // offset
1332         "pphi1="<<pphi[1]<<   // mean
1333         "pphi2="<<pphi[2]<<   // norm chi2
1334         //
1335         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1336         "pphiP0="<<pphiP[0]<< // offset
1337         "pphiP1="<<pphiP[1]<< // mean
1338         "pphiP2="<<pphiP[2]<< // norm chi2
1339         //
1340         "gz1="<<mZ[ientry]<<
1341         "pmZ0="<<pmZ[0]<<     // offset
1342         "pmZ1="<<pmZ[1]<<     // mean
1343         "pmZ2="<<pmZ[2]<<     // norm chi2
1344         "\n";
1345     }
1346   }
1347
1348   delete pcstream;
1349
1350 }
1351
1352
1353 void AliTPCcalibLaser::Analyze(){
1354   //
1355   //
1356   //
1357 }
1358
1359
1360 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
1361
1362   TIterator* iter = li->MakeIterator();
1363   AliTPCcalibLaser* cal = 0;
1364
1365   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1366     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1367       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1368       return -1;
1369     }
1370     //
1371    //  fHistNTracks->Add(cal->fHistNTracks);
1372 //     fClusters->Add(cal->fClusters);
1373 //     fModules->Add(cal->fModules);
1374 //     fHistPt->Add(cal->fHistPt);
1375 //     fPtResolution->Add(cal->fPtResolution);
1376 //     fDeDx->Add(cal->fDeDx);
1377
1378
1379     TH1F *h=0x0;
1380     TH1F *hm=0x0;
1381     TProfile *hp=0x0;
1382     TProfile *hpm=0x0;
1383
1384     for (Int_t id=0; id<336; id++){
1385       // merge fDeltaZ histograms
1386       hm = (TH1F*)cal->fDeltaZ.At(id);
1387       h  = (TH1F*)fDeltaZ.At(id);
1388       if (!h) {
1389         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1390         h->SetDirectory(0);
1391         fDeltaZ.AddAt(h,id);
1392       }
1393       if (hm) h->Add(hm);
1394       // merge fDeltaPhi histograms
1395       hm = (TH1F*)cal->fDeltaPhi.At(id);
1396       h  = (TH1F*)fDeltaPhi.At(id);
1397       if (!h) {
1398         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
1399         h->SetDirectory(0);
1400         fDeltaPhi.AddAt(h,id);
1401       }
1402       if (hm) h->Add(hm);
1403       // merge fDeltaPhiP histograms
1404       hm = (TH1F*)cal->fDeltaPhiP.At(id);
1405       h  = (TH1F*)fDeltaPhiP.At(id);
1406       if (!h) {
1407         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
1408         h->SetDirectory(0);
1409         fDeltaPhiP.AddAt(h,id);
1410       }
1411       if (hm) h->Add(hm);
1412       // merge fSignals histograms
1413       hm = (TH1F*)cal->fSignals.At(id);
1414       h  = (TH1F*)fSignals.At(id);
1415       if (!h) {
1416         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000);
1417         h->SetDirectory(0);
1418         fSignals.AddAt(h,id);
1419       }
1420       if (hm) h->Add(hm);
1421       //
1422       //
1423       // merge ProfileY histograms
1424       hpm = (TProfile*)cal->fDeltaYres.At(id);
1425       hp  = (TProfile*)fDeltaYres.At(id);
1426       if (!hp) {
1427         hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),115,80,250);
1428         hp->SetDirectory(0);
1429         fDeltaYres.AddAt(hp,id);
1430       }
1431       if (hpm) hp->Add(hpm);
1432       //
1433       hpm = (TProfile*)cal->fDeltaZres.At(id);
1434       hp  = (TProfile*)fDeltaZres.At(id);
1435       if (!hp) {
1436         hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),115,80,250);
1437         hp->SetDirectory(0);
1438         fDeltaZres.AddAt(hp,id);
1439       }
1440       if (hpm) hp->Add(hpm);
1441       //
1442       //
1443       //merge fit param histograms
1444       //local hists
1445       TH1F *pol2InnerY =(TH1F*)fPol2Par2InY.UncheckedAt(id);
1446       TH1F *diff1InnerY=(TH1F*)fDiffPar1InY.UncheckedAt(id);
1447       TH1F *pol2OuterY =(TH1F*)fPol2Par2OutY.UncheckedAt(id);
1448       TH1F *diff1OuterY=(TH1F*)fDiffPar1OutY.UncheckedAt(id);
1449       TH1F *pol2InnerZ =(TH1F*)fPol2Par2InZ.UncheckedAt(id);
1450       TH1F *diff1InnerZ=(TH1F*)fDiffPar1InZ.UncheckedAt(id);
1451       TH1F *pol2OuterZ =(TH1F*)fPol2Par2OutZ.UncheckedAt(id);
1452       TH1F *diff1OuterZ=(TH1F*)fDiffPar1OutZ.UncheckedAt(id);
1453       //hists to merge
1454       TH1F *pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1455       TH1F *diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1456       TH1F *pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1457       TH1F *diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1458       TH1F *pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1459       TH1F *diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1460       TH1F *pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1461       TH1F *diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1462       //create histos if they do not exist
1463       if (!pol2InnerY){
1464           pol2InnerY =new TH1F(Form("pol2par2inY%03d",id),
1465                                Form("2nd derivative from pol2 fit in Y for Laser beam %03d (inner sector)",id),
1466                                500,-.005,.005);
1467           pol2InnerY->SetDirectory(0);
1468
1469           pol2OuterY =new TH1F(Form("pol2par2outY%03d",id),
1470                                Form("2nd derivative from pol2 fit in Y for Laser beam %03d (outer sector)",id),
1471                                500,0.01,.01);
1472           pol2OuterY->SetDirectory(0);
1473           diff1InnerY=new TH1F(Form("diff1inY%03d",id),
1474                                Form("diff of 1st derivative from pol1 and pol2 in Y fit for Laser beam %03d (inner sector)",id),
1475                                500,-.5,.5);
1476           diff1OuterY=new TH1F(Form("diff1outY%03d",id),
1477                                Form("diff of 1st derivative from pol1 and pol2 in Yfit for Laser beam %03d (outer sector)",id),
1478                                500,-1,1);
1479           diff1InnerY->SetDirectory(0);
1480
1481
1482           pol2InnerZ =new TH1F(Form("pol2par2inZ%03d",id),
1483                                Form("2nd derivative from pol2 fit in Z for Laser beam %03d (inner sector)",id),
1484                                500,-.002,.002);
1485           pol2InnerZ->SetDirectory(0);
1486
1487           pol2OuterZ =new TH1F(Form("pol2par2outZ%03d",id),
1488                                Form("2nd derivative from pol2 fit in Z for Laser beam %03d (outer sector)",id),
1489                                500,-.005,.005);
1490           pol2OuterZ->SetDirectory(0);
1491           diff1InnerZ=new TH1F(Form("diff1inZ%03d",id),
1492                                Form("diff of 1st derivative from pol1 and pol2 in Z fit for Laser beam %03d (inner sector)",id),
1493                                500,-.02,.02);
1494           diff1InnerZ->SetDirectory(0);
1495           diff1OuterZ=new TH1F(Form("diff1outZ%03d",id),
1496                                Form("diff of 1st derivative from pol1 and pol2 in Zfit for Laser beam %03d (outer sector)",id),
1497                                500,-.03,.03);
1498           diff1OuterZ->SetDirectory(0);
1499       }
1500       pol2InnerYm =(TH1F*)cal->fPol2Par2InY.UncheckedAt(id);
1501       diff1InnerYm=(TH1F*)cal->fDiffPar1InY.UncheckedAt(id);
1502       pol2OuterYm =(TH1F*)cal->fPol2Par2OutY.UncheckedAt(id);
1503       diff1OuterYm=(TH1F*)cal->fDiffPar1OutY.UncheckedAt(id);
1504       pol2InnerZm =(TH1F*)cal->fPol2Par2InZ.UncheckedAt(id);
1505       diff1InnerZm=(TH1F*)cal->fDiffPar1InZ.UncheckedAt(id);
1506       pol2OuterZm =(TH1F*)cal->fPol2Par2OutZ.UncheckedAt(id);
1507       diff1OuterZm=(TH1F*)cal->fDiffPar1OutZ.UncheckedAt(id);
1508     }
1509   }
1510   return 0;
1511 }
1512
1513
1514
1515
1516 /*
1517  gSystem->Load("libSTAT.so")
1518  TStatToolkit toolkit;
1519  Double_t chi2;
1520  TVectorD fitParam;
1521  TMatrixD covMatrix;
1522  Int_t npoints;
1523  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
1524
1525
1526 TString fstring="";
1527 //
1528 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
1529 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
1530 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
1531 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4
1532 //
1533 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
1534 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
1535 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
1536 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
1537 //
1538 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
1539 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
1540 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
1541 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
1542
1543
1544
1545
1546  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1547
1548  treeT->SetAlias("fit",strq0->Data());
1549
1550
1551  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
1552
1553  treeT->SetAlias("fitP",strqP->Data());
1554
1555
1556  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
1557  treeT->SetAlias("fitD",strqDrift->Data());
1558
1559
1560 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
1561 {
1562 for (Int_t i=0; i<6;i++){
1563 treeT->SetLineColor(i+2);
1564 treeT->SetMarkerSize(1);
1565 treeT->SetMarkerStyle(22+i);
1566 treeT->SetMarkerColor(i+2);
1567
1568 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
1569 }
1570 }
1571  */
1572
1573
1574
1575 /*
1576   TTree * tree = (TTree*)f.Get("FitModels");
1577
1578   TEventList listLFit0("listLFit0","listLFit0");
1579   TEventList listLFit1("listLFit1","listLFit1");
1580   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
1581   tree->SetEventList(&listLFit0);
1582
1583
1584
1585
1586   gSystem->Load("libSTAT.so")
1587   TStatToolkit toolkit;
1588   Double_t chi2;
1589   TVectorD fitParam;
1590   TMatrixD covMatrix;
1591   Int_t npoints;
1592
1593   chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
1594   chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
1595
1596
1597   TString fstring="";
1598   fstring+="cos(dp)++";
1599   fstring+="sin(dp)++";
1600   fstring+="cos(dt)++";
1601   fstring+="sin(dt)++";
1602
1603   TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
1604
1605
1606 */