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   // 4. The tracks are refitted at different sectors
  30   //    Fit model
  31   //      4.a) line
  32   //      4.b) parabola
  33   //      4.c) parabola with common P2 for inner and outer
  34   //
  35   // To make laser scan the user interaction neccessary
  36   //
  37   .x ~/UliStyle.C
  38   gSystem->Load("libANALYSIS");
  39   gSystem->Load("libTPCcalib");
  40   TFile fcalib("CalibObjects.root");
  41   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  42   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
  43   laser->DumpMeanInfo(-0,0,10)
  44   TFile fmean("laserMean.root")
  45   //
  46   //  laser track clasification;
  47   //
  48   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
  49   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
  50   TCut cutN("cutN","fTPCncls>70");
  51   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
  52   TCut cutA = cutT+cutPt+cutP;
  53   TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
  54
  55   //
  56   //
  57   // Analyze  LASER scan
  58   //
  59
  60   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  61   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  62   AliXRDPROOFtoolkit tool;
  63   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
  64   chain->Lookup();
  65   AliTPCcalibLaser::DumpScanInfo(chain)
  66   TFile fscan("laserScan.root")
  67   TTree * treeT = (TTree*)fscan.Get("Mean")
  68   //
  69   // Analyze laser
  70   //
  71   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  72   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  73   AliXRDPROOFtoolkit tool;
  74   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
  75   chain->Lookup();
  76   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
  77   chainFit->Lookup();
  78
  79 */
  80
  81
  82
  83 #include "TLinearFitter.h"
  84 #include "AliTPCcalibLaser.h"
  85 #include "AliExternalTrackParam.h"
  86 #include "AliESDEvent.h"
  87 #include "AliESDfriend.h"
  88 #include "AliESDtrack.h"
  89 #include "AliTPCTracklet.h"
  90 #include "TH1D.h"
  91 #include "TH1F.h"
  92 #include "TProfile.h"
  93 #include "TVectorD.h"
  94 #include "TTreeStream.h"
  95 #include "TFile.h"
  96 #include "TF1.h"
  97 #include "TGraphErrors.h"
  98 #include "AliTPCclusterMI.h"
  99 #include "AliTPCseed.h"
 100 #include "AliTracker.h"
 101 #include "AliLog.h"
 102 #include "TClonesArray.h"
 103 #include "TPad.h"
 104 #include "TSystem.h"
 105 #include "TCut.h"
 106 #include "TChain.h"
 107 #include "TH2F.h"
 108 #include "TStatToolkit.h"
 109 #include "TROOT.h"
 110
 111
 112 #include "TTreeStream.h"
 113 #include <iostream>
 114 #include <sstream>
 115 #include "AliTPCLaserTrack.h"
 116
 117 using namespace std;
 118
 119 ClassImp(AliTPCcalibLaser)
 120
 121 AliTPCcalibLaser::AliTPCcalibLaser():
 122   AliTPCcalibBase(),
 123   fESD(0),
 124   fESDfriend(),
 125   fTracksMirror(336),
 126   fTracksEsd(336),
 127   fTracksEsdParam(336),
 128   fTracksTPC(336),
 129   fDeltaZ(336),
 130   fDeltaP3(336),
 131   fDeltaP4(336),
 132   fDeltaPhi(336),
 133   fDeltaPhiP(336),
 134   fSignals(336),
 135   //
 136   fHisNclIn(0),      //->Number of clusters inner
 137   fHisNclOut(0),     //->Number of clusters outer
 138   fHisNclIO(0),      //->Number of cluster inner outer
 139   fHisLclIn(0),      //->Level arm inner
 140   fHisLclOut(0),     //->Level arm outer
 141   fHisLclIO(0),      //->Number of cluster inner outer
 142   fHisdEdx(0),     //->dEdx histo
 143   fHisdZfit(0),     //->distance to the mirror after linear fit
 144   //
 145   //
 146   fHisChi2YIn1(0),      //->chi2 y inner - line
 147   fHisChi2YOut1(0),     //->chi2 y inner - line
 148   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 149   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 150   fHisChi2YIO1(0),      //->chi2 y IO    - common
 151   fHisChi2ZIn1(0),      //->chi2 z inner - line
 152   fHisChi2ZOut1(0),     //->chi2 z inner - line
 153   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 154   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 155   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 156   //
 157   //
 158   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 159   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 160   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 161   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 162   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 163   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 164   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 165   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 166   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 167   //
 168   //
 169   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 170   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 171   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 172   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 173   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 174   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 175   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 176   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 177   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 178   //
 179   fDeltaYres(336),
 180   fDeltaZres(336),
 181   fFitAside(new TVectorD(3)),
 182   fFitCside(new TVectorD(3)),
 183   fEdgeXcuts(3),
 184   fEdgeYcuts(3),
 185   fNClCuts(5),
 186   fNcuts(0),
 187   fRun(0),
 188   fEvent(0)
 189 {
 190   //
 191   // Constructor
 192   //
 193   fTracksEsdParam.SetOwner(kTRUE);
 194 }
 195
 196 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
 197   AliTPCcalibBase(),
 198   fESD(0),
 199   fESDfriend(0),
 200   fTracksMirror(336),
 201   fTracksEsd(336),
 202   fTracksEsdParam(336),
 203   fTracksTPC(336),
 204   //
 205   fDeltaZ(336),          // array of histograms of delta z for each track
 206   fDeltaP3(336),          // array of histograms of delta z for each track
 207   fDeltaP4(336),          // array of histograms of P3 for each track
 208   fDeltaPhi(336),          // array of histograms of P4 for each track
 209   fDeltaPhiP(336),          // array of histograms of delta z for each track
 210   fSignals(336),           // array of dedx signals
 211   //
 212   //
 213   fHisNclIn(0),      //->Number of clusters inner
 214   fHisNclOut(0),     //->Number of clusters outer
 215   fHisNclIO(0),      //->Number of cluster inner outer
 216   fHisLclIn(0),      //->Level arm inner
 217   fHisLclOut(0),     //->Level arm outer
 218   fHisLclIO(0),      //->Number of cluster inner outer
 219   fHisdEdx(0), //->dEdx histo
 220   fHisdZfit(0),     //->distance to the mirror after linear fit
 221   //
 222   //
 223   fHisChi2YIn1(0),      //->chi2 y inner - line
 224   fHisChi2YOut1(0),     //->chi2 y inner - line
 225   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 226   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 227   fHisChi2YIO1(0),      //->chi2 y IO    - common
 228   fHisChi2ZIn1(0),      //->chi2 z inner - line
 229   fHisChi2ZOut1(0),     //->chi2 z inner - line
 230   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 231   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 232   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 233   //
 234   //
 235   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 236   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 237   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 238   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 239   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 240   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 241   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 242   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 243   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 244   //
 245   //
 246   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 247   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 248   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 249   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 250   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 251   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 252   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 253   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 254   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 255   //
 256   //
 257   //
 258   fDeltaYres(336),
 259   fDeltaZres(336),
 260   fFitAside(new TVectorD(3)),        // drift fit - A side
 261   fFitCside(new TVectorD(3)),        // drift fit - C- side
 262   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 263   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 264   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 265   fNcuts(0),           // number of cuts
 266   fRun(0),
 267   fEvent(0)
 268 {
 269   SetName(name);
 270   SetTitle(title);
 271   //
 272   // Constructor
 273   //
 274   fTracksEsdParam.SetOwner(kTRUE);
 275 }
 276
 277 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
 278   AliTPCcalibBase(calibLaser),
 279   fESD(0),
 280   fESDfriend(0),
 281   fTracksMirror(336),
 282   fTracksEsd(336),
 283   fTracksEsdParam(336),
 284   fTracksTPC(336),
 285   //
 286   fDeltaZ(calibLaser.fDeltaZ),          // array of histograms of delta z for each track
 287   fDeltaP3(((calibLaser.fDeltaP3))),          // array of histograms of delta z for each track
 288   fDeltaP4(((calibLaser.fDeltaP4))),          // array of histograms of P3 for each track
 289   fDeltaPhi(((calibLaser.fDeltaPhi))),          // array of histograms of P4 for each track
 290   fDeltaPhiP(((calibLaser.fDeltaPhiP))),          // array of histograms of delta z for each track
 291   fSignals(((calibLaser.fSignals))),           // array of dedx signals
 292   //
 293   //
 294   fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))),      //->Number of clusters inner
 295   fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))),     //->Number of clusters outer
 296   fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))),      //->Number of cluster inner outer
 297   fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))),      //->Level arm inner
 298   fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))),     //->Level arm outer
 299   fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))),      //->Number of cluster inner outer
 300   fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))),      //
 301   fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))),     //->distance to the mirror after linear fit
 302   //
 303   //
 304   fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))),      //->chi2 y inner - line
 305   fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))),     //->chi2 y inner - line
 306   fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))),      //->chi2 y inner - parabola
 307   fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))),     //->chi2 y inner - parabola
 308   fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))),      //->chi2 y IO    - common
 309   fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))),      //->chi2 z inner - line
 310   fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))),     //->chi2 z inner - line
 311   fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))),      //->chi2 z inner - parabola
 312   fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))),     //->chi2 z inner - parabola
 313   fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))),      //->chi2 z IO    - common
 314   //
 315   //
 316   fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))),     //-> delta y   P0outer-P0inner - line
 317   fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))),     //-> delta y   P0outer-P0inner - parabola
 318   fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))),     //-> delta y   P0outer-P0inner - common parabola
 319   fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))),     //-> delta ky  P1outer-P1inner - line
 320   fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))),     //-> delta ky  P1outer-P1inner - parabola
 321   fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))),     //-> delta ky  P1outer-P1inner - common parabola
 322   fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))),   //-> Curv  P2inner - parabola
 323   fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))),  //-> Curv  P2outer - parabola
 324   fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))),   //-> Curv  P2outerinner - common parabola
 325   //
 326   //
 327   fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))),     //-> delta z   P0outer-P0inner - line
 328   fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))),     //-> delta z   P0outer-P0inner - parabola
 329   fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))),     //-> delta z   P0outer-P0inner - common parabola
 330   fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))),     //-> delta kz  P1outer-P1inner - line
 331   fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))),     //-> delta kz  P1outer-P1inner - parabola
 332   fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))),     //-> delta kz  P1outer-P1inner - common parabola
 333   fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))),   //-> Curv  P2inner - parabola
 334   fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))),  //-> Curv  P2outer - parabola
 335   fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))),   //-> Curv  P2outerinner - common parabola
 336   //
 337   //
 338   fDeltaYres(((calibLaser.fDeltaYres))),
 339   fDeltaZres(((calibLaser.fDeltaZres))),
 340   fFitAside(new TVectorD(3)),        // drift fit - A side
 341   fFitCside(new TVectorD(3)),        // drift fit - C- side
 342   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 343   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 344   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 345   fNcuts(0),           // number of cuts
 346   fRun(0),             // run number
 347   fEvent(0)            // current eventnumber
 348 {
 349   //
 350   // copy constructor
 351   //
 352 }
 353
 354
 355
 356 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
 357   //
 358   // assgnment operator
 359   //
 360   if (this != &calibLaser) {
 361     new (this) AliTPCcalibLaser(calibLaser);
 362   }
 363   return *this;
 364
 365 }
 366
 367
 368
 369
 370 AliTPCcalibLaser::~AliTPCcalibLaser() {
 371   //
 372   // destructor
 373   //
 374   if ( fHisNclIn){
 375     delete fHisNclIn;      //->Number of clusters inner
 376     delete fHisNclOut;     //->Number of clusters outer
 377     delete fHisNclIO;      //->Number of cluster inner outer
 378     delete fHisLclIn;      //->Level arm inner
 379     delete fHisLclOut;     //->Level arm outer
 380     delete fHisLclIO;      //->Number of cluster inner outer
 381     delete fHisdEdx;
 382     delete fHisdZfit;
 383     //
 384     //
 385     delete fHisChi2YIn1;      //->chi2 y inner - line
 386     delete fHisChi2YOut1;     //->chi2 y inner - line
 387     delete fHisChi2YIn2;      //->chi2 y inner - parabola
 388     delete fHisChi2YOut2;     //->chi2 y inner - parabola
 389     delete fHisChi2YIO1;      //->chi2 y IO    - common
 390     delete fHisChi2ZIn1;      //->chi2 z inner - line
 391     delete fHisChi2ZOut1;     //->chi2 z inner - line
 392     delete fHisChi2ZIn2;      //->chi2 z inner - parabola
 393     delete fHisChi2ZOut2;     //->chi2 z inner - parabola
 394     delete fHisChi2ZIO1;      //->chi2 z IO    - common
 395     //
 396     //
 397     delete fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
 398     delete fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
 399     delete fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
 400     delete fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
 401     delete fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
 402     delete fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
 403     delete fHisPy2vP2In;   //-> Curv  P2inner - parabola
 404     delete fHisPy2vP2Out;  //-> Curv  P2outer - parabola
 405     delete fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
 406     //
 407     //
 408     delete fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
 409     delete fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
 410     delete fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
 411     delete fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
 412     delete fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
 413     delete fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
 414     delete fHisPz2vP2In;   //-> Curv  P2inner - parabola
 415     delete fHisPz2vP2Out;  //-> Curv  P2outer - parabola
 416     delete fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
 417   }
 418 }
 419
 420
 421
 422 void AliTPCcalibLaser::Process(AliESDEvent * event) {
 423   //
 424   //
 425   // Loop over tracks and call  Process function
 426   //
 427   fESD = event;
 428   if (!fESD) {
 429     return;
 430   }
 431   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 432   if (!fESDfriend) {
 433     return;
 434   }
 435   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
 436   fTracksTPC.Clear();
 437   fTracksEsd.Clear();
 438   fTracksEsdParam.Delete();
 439   for (Int_t id=0; id<336;id++) {
 440     fCounter[id]=0;
 441     fClusterCounter[id]=0;
 442     fClusterSatur[id]=0;
 443   }
 444   static Bool_t init=kFALSE;
 445   if (!init){
 446     init = kTRUE;  // way around for PROOF - to be investigated
 447     MakeFitHistos();
 448   }
 449   //
 450   Int_t n=fESD->GetNumberOfTracks();
 451   Int_t run = fESD->GetRunNumber();
 452   fRun = run;
 453   for (Int_t i=0;i<n;++i) {
 454     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
 455     AliESDtrack *track=fESD->GetTrack(i);
 456     TObject *calibObject=0;
 457     AliTPCseed *seed=0;
 458     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
 459       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
 460         break;
 461     if (track&&seed) {
 462       FindMirror(track,seed);
 463     }
 464     //
 465   }
 466   FitDriftV();
 467   //
 468   for (Int_t id=0; id<336; id++){
 469     //
 470     if (!fTracksEsdParam.At(id)) continue;
 471     if (fClusterSatur[id]>0.3)   continue;  // tracks in saturation
 472     DumpLaser(id);
 473     if ( AcceptLaser(id) && fFitZ[id]<0.5){
 474       RefitLaserJW(id);
 475       MakeDistHisto(id);
 476     }
 477   }
 478
 479  fEvent++;
 480 }
 481
 482 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
 483   //
 484   //
 485   //
 486   //  for (Int_t id=0; id<336; id++){
 487     //
 488     //
 489     if (!fTracksEsdParam.At(id)) return;
 490     if (!AcceptLaser(id)) return;
 491     //
 492     //
 493     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
 494     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
 495     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
 496
 497     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
 498     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
 499     TH1F * hisSignal = (TH1F*)fSignals.At(id);
 500
 501     if (!hisdz){
 502       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
 503       hisdz->SetDirectory(0);
 504       fDeltaZ.AddAt(hisdz,id);
 505
 506       hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
 507       hisP3->SetDirectory(0);
 508       fDeltaP3.AddAt(hisP3,id);
 509       //
 510       hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
 511       hisP4->SetDirectory(0);
 512       fDeltaP4.AddAt(hisP4,id);
 513
 514       //
 515       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
 516       hisdphi->SetDirectory(0);
 517       fDeltaPhi.AddAt(hisdphi,id);
 518       //
 519       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
 520       hisdphiP->SetDirectory(0);
 521       fDeltaPhiP.AddAt(hisdphiP,id);
 522       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
 523       hisSignal->SetDirectory(0);
 524       fSignals.AddAt(hisSignal,id);
 525     }
 526
 527     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 528     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 529     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 530     if (!param) return;
 531     if (!ltrp) return;
 532     if (!track) return;
 533     Double_t xyz[3];
 534     Double_t pxyz[3];
 535     Double_t lxyz[3];
 536     Double_t lpxyz[3];
 537     param->GetXYZ(xyz);
 538     param->GetPxPyPz(pxyz);
 539     ltrp->GetXYZ(lxyz);
 540     ltrp->GetPxPyPz(lpxyz);
 541     //
 542     Float_t dz   = param->GetZ()-ltrp->GetZ();
 543     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
 544     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
 545     if (hisdz) hisdz->Fill(dz);
 546     if (hisP3) hisP3->Fill(param->GetParameter()[3]);
 547     if (hisP4) hisP4->Fill(param->GetParameter()[4]);
 548     if (hisdphi) hisdphi->Fill(dphi);
 549     if (hisdphiP) hisdphiP->Fill(dphiP);
 550     if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
 551     // }
 552 }
 553
 554 void AliTPCcalibLaser::FitDriftV(){
 555   //
 556   // Fit drift velocity - linear approximation in the z and global y
 557   //
 558   Float_t kDistCut           = 4;      // cut on distance to the fitted value
 559   const Float_t kZCut        = 200;    // remove the closest laser beam
 560   const Float_t kSaturCut    = 0.05;   // remove saturated lasers - cut on fraction of saturated
 561   const Float_t kMinClusters = 60;  // minimal amount of the clusters
 562   const Float_t kMinSignal   = 16;  // minimal mean height of the signal
 563   const Float_t kChi2Cut     = 0.1; // chi2 cut to accept drift fit
 564   static TLinearFitter fdriftA(3,"hyp2");
 565   static TLinearFitter fdriftC(3,"hyp2");
 566   static TLinearFitter fdriftA1(2,"hyp1");
 567   static TLinearFitter fdriftC1(2,"hyp1");
 568   TVectorD fitA(3),fitC(3);
 569
 570   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
 571
 572
 573   for (Int_t iter=0; iter<2; iter++){
 574     fdriftA.ClearPoints();
 575     fdriftC.ClearPoints();
 576     fdriftA1.ClearPoints();
 577     fdriftC1.ClearPoints();
 578     //
 579     for (Int_t id=0; id<336; id++){
 580       if (!fTracksEsdParam.At(id)) continue;
 581       if (!AcceptLaser(id)) continue;
 582       if (TMath::Abs(fFitZ[id])>kDistCut-iter) continue;
 583       if ( fClusterSatur[id]>kSaturCut)  continue;
 584       if ( fClusterCounter[id]<kMinClusters)  continue;
 585
 586       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 587       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 588       Double_t xyz[3];
 589       Double_t pxyz[3];
 590       Double_t lxyz[3];
 591       Double_t lpxyz[3];
 592       param->GetXYZ(xyz);
 593       param->GetPxPyPz(pxyz);
 594       ltrp->GetXYZ(lxyz);
 595       ltrp->GetPxPyPz(lpxyz);
 596       if (TMath::Abs(lxyz[2])>kZCut) continue;
 597       Double_t xxx[2] = {lxyz[2],lxyz[1]};
 598       if (ltrp->GetSide()==0){
 599         fdriftA1.AddPoint(xxx,xyz[2],1);
 600         fdriftA.AddPoint(xxx,xyz[2],1);
 601       }else{
 602         fdriftC1.AddPoint(xxx,xyz[2],1);
 603         fdriftC.AddPoint(xxx,xyz[2],1);
 604       }
 605     }
 606     Float_t chi2A = 0;
 607     Float_t chi2C = 0;
 608     Int_t npointsA=0;
 609     Int_t npointsC=0;
 610     //
 611     if (fdriftA.GetNpoints()>10){
 612       fdriftA.Eval();
 613       fdriftA.GetParameters(fitA);
 614       npointsA= fdriftA.GetNpoints();
 615       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
 616       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) (*fFitAside) = fitA;
 617     }
 618     if (fdriftC.GetNpoints()>10){
 619       fdriftC.Eval();
 620       fdriftC.GetParameters(fitC);
 621       npointsC= fdriftC.GetNpoints();
 622       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
 623       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) (*fFitCside) = fitC;
 624     }
 625
 626     for (Int_t id=0; id<336; id++){
 627       if (!fTracksEsdParam.At(id)) continue;
 628       //
 629       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 630       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 631       Double_t xyz[3];
 632       Double_t pxyz[3];
 633       Double_t lxyz[3];
 634       Double_t lpxyz[3];
 635       param->GetXYZ(xyz);
 636       param->GetPxPyPz(pxyz);
 637       ltrp->GetXYZ(lxyz);
 638       ltrp->GetPxPyPz(lpxyz);
 639       Float_t fz =0;
 640       if (ltrp->GetSide()==0){
 641         fz = (fitA)[0]+(fitA)[1]*lxyz[2]+(fitA)[2]*lxyz[1];
 642       }else{
 643         fz = (fitC)[0]+(fitC)[1]*lxyz[2]+(fitC)[2]*lxyz[1];
 644       }
 645       fFitZ[id]=xyz[2]-fz;
 646     }
 647     if (fStreamLevel>0){
 648       TTreeSRedirector *cstream = GetDebugStreamer();
 649       Int_t time = fESD->GetTimeStamp();
 650       if (cstream){
 651         (*cstream)<<"driftv"<<
 652           "iter="<<iter<<
 653           "driftA.="<<fFitAside<<
 654           "driftC.="<<fFitCside<<
 655           "chi2A="<<chi2A<<
 656           "chi2C="<<chi2C<<
 657           "nA="<<npointsA<<
 658           "nC="<<npointsC<<
 659           "time="<<time<<
 660           "\n";
 661       }
 662     }
 663   }
 664 }
 665
 666 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
 667   //
 668   // get distance between mirror and laser track
 669   //
 670   //
 671   Double_t xyz[3];
 672   Double_t lxyz[3];
 673   param->GetXYZ(xyz);
 674   ltrp->GetXYZ(lxyz);
 675   //
 676   //
 677   Double_t dist = 0;
 678   //radial distance
 679   dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
 680   //
 681   // z distance
 682   // apply drift correction if already exist
 683   //
 684   Float_t dz = 0;
 685   if (ltrp->GetSide()==0){
 686     if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
 687   }else{
 688     if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
 689   }
 690   if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
 691   dist+=TMath::Abs(dz);
 692   //
 693   // phi dist - divergence on 50 cm
 694   //
 695   dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
 696   return dist;
 697 }
 698
 699
 700 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
 701   //
 702   //
 703   //
 704   /*
 705   TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
 706   TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
 707   TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
 708   TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
 709   TCut cutP4("cutPt","abs(Tr.fP[4])<0.1");
 710
 711   TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
 712   */
 713   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
 714   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 715   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 716   Double_t xyz[3];
 717   Double_t lxyz[3];
 718   param->GetXYZ(xyz);
 719   ltrp->GetXYZ(lxyz);
 720   if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
 721   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE;    // cutZ -P1
 722   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE;  // cut -P2
 723   if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE;   // cut Tl -P3
 724   if (TMath::Abs(param->GetParameter()[4])>0.1) return kFALSE;   // cut Pt  -P4
 725   //
 726   //
 727
 728   return kTRUE;
 729 }
 730
 731 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
 732   //
 733   // Find corresponding mirror
 734   // add the corresponding tracks
 735   //
 736   Float_t kRadius0 = 252;
 737   Float_t kRadius  = 253.4;
 738   if (!track->GetOuterParam()) return -1;
 739   AliExternalTrackParam param(*(track->GetOuterParam()));
 740   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
 741   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
 742   AliTPCLaserTrack ltr;
 743   AliTPCLaserTrack *ltrp=0x0;
 744   //  AliTPCLaserTrack *ltrpjw=0x0;
 745   //
 746   Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param);
 747  // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
 748   //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
 749
 750   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
 751     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
 752   else
 753     ltrp=&ltr;
 754
 755   if (id<0) return id;
 756   fCounter[id]++;
 757   //
 758   //
 759   //
 760   Int_t countercl=0;
 761   Float_t counterSatur=0;
 762   for (Int_t irow=158;irow>-1;--irow) {
 763     AliTPCclusterMI *c=seed->GetClusterPointer(irow);
 764     if (!c) continue;
 765     Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
 766     Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
 767     if (pedgeY<3) continue;
 768     if (pedgeX<3) continue;
 769     countercl++;
 770     if (c->GetMax()>900) counterSatur++;
 771   }
 772   counterSatur/=(countercl+1);
 773   //
 774   if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
 775   //
 776   //
 777   Float_t radius=TMath::Abs(ltrp->GetX());
 778   AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
 779   //
 780   if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
 781   Bool_t accept=kTRUE;
 782   //
 783   // choose closer track
 784   //
 785   AliExternalTrackParam * param0  = (AliExternalTrackParam *)fTracksEsdParam.At(id);
 786   if (param0){
 787     Float_t dist0=GetDistance(param0,ltrp);
 788     Float_t dist1=GetDistance(&param,ltrp);
 789     if (dist0<dist1)    accept=kFALSE;
 790   }
 791
 792   if (accept){
 793     fClusterCounter[id]=countercl;
 794     fTracksEsdParam.AddAt(param.Clone(),id);
 795     fTracksEsd.AddAt(track,id);
 796     fTracksTPC.AddAt(seed,id);
 797   }
 798   return id;
 799 }
 800
 801
 802
 803 void AliTPCcalibLaser::DumpLaser(Int_t id) {
 804   //
 805   //  Dump Laser info to the tree
 806   //
 807   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 808   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 809   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 810   //
 811   // Fast laser ID
 812   //
 813   Double_t xyz[3];
 814   Double_t pxyz[3];
 815   Double_t lxyz[3];
 816   Double_t lpxyz[3];
 817   param->GetXYZ(xyz);
 818   param->GetPxPyPz(pxyz);
 819   ltrp->GetXYZ(lxyz);
 820   ltrp->GetPxPyPz(lpxyz);
 821   Float_t dist3D   = GetDistance(param,ltrp);
 822   Float_t dist0    = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
 823   Float_t distphi  = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
 824
 825
 826   if (fStreamLevel>0){
 827     TTreeSRedirector *cstream = GetDebugStreamer();
 828     Int_t time = fESD->GetTimeStamp();
 829     Bool_t accept = AcceptLaser(id);
 830     if (cstream){
 831       (*cstream)<<"Track"<<
 832         "run="<<fRun<<
 833         "id="<<id<<
 834         "accept="<<accept<<
 835         "driftA.="<<fFitAside<<
 836         "driftC.="<<fFitCside<<
 837         "time="<<time<<
 838         "dist3D="<<dist3D<<
 839         "dist0="<<dist0<<
 840         "distphi="<<distphi<<
 841         //
 842         //
 843         "counter="<<fCounter[id]<<
 844         "clcounter="<<fClusterCounter[id]<<
 845         "satur="<<fClusterSatur[id]<<
 846         "fz="<<fFitZ[id]<<
 847         //
 848         "LTr.="<<ltrp<<
 849         "Esd.="<<track<<
 850         "Tr.="<<param<<
 851         "x0="<<xyz[0]<<
 852         "x1="<<xyz[1]<<
 853         "x2="<<xyz[2]<<
 854         "px0="<<pxyz[0]<<
 855         "px1="<<pxyz[1]<<
 856         "px2="<<pxyz[2]<<
 857         //
 858         "lx0="<<lxyz[0]<<
 859         "lx1="<<lxyz[1]<<
 860         "lx2="<<lxyz[2]<<
 861         "lpx0="<<lpxyz[0]<<
 862         "lpx1="<<lpxyz[1]<<
 863         "lpx2="<<lpxyz[2]<<
 864         "\n";
 865     }
 866   }
 867 }
 868
 869 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
 870   //
 871   // Refit the track with different tracklet models:
 872   // 1. Per ROC using the kalman filter, different edge cuts
 873   // 2. Per ROC linear in y and z
 874   // 3. Per ROC quadratic in y and z
 875   // 4. Per track offset for each sector, linear for each sector, common quadratic
 876   // store x, y, z information for all models and the cluster to calculate the residuals
 877   //
 878
 879   AliTPCseed *track      = (AliTPCseed*)fTracksTPC.At(id);
 880   AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 881   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
 882
 883   AliTPCclusterMI dummyCl;
 884
 885   //two tracklets
 886   Int_t kMaxTracklets=2;
 887   //=============================================//
 888   // Linear Fitters for the Different Approaches //
 889   //=============================================//
 890   //linear fit model in y and z; inner - outer sector
 891   static TLinearFitter fy1I(2,"hyp1");
 892   static TLinearFitter fy1O(2,"hyp1");
 893   static TLinearFitter fz1I(2,"hyp1");
 894   static TLinearFitter fz1O(2,"hyp1");
 895   //quadratic fit model in y and z; inner - sector
 896   static TLinearFitter fy2I(3,"hyp2");
 897   static TLinearFitter fy2O(3,"hyp2");
 898   static TLinearFitter fz2I(3,"hyp2");
 899   static TLinearFitter fz2O(3,"hyp2");
 900   //common quadratic fit for IROC and OROC in y and z
 901   static TLinearFitter fy4(5,"hyp4");
 902   static TLinearFitter fz4(5,"hyp4");
 903
 904
 905   //set standard cuts
 906   if ( fNcuts==0 ){
 907       fNcuts=1;
 908       fEdgeXcuts[0]=4;
 909       fEdgeYcuts[0]=3;
 910       fNClCuts[0]=20;
 911   }
 912   //=============================//
 913   // Loop over all Tracklet Cuts //
 914   //=============================//
 915   for (Int_t icut=0; icut<fNcuts; icut++){
 916     Float_t xinMin = 250, xinMax=90;
 917     Float_t xoutMin=250, xoutMax=90;
 918
 919
 920       AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
 921       //cut parameters
 922       Double_t edgeCutX = fEdgeXcuts[icut];
 923       Double_t edgeCutY = fEdgeYcuts[icut];
 924       Int_t    nclCut   = (Int_t)fNClCuts[icut];
 925       //=========================//
 926       // Parameters to calculate //
 927       //=========================//
 928       //fit parameter inner, outer tracklet and combined fit
 929       TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
 930       TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
 931       //
 932       TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
 933       TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
 934       TVectorD vecy4res(5),vecz4res(5);
 935       // cluster and track positions for each row - used for residuals
 936       TVectorD vecX(159);        // x is the same for all (row center)
 937       TVectorD vecYkalman(159);  // y from kalman fit
 938       TVectorD vecZkalman(159);  // z from kalman fit
 939       TVectorD vecY1(159);       // y from pol1 fit per ROC
 940       TVectorD vecZ1(159);       // z from pol1 fit per ROC
 941       TVectorD vecY2(159);       // y from pol2 fit per ROC
 942       TVectorD vecZ2(159);       // z from pol2 fit per ROC
 943       TVectorD vecY4(159);       // y from sector fit
 944       TVectorD vecZ4(159);       // z from sector fit
 945       TVectorD vecClY(159);      // y cluster position
 946       TVectorD vecClZ(159);      // z cluster position
 947       TVectorD vecSec(159);      // sector for each row
 948       //chi2 of fits
 949       Double_t chi2I1z=0;       // chi2 of pol1 fit in z (inner)
 950       Double_t chi2I1y=0;       // chi2 of pol1 fit in y (inner)
 951       Double_t chi2O1z=0;       // chi2 of pol1 fit in z (outer)
 952       Double_t chi2O1y=0;       // chi2 of pol1 fit in y (outer)
 953       Double_t chi2I2z=0;       // chi2 of pol2 fit in z (inner)
 954       Double_t chi2I2y=0;       // chi2 of pol2 fit in y (inner)
 955       Double_t chi2O2z=0;       // chi2 of pol2 fit in z (outer)
 956       Double_t chi2O2y=0;       // chi2 of pol2 fit in y (outer)
 957       Double_t chi2IOz=0;       // chi2 of hyp4 fit in z (inner+outer)
 958       Double_t chi2IOy=0;       // chi2 of hyp4 fit in y (inner+outer)
 959       //more
 960       Int_t innerSector = -1;    // number of inner sector
 961       Int_t outerSector = -1;    // number of outer sector
 962       Int_t nclI=0;              // number of clusters (inner)
 963       Int_t nclO=0;              // number of clusters (outer)
 964       //=================================================//
 965       // Perform the Kalman Fit using the Tracklet Class //
 966       //=================================================//
 967       AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
 968       TObjArray tracklets=
 969           AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
 970                                           kFALSE,nclCut,kMaxTracklets);
 971       // tracklet pointers
 972       AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
 973       AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
 974       AliTPCTracklet *tr=0x0;
 975       AliTPCTracklet dummy;
 976       //continue if we didn't find a tracklet
 977       if ( !trInner && !trOuter ) continue;
 978       //================================================================================//
 979       // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
 980       //================================================================================//
 981       if ( trInner && trOuter ){
 982           if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
 983           if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
 984               tr = trInner;
 985               trInner=trOuter;
 986               trOuter=tr;
 987               AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
 988           }
 989       } else {
 990           if ( trInner ){
 991               if ( !trInner->GetInner() ) continue;
 992               trOuter=&dummy;
 993               if ( trInner->GetSector()>35 ){
 994                   trOuter=trInner;
 995                   trInner=&dummy;
 996               }
 997           } else { //trOuter
 998               if ( !trOuter->GetInner() ) continue;
 999               trInner=&dummy;
1000               if ( trOuter->GetSector()<36 ){
1001                   trInner=trOuter;
1002                   trOuter=&dummy;
1003               }
1004           }
1005       }
1006       innerSector = trInner->GetSector();
1007       if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1008       outerSector = trOuter->GetSector();
1009       if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1010
1011       // array of clusters
1012       TClonesArray arrCl("AliTPCclusterMI",159);
1013       arrCl.ExpandCreateFast(159);
1014       //=======================================//
1015       // fill fitters with cluster information //
1016       //=======================================//
1017       AliDebug(3,"Filing Cluster Information");
1018       for (Int_t irow=158;irow>-1;--irow) {
1019           AliTPCclusterMI *c=track->GetClusterPointer(irow);
1020           AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1021           cl=dummyCl;
1022           vecX[irow]   = 0;
1023           vecClY[irow] = 0;
1024           vecClZ[irow] = 0;
1025           Float_t meanY=0, sumY=0;
1026           for (Int_t drow=-1;drow<=1;drow++) {
1027             if (irow+drow<0) continue;
1028             if (irow+drow>158) continue;
1029             AliTPCclusterMI *c=track->GetClusterPointer(irow);
1030             if (!c) continue;
1031             Int_t roc = static_cast<Int_t>(c->GetDetector());
1032             if ( roc!=innerSector && roc!=outerSector ) continue;
1033             if (c->GetX()<10) continue;
1034             if (c->GetY()==0) continue;
1035             meanY+=c->GetY();
1036             sumY++;
1037           }
1038           if (sumY>0)  meanY/=sumY;
1039
1040           //
1041           vecSec[irow]=-1;
1042           if (!c) continue;
1043           Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(meanY);
1044           Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1045
1046           //
1047           Int_t roc = static_cast<Int_t>(c->GetDetector());
1048           if ( roc!=innerSector && roc!=outerSector ) continue;
1049           vecSec[irow]=roc;
1050           //store clusters in clones array
1051           cl=*c;
1052           //
1053           if (c->GetMax()<4)   continue;  // noise cluster?
1054           if (TMath::Abs(c->GetY())<0.0001)   continue;  // noise cluster?
1055           //cluster position
1056           vecX[irow]   = c->GetX();
1057           vecClY[irow] = c->GetY();
1058           vecClZ[irow] = c->GetZ();
1059           //
1060           Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1061           Double_t y = vecClY[irow];
1062           Double_t z = vecClZ[irow];
1063           //
1064           Double_t x2[2]={x,x*x};   //linear and parabolic parameters
1065           Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1066           if ( roc == innerSector ) {
1067               x4[0]=1; //offset inner - outer sector
1068               x4[1]=x; //slope parameter inner sector
1069           } else {
1070               x4[2]=x; //slope parameter outer sector
1071           }
1072           x4[3]=x*x;   //common parabolic shape
1073           if (pedgeX < fEdgeXcuts[icut]) continue;
1074           if (pedgeY < fEdgeYcuts[icut]) continue;
1075           if (c->GetMax()>900) continue;  // cluster in saturation
1076           //
1077           if ( roc==innerSector ){
1078               fy1I.AddPoint(x2,y);
1079               fz1I.AddPoint(x2,z);
1080               fy2I.AddPoint(x2,y);
1081               fz2I.AddPoint(x2,z);
1082               ++nclI;
1083               if (x<xinMin) xinMin=x;
1084               if (x>xinMax) xinMax=x;
1085           }
1086           if ( roc==outerSector ){
1087               fy1O.AddPoint(x2,y);
1088               fz1O.AddPoint(x2,z);
1089               fy2O.AddPoint(x2,y);
1090               fz2O.AddPoint(x2,z);
1091               ++nclO;
1092               if (x<xoutMin) xoutMin=x;
1093               if (x>xoutMax) xoutMax=x;
1094           }
1095           fy4.AddPoint(x4,y);
1096           fz4.AddPoint(x4,z);
1097       }
1098       //======================================//
1099       // Evaluate and retrieve fit parameters //
1100       //======================================//
1101       AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1102       //inner sector
1103       if (  (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1104           fy1I.Eval();
1105           fz1I.Eval();
1106           fy2I.Eval();
1107           fz2I.Eval();
1108           fy1I.GetParameters(vecy1resInner);
1109           fz1I.GetParameters(vecz1resInner);
1110           fy2I.GetParameters(vecy2resInner);
1111           fz2I.GetParameters(vecz2resInner);
1112           chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1113           chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1114           chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1115           chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1116       }
1117       //outer sector
1118       if (  (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1119           fy1O.Eval();
1120           fz1O.Eval();
1121           fy2O.Eval();
1122           fz2O.Eval();
1123           fy1O.GetParameters(vecy1resOuter);
1124           fz1O.GetParameters(vecz1resOuter);
1125           fy2O.GetParameters(vecy2resOuter);
1126           fz2O.GetParameters(vecz2resOuter);
1127           chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1128           chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1129           chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1130           chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1131       }
1132       //combined hyp4 fit
1133       if ( innerSector>0 && outerSector>0 ){
1134           if (fy4.GetNpoints()>0) {
1135               fy4.Eval();
1136               fy4.GetParameters(vecy4res);
1137               chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1138           }
1139           if (fz4.GetNpoints()>0) {
1140               fz4.Eval();
1141               fz4.GetParameters(vecz4res);
1142               chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1143           }
1144       }
1145       //clear points
1146       fy4.ClearPoints();  fz4.ClearPoints();
1147       fy1I.ClearPoints(); fy1O.ClearPoints();
1148       fz1I.ClearPoints(); fz1O.ClearPoints();
1149       fy2I.ClearPoints(); fy2O.ClearPoints();
1150       fz2I.ClearPoints(); fz2O.ClearPoints();
1151       //==============================//
1152       // calculate tracklet positions //
1153       //==============================//
1154       AliDebug(4,"Calculate tracklet positions");
1155       for (Int_t irow=158;irow>-1;--irow) {
1156           if ( vecSec[irow]==-1 ) continue;  //no cluster info
1157           if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1158           tr=&dummy;
1159           Double_t x=vecX[irow];
1160           Double_t xref=x-133.4;
1161           //
1162           Double_t yoffInner=0;
1163           Double_t zoffInner=0;
1164           Double_t yslopeInner=0;
1165           Double_t yslopeOuter=0;
1166           Double_t zslopeInner=0;
1167           Double_t zslopeOuter=0;
1168           //positions of hyperplane fits
1169           if ( vecSec[irow] == outerSector ) {
1170               tr=trOuter;
1171               vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1172               vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1173               vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1174               vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1175               yslopeOuter=vecy4res[3];
1176               zslopeOuter=vecz4res[3];
1177           } else {
1178               tr=trInner;
1179               vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1180               vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1181               vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1182               vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1183               yoffInner=vecy4res[1];
1184               zoffInner=vecz4res[1];
1185               yslopeInner=vecy4res[2];
1186               zslopeInner=vecz4res[2];
1187           }
1188           vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1189           vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1190           //positions of kalman fits
1191           Double_t gxyz[3],xyz[3];
1192           AliExternalTrackParam *param = 0x0;
1193           //
1194           param=tr->GetInner();
1195           if (param){
1196               param->GetXYZ(gxyz);
1197               Float_t bz = AliTracker::GetBz(gxyz);
1198               param->GetYAt(x, bz, xyz[1]);
1199               param->GetZAt(x, bz, xyz[2]);
1200               vecYkalman[irow]=xyz[1];
1201               vecZkalman[irow]=xyz[2];
1202           }
1203       }
1204       //=====================================================================//
1205       // write results from the different tracklet fits with debug streamers //
1206       //=====================================================================//
1207       if (fStreamLevel>4){
1208           TTreeSRedirector *cstream = GetDebugStreamer();
1209           if (cstream){
1210             Float_t dedx = track->GetdEdx();
1211             (*cstream)<<"FitModels"<<
1212               "cutNr="      << icut <<
1213               "edgeCutX="   << edgeCutX <<
1214               "edgeCutY="   << edgeCutY <<
1215               "nclCut="     << nclCut <<
1216               "innerSector="<< innerSector <<
1217               "outerSector="<< outerSector <<
1218               "dEdx="       << dedx <<
1219               "LTr.="       << ltrp <<
1220               "Tr.="        << extparam <<
1221               "yPol1In.="   << &vecy1resInner <<
1222               "zPol1In.="   << &vecz1resInner <<
1223               "yPol2In.="   << &vecy2resInner <<
1224               "zPol2In.="   << &vecz2resInner <<
1225               "yPol1Out.="  << &vecy1resOuter <<
1226               "zPol1Out.="  << &vecz1resOuter <<
1227               "yPol2Out.="  << &vecy2resOuter <<
1228               "zPol2Out.="  << &vecz2resOuter <<
1229               "yInOut.="    << &vecy4res <<
1230               "zInOut.="    << &vecz4res <<
1231               "chi2y1In="   << chi2I1y <<
1232               "chi2z1In="   << chi2I1z <<
1233               "chi2y1Out="  << chi2O1y <<
1234               "chi2z1Out="  << chi2O1z <<
1235               "chi2y2In="   << chi2I2y <<
1236               "chi2z2In="   << chi2I2z <<
1237               "chi2y2Out="  << chi2O2y <<
1238               "chi2z2Out="  << chi2O2z <<
1239               "chi2yInOut=" << chi2IOy <<
1240               "chi2zInOut=" << chi2IOz <<
1241               "trletIn.="   << trInner <<
1242               "trletOut.="  << trOuter <<
1243               "nclI="       << nclI <<
1244               "nclO="       << nclO <<
1245               "xinMin="     << xinMin<<
1246               "xinMax="     << xinMax<<
1247               "xoutMin="    << xoutMin<<
1248               "xoutMax="     << xoutMax<<
1249               "\n";
1250           }
1251       }
1252
1253       // wirte residuals information
1254       if (fStreamLevel>5){
1255           TTreeSRedirector *cstream = GetDebugStreamer();
1256           if (cstream){
1257             Float_t dedx = track->GetdEdx();
1258             (*cstream)<<"Residuals"<<
1259               "cutNr="      << icut <<
1260               "edgeCutX="   << edgeCutX <<
1261               "edgeCutY="   << edgeCutY <<
1262               "nclCut="     << nclCut   <<
1263               "LTr.="       << ltrp <<
1264               "Tr.="        << extparam<<
1265               "dEdx="       << dedx <<
1266               "Cl.="        << &arrCl <<
1267               "TrX.="       << &vecX <<
1268               "TrYpol1.="   << &vecY1 <<
1269               "TrZpol1.="   << &vecZ1 <<
1270               "TrYpol2.="   << &vecY2 <<
1271               "TrZpol2.="   << &vecZ2 <<
1272               "TrYInOut.="  << &vecY4 <<
1273               "TrZInOut.="  << &vecZ4 <<
1274               "ClY.="       << &vecClY <<
1275               "ClZ.="       << &vecClZ <<
1276               "sec.="       << &vecSec <<
1277               "nclI="       << nclI <<
1278               "nclO="       << nclO <<
1279               "xinMin="     << xinMin<<
1280               "xinMax="     << xinMax<<
1281               "xoutMin="    << xoutMin<<
1282               "xoutMax="     << xoutMax<<
1283               "yInOut.="    << &vecy4res <<
1284               "zInOut.="    << &vecz4res <<
1285               "chi2y1In="   << chi2I1y <<    //
1286               "chi2z1In="   << chi2I1z <<
1287               "chi2y1Out="  << chi2O1y <<
1288               "chi2z1Out="  << chi2O1z <<
1289               "chi2y2In="   << chi2I2y <<
1290               "chi2z2In="   << chi2I2z <<
1291               "chi2y2Out="  << chi2O2y <<
1292               "chi2z2Out="  << chi2O2z <<
1293               "chi2yInOut=" << chi2IOy <<
1294               "chi2zInOut=" << chi2IOz <<
1295               "\n";
1296
1297           }
1298       }
1299       //==========================//
1300       // Fill Residual Histograms //
1301       //==========================//
1302       TProfile *profy = (TProfile*)fDeltaYres.UncheckedAt(id);
1303       TProfile *profz = (TProfile*)fDeltaZres.UncheckedAt(id);
1304       if (!profy){
1305         profy=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),160,0,160);
1306         profy->SetDirectory(0);
1307         fDeltaYres.AddAt(profy,id);
1308       }
1309       if (!profz){
1310         profz=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),160,0,160);
1311         profz->SetDirectory(0);
1312         fDeltaZres.AddAt(profz,id);
1313       }
1314       for (Int_t irow=158;irow>-1;--irow) {
1315         if (vecSec[irow]==-1)continue; //no cluster info
1316         Double_t x    = vecX[irow];
1317         Double_t ycl  = vecClY[irow];
1318         Double_t yfit = vecY1[irow];
1319         Double_t zcl  = vecClZ[irow];
1320         Double_t zfit = vecZ1[irow];
1321         if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
1322           if (profy)
1323             if (profy->GetEntries()<1000000)
1324               profy->Fill(irow,yfit-ycl);
1325           if (profz)
1326             if (profz->GetEntries()<1000000)
1327               profz->Fill(irow,zfit-zcl);
1328         }
1329       }
1330       //
1331       //
1332       // Fill laser fit histograms
1333       //
1334       if (!fHisNclIn) MakeFitHistos();
1335       Float_t dedx = track->GetdEdx();
1336       if (nclI>20&&nclO>20){
1337         fHisNclIn->Fill(id,nclI);      //->Number of clusters inner
1338         fHisNclOut->Fill(id,nclO);     //->Number of clusters outer
1339         fHisNclIO->Fill(id,nclI+nclO);      //->Number of cluster inner outer
1340         //
1341         fHisLclIn->Fill(id,xinMax-xinMin);        //->Level arm inner
1342         fHisLclOut->Fill(id,xoutMax-xoutMin);     //->Level arm outer
1343         fHisLclIO->Fill(id,xoutMax-xinMin);       //->Number of cluster inner outer
1344         //
1345         fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
1346         fHisdZfit->Fill(id,fFitZ[id]);
1347         //
1348         //
1349         fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y));      //->chi2 y inner - line
1350         fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y));     //->chi2 y inner - line
1351         fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y));      //->chi2 y inner - parabola
1352         fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y));     //->chi2 y inner - parabola
1353         fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy));      //->chi2 y IO    - common
1354
1355
1356         fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z));      //->chi2 z inner - line
1357         fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z));     //->chi2 z inner - line
1358         fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y));      //->chi2 z inner - parabola
1359         fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z));     //->chi2 z inner - parabola
1360         fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz));      //->chi2 z IO    - common
1361         //
1362         //
1363         fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]);     //-> delta y   P0outer-P0inner - line
1364         fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]);     //-> delta y   P0outer-P0inner - parabola
1365         fHisPy3vP0->Fill(id,vecy4res[1]);                           //-> delta y   P0outer-P0inner - common parabola
1366         //
1367         fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]);     //-> delta ky  P1outer-P1inner - line
1368         fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]);     //-> delta ky  P1outer-P1inner - parabola
1369         fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]);               //-> delta ky  P1outer-P1inner - common parabola
1370         //
1371         fHisPy3vP2IO->Fill(id,vecy4res[4]);        //-> Curv  P2outerinner - common parabola
1372         fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]);     //-> delta z   P0outer-P0inner - line
1373         fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]);     //-> delta z   P0outer-P0inner - parabola
1374         fHisPz3vP0->Fill(id,vecz4res[1]);                           //-> delta z   P0outer-P0inner - common parabola
1375         //
1376         fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]);     //-> delta kz  P1outer-P1inner - line
1377         fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]);     //-> delta kz  P1outer-P1inner - parabola
1378         fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]);               //-> delta kz  P1outer-P1inner - common parabola
1379         fHisPz3vP2IO->Fill(id,vecz4res[4]);        //-> Curv  P2outerinner - common parabola
1380       }
1381       if(nclI>20){
1382         fHisPy2vP2In->Fill(id,vecy2resInner[2]);   //-> Curv  P2inner - parabola
1383         fHisPz2vP2In->Fill(id,vecz2resInner[2]);   //-> Curv  P2inner - parabola
1384       }
1385       //
1386       if (nclO>20){
1387         fHisPz2vP2Out->Fill(id,vecz2resOuter[2]);  //-> Curv  P2outer - parabola
1388         fHisPy2vP2Out->Fill(id,vecy2resOuter[2]);  //-> Curv  P2outer - parabola
1389       }
1390
1391   }
1392 }
1393
1394
1395
1396 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run, Int_t minEntries){
1397   //
1398   //  Dump information about laser beams
1399   //  isOK variable indicates usability of the beam
1400   //  Beam is not usable if:
1401   //     a.  No entries in range (krmsCut0)
1402   //     b.  Big sperad          (krmscut1)
1403   //     c.  RMSto fit sigma bigger then (kmultiCut)
1404   //     d.  Too big angular spread
1405   //
1406
1407   const Float_t krmsCut0=0.001;
1408   const Float_t krmsCut1=0.16;
1409   const Float_t kmultiCut=2;
1410   const Float_t kcutP0=0.002;
1411   //
1412   AliTPCcalibLaser *laser = this;
1413   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1414   TF1 fg("fg","gaus");
1415   //
1416   //
1417   for (Int_t id=0; id<336; id++){
1418     Bool_t isOK=kTRUE;
1419     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
1420     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1421     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
1422     TH1F * hisP3    = (TH1F*)laser->fDeltaP3.At(id);
1423     TH1F * hisP4    = (TH1F*)laser->fDeltaP4.At(id);
1424     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
1425     //if (!hisphi) continue;
1426     Double_t entries = hisphi->GetEntries();
1427     //if (entries<minEntries) continue;
1428     //
1429     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1430     if (!ltrp) {
1431      AliTPCLaserTrack::LoadTracks();
1432       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1433     }
1434     pcstream->GetFile()->cd();
1435     if (hisphi)  hisphi->Write();
1436     if (hisphiP) hisphiP->Write();
1437     if (hisZ)    hisZ->Write();
1438     if (hisP3)    hisP3->Write();
1439     if (hisP4)    hisP4->Write();
1440
1441     Float_t meanphi = hisphi->GetMean();
1442     Float_t rmsphi = hisphi->GetRMS();
1443     //
1444     Float_t meanphiP = hisphiP->GetMean();
1445     Float_t rmsphiP = hisphiP->GetRMS();
1446     Float_t meanZ = hisZ->GetMean();
1447     Float_t rmsZ = hisZ->GetRMS();
1448     hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1449     Double_t gphi1 = fg.GetParameter(1);
1450     Double_t gphi2 = fg.GetParameter(2);
1451     hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1452     Double_t gphiP1 = fg.GetParameter(1);
1453     Double_t gphiP2 = fg.GetParameter(2);
1454     //
1455     hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS());
1456     Double_t gz1 = fg.GetParameter(1);
1457     Double_t gz2 = fg.GetParameter(2);
1458     //
1459     hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
1460     Double_t gp31 = fg.GetParameter(1);
1461     Double_t gp32 = fg.GetParameter(2);
1462     //
1463     hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
1464     Double_t gp41 = fg.GetParameter(1);
1465     Double_t gp42 = fg.GetParameter(2);
1466     //
1467     Float_t meanS=hisS->GetMean();
1468     //
1469     Double_t lxyz[3];
1470     Double_t lpxyz[3];
1471     ltrp->GetXYZ(lxyz);
1472     ltrp->GetPxPyPz(lpxyz);
1473
1474     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1475     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1476     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
1477     if (gphiP2>kcutP0) isOK=kFALSE;
1478     //
1479     //
1480     //
1481     TH1 *his =0;
1482     //
1483     his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
1484     Float_t mnclIn =  his->GetMean();
1485     delete his;
1486     his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
1487     Float_t mnclOut =  his->GetMean();
1488     delete his;
1489     his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
1490     Float_t mnclIO =  his->GetMean();
1491     delete his;
1492     his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
1493     Float_t mLclIn =  his->GetMean();
1494     delete his;
1495     his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
1496     Float_t mLclOut =  his->GetMean();
1497     delete his;
1498     his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
1499     Float_t mLclIO =  his->GetMean();
1500     delete his;
1501     //
1502     his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
1503     Float_t mdEdx =  his->GetMean();
1504     delete his;
1505     //
1506     his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
1507     Float_t mdZfit =  his->GetMean();
1508     delete his;
1509     //
1510     //
1511     //
1512     his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - line
1513     Float_t  mChi2YIn1=  his->GetMean();
1514     delete his;
1515     his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - line
1516     Float_t mChi2YOut1 =  his->GetMean();
1517     delete his;
1518     his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - parabola
1519     Float_t mChi2YIn2 =  his->GetMean();
1520     delete his;
1521     his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - parabola
1522     Float_t mChi2YOut2 =  his->GetMean();
1523     delete his;
1524     his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 y IO    - common
1525     Float_t mChi2YIO1 =  his->GetMean();
1526     delete his;
1527     his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - line
1528     Float_t mChi2ZIn1 =  his->GetMean();
1529     delete his;
1530     his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - line
1531     Float_t mChi2ZOut1 =  his->GetMean();
1532     delete his;
1533     his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - parabola
1534     Float_t mChi2ZIn2 =  his->GetMean();
1535     delete his;
1536     his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - parabola
1537     Float_t mChi2ZOut2 =  his->GetMean();
1538     delete his;
1539     his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 z IO    - common
1540     Float_t mChi2ZIO1 =  his->GetMean();
1541     delete his;
1542     //
1543     // fit res. histos
1544     //
1545     his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - line
1546     Float_t  ePy1vP0  = his->GetEntries();
1547     Float_t  mPy1vP0  = his->GetMean();
1548     Float_t  rPy1vP0  = his->GetRMS();
1549     delete his;
1550
1551     his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - parabola
1552     Float_t  ePy2vP0  = his->GetEntries();
1553     Float_t  mPy2vP0  = his->GetMean();
1554     Float_t  rPy2vP0  = his->GetRMS();
1555     delete his;
1556
1557     his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - common parabola
1558     Float_t  ePy3vP0  = his->GetEntries();
1559     Float_t  mPy3vP0  = his->GetMean();
1560     Float_t  rPy3vP0  = his->GetRMS();
1561     delete his;
1562
1563     his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - line
1564     Float_t  ePy1vP1  = his->GetEntries();
1565     Float_t  mPy1vP1  = his->GetMean();
1566     Float_t  rPy1vP1  = his->GetRMS();
1567     delete his;
1568
1569     his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - parabola
1570     Float_t  ePy2vP1  = his->GetEntries();
1571     Float_t  mPy2vP1  = his->GetMean();
1572     Float_t  rPy2vP1  = his->GetRMS();
1573     delete his;
1574
1575     his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - common parabola
1576     Float_t  ePy3vP1  = his->GetEntries();
1577     Float_t  mPy3vP1  = his->GetMean();
1578     Float_t  rPy3vP1  = his->GetRMS();
1579     delete his;
1580
1581     his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
1582     Float_t  ePy2vP2In  = his->GetEntries();
1583     Float_t  mPy2vP2In  = his->GetMean();
1584     Float_t  rPy2vP2In  = his->GetRMS();
1585     delete his;
1586
1587     his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
1588     Float_t  ePy2vP2Out  = his->GetEntries();
1589     Float_t  mPy2vP2Out  = his->GetMean();
1590     Float_t  rPy2vP2Out  = his->GetRMS();
1591     delete his;
1592
1593     his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
1594     Float_t  ePy3vP2IO  = his->GetEntries();
1595     Float_t  mPy3vP2IO  = his->GetMean();
1596     Float_t  rPy3vP2IO  = his->GetRMS();
1597     delete his;
1598
1599     //
1600     //
1601     his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - line
1602     Float_t  ePz1vP0  = his->GetEntries();
1603     Float_t  mPz1vP0  = his->GetMean();
1604     Float_t  rPz1vP0  = his->GetRMS();
1605     delete his;
1606
1607     his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - parabola
1608     Float_t  ePz2vP0  = his->GetEntries();
1609     Float_t  mPz2vP0  = his->GetMean();
1610     Float_t  rPz2vP0  = his->GetRMS();
1611     delete his;
1612
1613     his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - common parabola
1614     Float_t  ePz3vP0  = his->GetEntries();
1615     Float_t  mPz3vP0  = his->GetMean();
1616     Float_t  rPz3vP0  = his->GetRMS();
1617     delete his;
1618
1619     his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - line
1620     Float_t  ePz1vP1  = his->GetEntries();
1621     Float_t  mPz1vP1  = his->GetMean();
1622     Float_t  rPz1vP1  = his->GetRMS();
1623     delete his;
1624
1625     his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - parabola
1626     Float_t  ePz2vP1  = his->GetEntries();
1627     Float_t  mPz2vP1  = his->GetMean();
1628     Float_t  rPz2vP1  = his->GetRMS();
1629     delete his;
1630
1631     his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - common parabola
1632     Float_t  ePz3vP1  = his->GetEntries();
1633     Float_t  mPz3vP1  = his->GetMean();
1634     Float_t  rPz3vP1  = his->GetRMS();
1635     delete his;
1636
1637     his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
1638     Float_t  ePz2vP2In  = his->GetEntries();
1639     Float_t  mPz2vP2In  = his->GetMean();
1640     Float_t  rPz2vP2In  = his->GetRMS();
1641     delete his;
1642
1643     his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
1644     Float_t  ePz2vP2Out  = his->GetEntries();
1645     Float_t  mPz2vP2Out  = his->GetMean();
1646     Float_t  rPz2vP2Out  = his->GetRMS();
1647     delete his;
1648
1649     his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
1650     Float_t  ePz3vP2IO  = his->GetEntries();
1651     Float_t  mPz3vP2IO  = his->GetMean();
1652     Float_t  rPz3vP2IO  = his->GetRMS();
1653     delete his;
1654
1655     //
1656     if (run<=0) run=fRun;
1657     (*pcstream)<<"Mean"<<
1658       "run="<<run<<               //
1659       "isOK="<<isOK<<             //
1660       "id="<<id<<                 // track id
1661       "entries="<<entries<<       // number of entries
1662       "bz="<<bfield<<             // bfield
1663       "LTr.="<<ltrp<<             // refernece track
1664       //
1665       "mnclIn="<<mnclIn<<         // mean number of clusters in inner
1666       "mnclOut="<<mnclOut<<       // mean number of clusters in outer
1667       "mnclIO="<<mnclIO<<         // mean number of clusters in inner+outer
1668       "mLclIn="<<mLclIn<<         // mean number of clusters in inner
1669       "mLclOut="<<mLclOut<<       // mean number of clusters in outer
1670       "mLclIO="<<mLclIO<<         // mean number of clusters in inner+outer
1671       "mdEdx="<<mdEdx<<           // mean dEdx
1672       "mdZfit="<<mdZfit<<           // mean z fit
1673       //
1674       //
1675       "mChi2YIn1="<<mChi2YIn1<<       //->chi2 y inner - line
1676       "mChi2YOut1="<<mChi2YOut1<<     //->chi2 y inner - line
1677       "mChi2YIn2="<<mChi2YIn2<<       //->chi2 y inner - parabola
1678       "mChi2YOut2="<<mChi2YOut2<<     //->chi2 y inner - parabola
1679       "mChi2YIO1="<<mChi2YIO1<<       //->chi2 y IO    - common
1680       "mChi2ZIn1="<<mChi2ZIn1<<       //->chi2 z inner - line
1681       "mChi2ZOut1="<<mChi2ZOut1<<     //->chi2 z inner - line
1682       "mChi2ZIn2="<<mChi2ZIn2<<       //->chi2 z inner - parabola
1683       "mChi2ZOut2="<<mChi2ZOut2<<     //->chi2 z inner - parabola
1684       "mChi2ZIO1="<<mChi2ZIO1<<       //->chi2 z IO    - common
1685       //
1686       //
1687       "ePy1vP0="<<ePy1vP0<<
1688       "mPy1vP0="<<mPy1vP0<<
1689       "rPy1vP0="<<rPy1vP0<<
1690       "ePy2vP0="<<ePy2vP0<<
1691       "mPy2vP0="<<mPy2vP0<<
1692       "rPy2vP0="<<rPy2vP0<<
1693       "ePy3vP0="<<ePy3vP0<<
1694       "mPy3vP0="<<mPy3vP0<<
1695       "rPy3vP0="<<rPy3vP0<<
1696       "ePy1vP1="<<ePy1vP1<<
1697       "mPy1vP1="<<mPy1vP1<<
1698       "rPy1vP1="<<rPy1vP1<<
1699       "ePy2vP1="<<ePy2vP1<<
1700       "mPy2vP1="<<mPy2vP1<<
1701       "rPy2vP1="<<rPy2vP1<<
1702       "ePy3vP1="<<ePy3vP1<<
1703       "mPy3vP1="<<mPy3vP1<<
1704       "rPy3vP1="<<rPy3vP1<<
1705       "ePy2vP2In="<<ePy2vP2In<<
1706       "mPy2vP2In="<<mPy2vP2In<<
1707       "rPy2vP2In="<<rPy2vP2In<<
1708       "ePy2vP2Out="<<ePy2vP2Out<<
1709       "mPy2vP2Out="<<mPy2vP2Out<<
1710       "rPy2vP2Out="<<rPy2vP2Out<<
1711       "ePy3vP2IO="<<ePy3vP2IO<<
1712       "mPy3vP2IO="<<mPy3vP2IO<<
1713       "rPy3vP2IO="<<rPy3vP2IO<<
1714       //
1715       //
1716       "ePz1vP0="<<ePz1vP0<<
1717       "mPz1vP0="<<mPz1vP0<<
1718       "rPz1vP0="<<rPz1vP0<<
1719       "ePz2vP0="<<ePz2vP0<<
1720       "mPz2vP0="<<mPz2vP0<<
1721       "rPz2vP0="<<rPz2vP0<<
1722       "ePz3vP0="<<ePz3vP0<<
1723       "mPz3vP0="<<mPz3vP0<<
1724       "rPz3vP0="<<rPz3vP0<<
1725       "ePz1vP1="<<ePz1vP1<<
1726       "mPz1vP1="<<mPz1vP1<<
1727       "rPz1vP1="<<rPz1vP1<<
1728       "ePz2vP1="<<ePz2vP1<<
1729       "mPz2vP1="<<mPz2vP1<<
1730       "rPz2vP1="<<rPz2vP1<<
1731       "ePz3vP1="<<ePz3vP1<<
1732       "mPz3vP1="<<mPz3vP1<<
1733       "rPz3vP1="<<rPz3vP1<<
1734       "ePz2vP2In="<<ePz2vP2In<<
1735       "mPz2vP2In="<<mPz2vP2In<<
1736       "rPz2vP2In="<<rPz2vP2In<<
1737       "ePz2vP2Out="<<ePz2vP2Out<<
1738       "mPz2vP2Out="<<mPz2vP2Out<<
1739       "rPz2vP2Out="<<rPz2vP2Out<<
1740       "ePz3vP2IO="<<ePz3vP2IO<<
1741       "mPz3vP2IO="<<mPz3vP2IO<<
1742       "rPz3vP2IO="<<rPz3vP2IO<<
1743       //
1744       //
1745       //
1746       "lx0="<<lxyz[0]<<            // reference x
1747       "lx1="<<lxyz[1]<<            // reference y
1748       "lx2="<<lxyz[2]<<            // refernece z
1749       "lpx0="<<lpxyz[0]<<           // reference x
1750       "lpx1="<<lpxyz[1]<<          // reference y
1751       "lpx2="<<lpxyz[2]<<          // refernece z
1752       //
1753       "msig="<<meanS<<
1754       //
1755       "mphi="<<meanphi<<         //
1756       "rmsphi="<<rmsphi<<        //
1757       "gphi1="<<gphi1<<
1758       "gphi2="<<gphi2<<
1759       //
1760       "mphiP="<<meanphiP<<       //
1761       "rmsphiP="<<rmsphiP<<      //
1762       "gphiP1="<<gphiP1<<
1763       "gphiP2="<<gphiP2<<
1764       //
1765       "meanZ="<<meanZ<<
1766       "rmsZ="<<rmsZ<<
1767       "gz1="<<gz1<<
1768       "gz2="<<gz2<<
1769       //
1770       "gp31="<<gp31<<            //gaus mean - tgl
1771       "gp32="<<gp32<<            //gaus rms  -tgl
1772       "gp41="<<gp41<<            //gaus mean - P4
1773       "gp42="<<gp42<<            //gaus rms  - P4
1774
1775       "\n";
1776   }
1777   delete pcstream;
1778   /*
1779     Browse the content
1780     TFile fmean("laserMean.root");
1781
1782
1783    */
1784
1785
1786 }
1787
1788
1789
1790 void AliTPCcalibLaser::DumpScanInfo(TTree * chain){
1791   //
1792   //
1793   //
1794   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
1795   TFile * f = pcstream->GetFile();
1796   f->mkdir("dirphi");
1797   f->mkdir("dirphiP");
1798   f->mkdir("dirZ");
1799   TF1 fp("p1","pol1");
1800   //
1801   //
1802   char    cut[1000];
1803   char grname[1000];
1804   char grnamefull[1000];
1805
1806   Double_t mphi[100];
1807   Double_t mphiP[100];
1808   Double_t smphi[100];
1809   Double_t smphiP[100];
1810   Double_t mZ[100];
1811   Double_t smZ[100];
1812   Double_t bz[100];
1813   Double_t sbz[100];
1814   // fit parameters
1815   Double_t pphi[3];
1816   Double_t pphiP[3];
1817   Double_t pmZ[3];
1818   //
1819   for (Int_t id=0; id<336; id++){
1820     // id =205;
1821     sprintf(cut,"isOK&&fId==%d",id);
1822     Int_t entries = chain->Draw("bz",cut,"goff");
1823     if (entries<3) continue;
1824     AliTPCLaserTrack *ltrp = 0;
1825     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
1826     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1827     Double_t lxyz[3];
1828     Double_t lpxyz[3];
1829     ltrp->GetXYZ(lxyz);
1830     ltrp->GetPxPyPz(lpxyz);
1831
1832     chain->Draw("bz",cut,"goff");
1833     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
1834     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
1835     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
1836     //
1837     chain->Draw("gphi1",cut,"goff");
1838     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
1839     chain->Draw("0.05*abs(mphi)+gphi2",cut,"goff");
1840     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
1841     //
1842     chain->Draw("gphiP1",cut,"goff");
1843     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
1844     chain->Draw("0.05*abs(mphiP)+gphiP2",cut,"goff");
1845     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
1846     //
1847     chain->Draw("gz1",cut,"goff");
1848     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
1849     chain->Draw("0.01*abs(meanZ)+gz2",cut,"goff");
1850     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
1851     //
1852     //
1853     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
1854             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
1855     // store data
1856     // phi
1857     f->cd("dirphi");
1858     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
1859     grphi->Draw("a*");
1860     grphi->Fit(&fp);
1861     pphi[0] = fp.GetParameter(0);                          // offset
1862     pphi[1] = fp.GetParameter(1);                          // slope
1863     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1864     sprintf(grname,"phi_id%d",id);
1865     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
1866     grphi->GetXaxis()->SetTitle("b_{z} (T)");
1867     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
1868     grphi->SetMaximum(1.2);
1869     grphi->SetMinimum(-1.2);
1870     grphi->Draw("a*");
1871
1872     grphi->Write();
1873     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
1874     // phiP
1875     f->cd("dirphiP)");
1876     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
1877     grphiP->Draw("a*");
1878     grphiP->Fit(&fp);
1879     pphiP[0] = fp.GetParameter(0);                          // offset
1880     pphiP[1] = fp.GetParameter(1);                          // slope
1881     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1882     sprintf(grname,"phiP_id%d",id);
1883     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
1884     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
1885     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
1886     grphiP->SetMaximum(pphiP[0]+0.005);
1887     grphiP->SetMinimum(pphiP[0]-0.005);
1888
1889     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
1890     grphiP->Write();
1891     //
1892     //Z
1893     f->cd("dirZ");
1894     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
1895     grmZ->Draw("a*");
1896     grmZ->Fit(&fp);
1897     pmZ[0] = fp.GetParameter(0);                          // offset
1898     pmZ[1] = fp.GetParameter(1);                          // slope
1899     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
1900     sprintf(grname,"mZ_id%d",id);
1901     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
1902     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
1903     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
1904
1905     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
1906     grmZ->Write();
1907
1908
1909     for (Int_t ientry=0; ientry<entries; ientry++){
1910       (*pcstream)<<"Mean"<<
1911         "id="<<id<<
1912         "LTr.="<<ltrp<<
1913         "entries="<<entries<<
1914         "bz="<<bz[ientry]<<
1915         "lx0="<<lxyz[0]<<          // reference x
1916         "lx1="<<lxyz[1]<<          // reference y
1917         "lx2="<<lxyz[2]<<          // refernece z
1918         "lpx0="<<lpxyz[0]<<          // reference x
1919         "lpx1="<<lpxyz[1]<<          // reference y
1920         "lpx2="<<lpxyz[2]<<          // refernece z
1921         //values
1922         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
1923         "pphi0="<<pphi[0]<<   // offset
1924         "pphi1="<<pphi[1]<<   // mean
1925         "pphi2="<<pphi[2]<<   // norm chi2
1926         //
1927         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
1928         "pphiP0="<<pphiP[0]<< // offset
1929         "pphiP1="<<pphiP[1]<< // mean
1930         "pphiP2="<<pphiP[2]<< // norm chi2
1931         //
1932         "gz1="<<mZ[ientry]<<
1933         "pmZ0="<<pmZ[0]<<     // offset
1934         "pmZ1="<<pmZ[1]<<     // mean
1935         "pmZ2="<<pmZ[2]<<     // norm chi2
1936         "\n";
1937     }
1938   }
1939
1940   delete pcstream;
1941
1942 }
1943
1944
1945 void AliTPCcalibLaser::Analyze(){
1946   //
1947   //
1948   //
1949 }
1950
1951
1952 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
1953
1954   TIterator* iter = li->MakeIterator();
1955   AliTPCcalibLaser* cal = 0;
1956
1957   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
1958     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
1959       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
1960       return -1;
1961     }
1962     //
1963     MergeFitHistos(cal);
1964     TH1F *h=0x0;
1965     TH1F *hm=0x0;
1966     TProfile *hp=0x0;
1967     TProfile *hpm=0x0;
1968
1969     for (Int_t id=0; id<336; id++){
1970       // merge fDeltaZ histograms
1971       hm = (TH1F*)cal->fDeltaZ.At(id);
1972       h  = (TH1F*)fDeltaZ.At(id);
1973       if (!h) {
1974         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
1975         h->SetDirectory(0);
1976         fDeltaZ.AddAt(h,id);
1977       }
1978       if (hm) h->Add(hm);
1979       // merge fP3 histograms
1980       hm = (TH1F*)cal->fDeltaP3.At(id);
1981       h  = (TH1F*)fDeltaP3.At(id);
1982       if (!h) {
1983         h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
1984         h->SetDirectory(0);
1985         fDeltaP3.AddAt(h,id);
1986       }
1987       if (hm) h->Add(hm);
1988       // merge fP4 histograms
1989       hm = (TH1F*)cal->fDeltaP4.At(id);
1990       h  = (TH1F*)fDeltaP4.At(id);
1991       if (!h) {
1992         h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
1993         h->SetDirectory(0);
1994         fDeltaP4.AddAt(h,id);
1995       }
1996       if (hm) h->Add(hm);
1997
1998       //
1999       // merge fDeltaPhi histograms
2000       hm = (TH1F*)cal->fDeltaPhi.At(id);
2001       h  = (TH1F*)fDeltaPhi.At(id);
2002       if (!h) {
2003         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
2004         h->SetDirectory(0);
2005         fDeltaPhi.AddAt(h,id);
2006       }
2007       if (hm) h->Add(hm);
2008       // merge fDeltaPhiP histograms
2009       hm = (TH1F*)cal->fDeltaPhiP.At(id);
2010       h  = (TH1F*)fDeltaPhiP.At(id);
2011       if (!h) {
2012         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
2013         h->SetDirectory(0);
2014         fDeltaPhiP.AddAt(h,id);
2015       }
2016       if (hm) h->Add(hm);
2017       // merge fSignals histograms
2018       hm = (TH1F*)cal->fSignals.At(id);
2019       h  = (TH1F*)fSignals.At(id);
2020       if (!h) {
2021         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
2022         h->SetDirectory(0);
2023         fSignals.AddAt(h,id);
2024       }
2025       if (hm) h->Add(hm);
2026       //
2027       //
2028       // merge ProfileY histograms
2029       hpm = (TProfile*)cal->fDeltaYres.At(id);
2030       hp  = (TProfile*)fDeltaYres.At(id);
2031       if (!hp) {
2032         hp=new TProfile(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d",id),160,0,160);
2033         hp->SetDirectory(0);
2034         fDeltaYres.AddAt(hp,id);
2035       }
2036       if (hpm) hp->Add(hpm);
2037       //
2038       hpm = (TProfile*)cal->fDeltaZres.At(id);
2039       hp  = (TProfile*)fDeltaZres.At(id);
2040       if (!hp) {
2041         hp=new TProfile(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d",id),160,0,160);
2042         hp->SetDirectory(0);
2043         fDeltaZres.AddAt(hp,id);
2044       }
2045       if (hpm) hp->Add(hpm);
2046       //
2047       //
2048     }
2049   }
2050   return 0;
2051 }
2052
2053 void   AliTPCcalibLaser::MakeFitHistos(){
2054   //
2055   // Make a fit histograms
2056   //
2057   // Number of clusters
2058   //
2059   //TH2F           *fHisNclIn;      //->Number of clusters inner
2060   //TH2F           *fHisNclOut;     //->Number of clusters outer
2061   //TH2F           *fHisNclIO;      //->Number of cluster inner outer
2062   // TH2F           *fHisdEdx;       //->dEdx histo
2063   fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
2064   fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
2065   fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
2066   //
2067   fHisLclIn  = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
2068   fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
2069   fHisLclIO  = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
2070   //
2071   fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
2072   fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-1.,1.);
2073
2074   //
2075   //  Chi2
2076   //
2077   //   TH2F           *fHisChi2YIn1;      //->chi2 y inner - line
2078   //   TH2F           *fHisChi2YOut1;     //->chi2 y inner - line
2079   //   TH2F           *fHisChi2YIn2;      //->chi2 y inner - parabola
2080   //   TH2F           *fHisChi2YOut2;     //->chi2 y inner - parabola
2081   //   TH2F           *fHisChi2YIO1;      //->chi2 y IO    - common
2082   fHisChi2YIn1   = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
2083   fHisChi2YOut1  = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
2084   fHisChi2YIn2   = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
2085   fHisChi2YOut2  = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
2086   fHisChi2YIO1   = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
2087   //   TH2F           *fHisChi2ZIn1;      //->chi2 z inner - line
2088   //   TH2F           *fHisChi2ZOut1;     //->chi2 z inner - line
2089   //   TH2F           *fHisChi2ZIn2;      //->chi2 z inner - parabola
2090   //   TH2F           *fHisChi2ZOut2;     //->chi2 z inner - parabola
2091   //   TH2F           *fHisChi2ZIO1;      //->chi2 z IO    - common
2092   fHisChi2ZIn1   = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
2093   fHisChi2ZOut1  = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
2094   fHisChi2ZIn2   = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
2095   fHisChi2ZOut2  = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
2096   fHisChi2ZIO1   = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
2097   //
2098   // Fit
2099   //
2100   //
2101   //   TH2F           *fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
2102   //   TH2F           *fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
2103   //   TH2F           *fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
2104   //   TH2F           *fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
2105   //   TH2F           *fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
2106   //   TH2F           *fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
2107   //   TH2F           *fHisPy2vP2In;   //-> Curv  P2inner - parabola
2108   //   TH2F           *fHisPy2vP2Out;  //-> Curv  P2outer - parabola
2109   //   TH2F           *fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
2110   fHisPy1vP0   = new TH2F("HisPy1vP0",   "HisPy1vP0",336,0,336,500,-0.3,0.3);
2111   fHisPy2vP0   = new TH2F("HisPy2vP0",   "HisPy2vP0",336,0,336,500,-0.3,0.3);
2112   fHisPy3vP0   = new TH2F("HisPy3vP0",   "HisPy3vP0",336,0,336,500,-0.3,0.3);
2113   fHisPy1vP1   = new TH2F("HisPy1vP1",   "HisPy1vP1",336,0,336,500,-0.01,0.01);
2114   fHisPy2vP1   = new TH2F("HisPy2vP1",   "HisPy2vP1",336,0,336,500,-0.01,0.01);
2115   fHisPy3vP1   = new TH2F("HisPy3vP1",   "HisPy3vP1",336,0,336,500,-0.01,0.01);
2116   fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
2117   fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
2118   fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
2119   //
2120   //
2121   //   TH2F           *fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
2122   //   TH2F           *fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
2123   //   TH2F           *fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
2124   //   TH2F           *fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
2125   //   TH2F           *fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
2126   //   TH2F           *fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
2127   //   TH2F           *fHisPz2vP2In;   //-> Curv  P2inner - parabola
2128   //   TH2F           *fHisPz2vP2Out;  //-> Curv  P2outer - parabola
2129   //   TH2F           *fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
2130   fHisPz1vP0   = new TH2F("HisPz1vP0",   "HisPz1vP0",336,0,336,500,-0.3,0.3);
2131   fHisPz2vP0   = new TH2F("HisPz2vP0",   "HisPz2vP0",336,0,336,500,-0.3,0.3);
2132   fHisPz3vP0   = new TH2F("HisPz3vP0",   "HisPz3vP0",336,0,336,500,-0.3,0.3);
2133   fHisPz1vP1   = new TH2F("HisPz1vP1",   "HisPz1vP1",336,0,336,500,-0.01,0.01);
2134   fHisPz2vP1   = new TH2F("HisPz2vP1",   "HisPz2vP1",336,0,336,500,-0.01,0.01);
2135   fHisPz3vP1   = new TH2F("HisPz3vP1",   "HisPz3vP1",336,0,336,500,-0.01,0.01);
2136   fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
2137   fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
2138   fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
2139
2140   fHisNclIn->SetDirectory(0);      //->Number of clusters inner
2141   fHisNclOut->SetDirectory(0);     //->Number of clusters outer
2142   fHisNclIO->SetDirectory(0);      //->Number of cluster inner outer
2143   fHisLclIn->SetDirectory(0);      //->Level arm inner
2144   fHisLclOut->SetDirectory(0);     //->Level arm outer
2145   fHisLclIO->SetDirectory(0);      //->Number of cluster inner outer
2146   fHisdEdx->SetDirectory(0);      //->Number of cluster inner outer
2147   fHisdZfit->SetDirectory(0);      //->Number of cluster inner outer
2148   //
2149   //
2150   fHisChi2YIn1->SetDirectory(0);      //->chi2 y inner - line
2151   fHisChi2YOut1->SetDirectory(0);     //->chi2 y inner - line
2152   fHisChi2YIn2->SetDirectory(0);      //->chi2 y inner - parabola
2153   fHisChi2YOut2->SetDirectory(0);     //->chi2 y inner - parabola
2154   fHisChi2YIO1->SetDirectory(0);      //->chi2 y IO    - common
2155   fHisChi2ZIn1->SetDirectory(0);      //->chi2 z inner - line
2156   fHisChi2ZOut1->SetDirectory(0);     //->chi2 z inner - line
2157   fHisChi2ZIn2->SetDirectory(0);      //->chi2 z inner - parabola
2158   fHisChi2ZOut2->SetDirectory(0);     //->chi2 z inner - parabola
2159   fHisChi2ZIO1->SetDirectory(0);      //->chi2 z IO    - common
2160   //
2161   //
2162   fHisPy1vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - line
2163   fHisPy2vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - parabola
2164   fHisPy3vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - common parabola
2165   fHisPy1vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - line
2166   fHisPy2vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - parabola
2167   fHisPy3vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - common parabola
2168   fHisPy2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
2169   fHisPy2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
2170   fHisPy3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
2171   //
2172   //
2173   fHisPz1vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - line
2174   fHisPz2vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - parabola
2175   fHisPz3vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - common parabola
2176   fHisPz1vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - line
2177   fHisPz2vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - parabola
2178   fHisPz3vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - common parabola
2179   fHisPz2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
2180   fHisPz2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
2181   fHisPz3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
2182
2183
2184
2185 }
2186
2187 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
2188   //
2189   // Merge content of histograms
2190   //
2191   // Only first histogram is checked - all other should be the same
2192   if (!laser->fHisNclIn) return;  // empty histograms
2193   if (!fHisNclIn) MakeFitHistos();
2194   //
2195
2196   fHisNclIn->Add(laser->fHisNclIn  );      //->Number of clusters inner
2197   fHisNclOut->Add(laser->fHisNclOut  );     //->Number of clusters outer
2198   fHisNclIO->Add(laser->fHisNclIO  );      //->Number of cluster inner outer
2199   fHisLclIn->Add(laser->fHisLclIn  );      //->Level arm inner
2200   fHisLclOut->Add(laser->fHisLclOut  );     //->Level arm outer
2201   fHisLclIO->Add(laser->fHisLclIO  );      //->Number of cluster inner outer
2202   fHisdEdx->Add(laser->fHisdEdx  );      //->dedx
2203   fHisdZfit->Add(laser->fHisdZfit  );      //->dedx
2204   //
2205   //
2206   fHisChi2YIn1->Add(laser->fHisChi2YIn1  );      //->chi2 y inner - line
2207   fHisChi2YOut1->Add(laser->fHisChi2YOut1  );     //->chi2 y inner - line
2208   fHisChi2YIn2->Add(laser->fHisChi2YIn2  );      //->chi2 y inner - parabola
2209   fHisChi2YOut2->Add(laser->fHisChi2YOut2  );     //->chi2 y inner - parabola
2210   fHisChi2YIO1->Add(laser->fHisChi2YIO1  );      //->chi2 y IO    - common
2211   fHisChi2ZIn1->Add(laser->fHisChi2ZIn1  );      //->chi2 z inner - line
2212   fHisChi2ZOut1->Add(laser->fHisChi2ZOut1  );     //->chi2 z inner - line
2213   fHisChi2ZIn2->Add(laser->fHisChi2ZIn2  );      //->chi2 z inner - parabola
2214   fHisChi2ZOut2->Add(laser->fHisChi2ZOut2  );     //->chi2 z inner - parabola
2215   fHisChi2ZIO1->Add(laser->fHisChi2ZIO1  );      //->chi2 z IO    - common
2216   //
2217   //
2218   fHisPy1vP0->Add(laser->fHisPy1vP0  );     //-> delta y   P0outer-P0inner - line
2219   fHisPy2vP0->Add(laser->fHisPy2vP0  );     //-> delta y   P0outer-P0inner - parabola
2220   fHisPy3vP0->Add(laser->fHisPy3vP0  );     //-> delta y   P0outer-P0inner - common parabola
2221   fHisPy1vP1->Add(laser->fHisPy1vP1  );     //-> delta ky  P1outer-P1inner - line
2222   fHisPy2vP1->Add(laser->fHisPy2vP1  );     //-> delta ky  P1outer-P1inner - parabola
2223   fHisPy3vP1->Add(laser->fHisPy3vP1  );     //-> delta ky  P1outer-P1inner - common parabola
2224   fHisPy2vP2In->Add(laser-> fHisPy2vP2In  );   //-> Curv  P2inner - parabola
2225   fHisPy2vP2Out->Add(laser->fHisPy2vP2Out  );  //-> Curv  P2outer - parabola
2226   fHisPy3vP2IO->Add(laser->fHisPy3vP2IO  );   //-> Curv  P2outerinner - common parabola
2227   //
2228   //
2229   fHisPz1vP0->Add(laser->fHisPz1vP0  );     //-> delta z   P0outer-P0inner - line
2230   fHisPz2vP0->Add(laser->fHisPz2vP0  );     //-> delta z   P0outer-P0inner - parabola
2231   fHisPz3vP0->Add(laser->fHisPz3vP0  );     //-> delta z   P0outer-P0inner - common parabola
2232   fHisPz1vP1->Add(laser->fHisPz1vP1  );     //-> delta kz  P1outer-P1inner - line
2233   fHisPz2vP1->Add(laser->fHisPz2vP1  );     //-> delta kz  P1outer-P1inner - parabola
2234   fHisPz3vP1->Add(laser->fHisPz3vP1  );     //-> delta kz  P1outer-P1inner - common parabola
2235   fHisPz2vP2In->Add(laser->fHisPz2vP2In  );   //-> Curv  P2inner - parabola
2236   fHisPz2vP2Out->Add(laser->fHisPz2vP2Out  );  //-> Curv  P2outer - parabola
2237   fHisPz3vP2IO->Add(laser->fHisPz3vP2IO  );   //-> Curv  P2outerinner - common parabola
2238   //
2239   //
2240   //
2241
2242
2243
2244
2245 }
2246
2247
2248
2249
2250 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
2251   //
2252   // Dump fit information - collect information from the streamers
2253   //
2254   /*
2255     TChain * chainFit=0;
2256     TChain * chainTrack=0;
2257     TChain * chain=0;
2258     //
2259     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2260     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
2261     AliXRDPROOFtoolkit tool;
2262     chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
2263     chainTrack->Lookup();
2264     chainTrack->SetProof(kTRUE);
2265     chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
2266     chainDrift->Lookup();
2267     chainDrift->SetProof(kTRUE);
2268
2269     chain = tool.MakeChain("laser.txt","Residuals",0,10200);
2270     chain->Lookup();
2271     chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
2272     chainFit->Lookup();
2273     chainFit->SetProof(kTRUE);
2274     chain->SetProof(kTRUE);
2275     AliTPCLaserTrack::LoadTracks();
2276     //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
2277
2278   */
2279   //
2280   // Fit cuts
2281   //
2282   TCut cutP3("abs(Tr.fP[3])<0.1");
2283   TCut cutP4("abs(Tr.fP[4])<0.5");
2284   TCut cutPx = cutP3+cutP4;
2285   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
2286   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
2287   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
2288   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
2289   //
2290   TCut cutdEdx("sqrt(dEdx)>3");
2291   TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
2292   TCut cutN("nclO>20&&nclI>20");
2293   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
2294   //
2295   // Cluster cuts
2296   //
2297   TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
2298   TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
2299   TCut cutClX("abs(Cl[].fX)>10");
2300   TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
2301   TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
2302   TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
2303   TCut cutQ("sqrt(Cl[].fMax)>4");
2304   TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
2305
2306
2307   TH1F * phisAl     = 0;
2308   TH1F * phisAccept = 0;
2309   TH1F * phisOut    = 0;
2310   TProfile * pdEdx  = 0;
2311
2312   TProfile * pP0    = 0;
2313   TProfile * pP1    = 0;
2314   TProfile * pP2    = 0;
2315   TProfile * pP3    = 0;
2316   TProfile * pP4    = 0;
2317   //
2318   TProfile * pNclI  = 0;
2319   TProfile * pNclO  = 0;
2320   //
2321   TProfile * pchi2YIn   =0;
2322   TProfile * pchi2ZIn   =0;
2323   TProfile * pchi2YOut  =0;
2324   TProfile * pchi2ZOut  =0;
2325   TProfile * pchi2YInOut =0;
2326   TProfile * pchi2ZInOut =0;;
2327   // laser counters
2328   chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
2329   phisAl = (TH1F*)gROOT->FindObject("hisAl");
2330   chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
2331   phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
2332   chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
2333   phisOut = (TH1F*)gROOT->FindObject("hisOut");
2334   //
2335   chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
2336   pdEdx   = (TProfile*)gROOT->FindObject("hdEdx");
2337   // track param
2338   //
2339   chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2340   pP0   = (TProfile*)gROOT->FindObject("hP0");
2341   chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2342   pP1   = (TProfile*)gROOT->FindObject("hP1");
2343   chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2344   pP2   = (TProfile*)gROOT->FindObject("hP2");
2345   chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2346   pP3   = (TProfile*)gROOT->FindObject("hP3");
2347   chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2348   pP4   = (TProfile*)gROOT->FindObject("hP4");
2349   //
2350   chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
2351   pNclI   = (TProfile*)gROOT->FindObject("hNclI");
2352   chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
2353   pNclO   = (TProfile*)gROOT->FindObject("hNclO");
2354   //
2355   //
2356   chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
2357   pchi2YIn   = (TProfile*)gROOT->FindObject("hChi2YIn");
2358   chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
2359   pchi2YOut   = (TProfile*)gROOT->FindObject("hChi2YOut");
2360   chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
2361   pchi2YInOut   = (TProfile*)gROOT->FindObject("hChi2YInOut");
2362   chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
2363   pchi2ZIn   = (TProfile*)gROOT->FindObject("hChi2ZIn");
2364   chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
2365   pchi2ZOut   = (TProfile*)gROOT->FindObject("hChi2ZOut");
2366   chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
2367   pchi2ZInOut   = (TProfile*)gROOT->FindObject("hChi2ZInOut");
2368   //
2369   // second derivatives
2370   //
2371   TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
2372   chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
2373   TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
2374   chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
2375   TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
2376   chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
2377   //
2378   phisPy2In->FitSlicesY();
2379   TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
2380   TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
2381   TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
2382   //
2383   phisPy2Out->FitSlicesY();
2384   TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
2385   TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
2386   TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
2387   //
2388   phisPy2InOut->FitSlicesY();
2389   TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
2390   TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
2391   TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
2392   //
2393   TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
2394   chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
2395   TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
2396   chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
2397   TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
2398   chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
2399   //
2400   phisPz2In->FitSlicesY();
2401   TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
2402   TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
2403   TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
2404   //
2405   phisPz2Out->FitSlicesY();
2406   TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
2407   TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
2408   TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
2409   //
2410   phisPz2InOut->FitSlicesY();
2411   TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
2412   TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
2413   TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
2414   //
2415   //
2416   //
2417
2418
2419   {
2420     TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
2421     for (Int_t ilaser=0; ilaser<336; ilaser++){
2422       Float_t all=phisAl->GetBinContent(ilaser+1);
2423       Float_t accept=phisAccept->GetBinContent(ilaser+1);
2424       Float_t out=phisOut->GetBinContent(ilaser+1);
2425       Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
2426       Float_t mP0 = pP0->GetBinContent(ilaser+1);
2427       Float_t mP1 = pP1->GetBinContent(ilaser+1);
2428       Float_t mP2 = pP2->GetBinContent(ilaser+1);
2429       Float_t mP3 = pP3->GetBinContent(ilaser+1);
2430       Float_t mP4 = pP4->GetBinContent(ilaser+1);
2431
2432
2433       Float_t nclI  = pNclI->GetBinContent(ilaser+1);
2434       Float_t nclO  = pNclO->GetBinContent(ilaser+1);
2435       //
2436       Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
2437       Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
2438       Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
2439       Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
2440       Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
2441       Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
2442       //
2443       Float_t entriesPy2In  = phisPy2InEntries->GetBinContent(ilaser+1);
2444       Float_t meanPy2In     = phisPy2InMean->GetBinContent(ilaser+1);
2445       Float_t sigmaPy2In    = phisPy2InSigma->GetBinContent(ilaser+1);
2446       //
2447       Float_t entriesPy2Out  = phisPy2OutEntries->GetBinContent(ilaser+1);
2448       Float_t meanPy2Out     = phisPy2OutMean->GetBinContent(ilaser+1);
2449       Float_t sigmaPy2Out    = phisPy2OutSigma->GetBinContent(ilaser+1);
2450       //
2451       Float_t entriesPy2InOut  = phisPy2InOutEntries->GetBinContent(ilaser+1);
2452       Float_t meanPy2InOut     = phisPy2InOutMean->GetBinContent(ilaser+1);
2453       Float_t sigmaPy2InOut    = phisPy2InOutSigma->GetBinContent(ilaser+1);
2454       //
2455       Float_t entriesPz2In  = phisPz2InEntries->GetBinContent(ilaser+1);
2456       Float_t meanPz2In     = phisPz2InMean->GetBinContent(ilaser+1);
2457       Float_t sigmaPz2In    = phisPz2InSigma->GetBinContent(ilaser+1);
2458       //
2459       Float_t entriesPz2Out  = phisPz2OutEntries->GetBinContent(ilaser+1);
2460       Float_t meanPz2Out     = phisPz2OutMean->GetBinContent(ilaser+1);
2461       Float_t sigmaPz2Out    = phisPz2OutSigma->GetBinContent(ilaser+1);
2462       //
2463       Float_t entriesPz2InOut  = phisPz2InOutEntries->GetBinContent(ilaser+1);
2464       Float_t meanPz2InOut     = phisPz2InOutMean->GetBinContent(ilaser+1);
2465       Float_t sigmaPz2InOut    = phisPz2InOutSigma->GetBinContent(ilaser+1);
2466
2467       AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
2468       (*pcstream)<<"Scan"<<
2469         "Run="<<id<<
2470         "LTr.="<<ltrp<<
2471         "all="<<all<<
2472         "accept="<<accept<<
2473         "out="<<out<<
2474         "sdedx="<<sdedx<<
2475         "mP0="<<mP0<<
2476         "mP1="<<mP1<<
2477         "mP2="<<mP2<<
2478         "mP3="<<mP3<<
2479         "mP4="<<mP4<<
2480         "nclI="<<nclI<<
2481         "nclO="<<nclO<<
2482         "chi2YIn="<<chi2YIn<<
2483         "chi2YOut="<<chi2YOut<<
2484         "chi2YInOut="<<chi2YInOut<<
2485         "chi2ZIn="<<chi2ZIn<<
2486         "chi2ZOut="<<chi2ZOut<<
2487         "chi2ZInOut="<<chi2ZInOut<<
2488         //
2489         "nPy2In="<<entriesPy2In<<
2490         "mPy2In="<<meanPy2In<<
2491         "sPy2In="<<sigmaPy2In<<
2492         //
2493         "nPy2Out="<<entriesPy2Out<<
2494         "mPy2Out="<<meanPy2Out<<
2495         "sPy2Out="<<sigmaPy2Out<<
2496         //
2497         "nPy2InOut="<<entriesPy2InOut<<
2498         "mPy2InOut="<<meanPy2InOut<<
2499         "sPy2InOut="<<sigmaPy2InOut<<
2500         //
2501         "nPz2In="<<entriesPz2In<<
2502         "mPz2In="<<meanPz2In<<
2503         "sPz2In="<<sigmaPz2In<<
2504         //
2505         "nPz2Out="<<entriesPz2Out<<
2506         "mPz2Out="<<meanPz2Out<<
2507         "sPz2Out="<<sigmaPz2Out<<
2508         //
2509         "nPz2InOut="<<entriesPz2InOut<<
2510         "mPz2InOut="<<meanPz2InOut<<
2511         "sPz2InOut="<<sigmaPz2InOut<<
2512         "\n";
2513     }
2514
2515     delete pcstream;
2516   }
2517   /*
2518     TFile f("vscan.root");
2519    */
2520
2521   /*
2522     pad binning effect
2523    chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
2524    //
2525    chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
2526    //
2527
2528 chain->Draw("Cl.fY-TrYpol1.fElements-AliTPCClusterParam::SPosCorrection(0,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
2529
2530
2531 chain->Draw("Cl[].fZ-TrZpol1.fElements-0*AliTPCClusterParam::SPosCorrection(1,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fTimeBin-int(Cl[].fTimeBin)",cutA+cutCl+"Cl[].fZ>0","prof",10000)
2532
2533   */
2534
2535
2536
2537
2538
2539   /*
2540   // check edge effects
2541   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
2542   //
2543   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
2544
2545   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
2546
2547
2548
2549   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
2550   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
2551
2552 */
2553
2554
2555
2556
2557
2558   /*
2559     Edge y effect
2560
2561     dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
2562
2563
2564     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):pi/18-abs(Cl.fY/Cl.fX)>>hisYdphi(100,0,0.03)",""+cutA+cutCl,"prof",10000)
2565
2566     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdy(100,0,5)",""+cutA+cutCl,"prof",10000)
2567
2568
2569
2570
2571
2572     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyIROC(100,0,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
2573
2574     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyOROC(100,0,5)","Cl.fDetector>36"+cutA+cutCl,"prof",100000)
2575
2576
2577
2578     chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>his(100,-5,5)",""+cutA+cutCl,"prof",100000)
2579
2580     chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>hisdyInner(100,-5,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
2581
2582
2583
2584 */
2585
2586
2587 /*
2588
2589 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
2590
2591 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
2592
2593
2594
2595 chainFit->Draw("LTr.fId","nclI>10",100000)
2596
2597 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
2598
2599 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
2600
2601 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
2602
2603 */
2604
2605 }
2606
2607
2608
2609
2610
2611
2612 /*
2613  gSystem->Load("libSTAT.so")
2614  TStatToolkit toolkit;
2615  Double_t chi2;
2616  TVectorD fitParam;
2617  TMatrixD covMatrix;
2618  Int_t npoints;
2619
2620  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
2621
2622
2623 TString fstring="";
2624 //
2625 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
2626 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
2627 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
2628 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4
2629 //
2630 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
2631 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
2632 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
2633 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
2634 //
2635 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
2636 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
2637 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
2638 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
2639
2640
2641
2642
2643  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
2644
2645  treeT->SetAlias("fit",strq0->Data());
2646
2647
2648  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
2649
2650  treeT->SetAlias("fitP",strqP->Data());
2651
2652
2653  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
2654  treeT->SetAlias("fitD",strqDrift->Data());
2655
2656
2657 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
2658 {
2659 for (Int_t i=0; i<6;i++){
2660 treeT->SetLineColor(i+2);
2661 treeT->SetMarkerSize(1);
2662 treeT->SetMarkerStyle(22+i);
2663 treeT->SetMarkerColor(i+2);
2664
2665 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
2666 }
2667 }
2668  */
2669
2670
2671
2672 /*
2673   TTree * tree = (TTree*)f.Get("FitModels");
2674
2675   TEventList listLFit0("listLFit0","listLFit0");
2676   TEventList listLFit1("listLFit1","listLFit1");
2677   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
2678   tree->SetEventList(&listLFit0);
2679
2680
2681
2682
2683   gSystem->Load("libSTAT.so")
2684   TStatToolkit toolkit;
2685   Double_t chi2;
2686   TVectorD fitParam;
2687   TMatrixD covMatrix;
2688   Int_t npoints;
2689
2690   chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
2691   chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
2692
2693
2694   TString fstring="";
2695   fstring+="cos(dp)++";
2696   fstring+="sin(dp)++";
2697   fstring+="cos(dt)++";
2698   fstring+="sin(dt)++";
2699
2700   TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
2701
2702
2703
2704 */
2705
2706
2707
2708 /*
2709   Edge effects
2710   //
2711   //
2712   //
2713   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2714   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
2715   AliXRDPROOFtoolkit tool;
2716   TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
2717   chainTrack->Lookup();
2718   chainTrack->SetProof(kTRUE);
2719
2720   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
2721   chain->Lookup();
2722   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
2723   chainFit->Lookup();
2724   chainFit->SetProof(kTRUE);
2725   chain->SetProof(kTRUE);
2726   //
2727   // Fit cuts
2728   //
2729   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
2730   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
2731   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
2732   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
2733   //
2734   TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
2735   TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
2736   TCut cutN("nclO>20&&nclI>20");
2737   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
2738   //
2739   // Cluster cuts
2740   //
2741   TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
2742   TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
2743   TCut cutClX("abs(Cl.fX)>10");
2744   TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
2745   TCut cutCl=cutClY+cutClZ+cutClX;
2746
2747
2748   // check edge effects
2749   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
2750   //
2751   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
2752
2753   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
2754
2755
2756
2757   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
2758   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
2759
2760 */
2761