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