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 ~/NimStyle.C
  38   gSystem->Load("libANALYSIS");
  39   gSystem->Load("libTPCcalib");
  40   TFile fcalib("CalibObjectsTrain2.root");
  41   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
  42   laser->DumpMeanInfo(0)
  43   TFile fmean("laserMean.root")
  44   //
  45   //  laser track clasification;
  46   //
  47   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
  48   TCut cutN("cutN","fTPCncls>70");
  49   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
  50   TCut cutA = cutT+cutPt+cutP;
  51   TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
  52
  53   //
  54   //
  55   // Analyze  LASER scan
  56   //
  57
  58   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  59   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  60   AliXRDPROOFtoolkit tool;
  61   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
  62   chain->Lookup();
  63   AliTPCcalibLaser::DumpScanInfo(chain)
  64   TFile fscan("laserScan.root")
  65   TTree * treeT = (TTree*)fscan.Get("Mean")
  66   //
  67   // Analyze laser
  68   //
  69   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  70   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  71   AliXRDPROOFtoolkit tool;
  72   AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1)
  73   TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
  74   chainDrift->Lookup();
  75   TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
  76   chainDriftN->Lookup();
  77
  78
  79   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
  80   chain->Lookup();
  81   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
  82   chainFit->Lookup();
  83   TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
  84   chainTrack->Lookup();
  85
  86 */
  87
  88
  89
  90 #include "TLinearFitter.h"
  91 #include "AliTPCcalibLaser.h"
  92 #include "AliExternalTrackParam.h"
  93 #include "AliESDEvent.h"
  94 #include "AliESDfriend.h"
  95 #include "AliESDtrack.h"
  96 #include "AliTPCTracklet.h"
  97 #include "TH1D.h"
  98 #include "TH1F.h"
  99 #include "TProfile.h"
 100 #include "TVectorD.h"
 101 #include "TTreeStream.h"
 102 #include "TFile.h"
 103 #include "TF1.h"
 104 #include "TGraphErrors.h"
 105 #include "AliTPCclusterMI.h"
 106 #include "AliTPCseed.h"
 107 #include "AliTracker.h"
 108 #include "AliLog.h"
 109 #include "TClonesArray.h"
 110 #include "TPad.h"
 111 #include "TSystem.h"
 112 #include "TCut.h"
 113 #include "TChain.h"
 114 #include "TH2F.h"
 115 #include "TStatToolkit.h"
 116 #include "TROOT.h"
 117
 118
 119 #include "TTreeStream.h"
 120 #include <iostream>
 121 #include <sstream>
 122 #include "AliTPCLaserTrack.h"
 123 #include "AliTPCcalibDB.h"
 124 #include "AliTPCParam.h"
 125 #include "AliTPCParamSR.h"
 126 #include "TTimeStamp.h"
 127 #include "AliDCSSensorArray.h"
 128 #include "AliDCSSensor.h"
 129 #include "AliGRPObject.h"
 130
 131 using namespace std;
 132
 133 ClassImp(AliTPCcalibLaser)
 134
 135 AliTPCcalibLaser::AliTPCcalibLaser():
 136   AliTPCcalibBase(),
 137   fESD(0),
 138   fESDfriend(0),
 139   fNtracks(0),
 140   fTracksMirror(336),
 141   fTracksEsd(336),
 142   fTracksEsdParam(336),
 143   fTracksTPC(336),
 144   fFullCalib(kTRUE),
 145   fDeltaZ(336),
 146   fDeltaP3(336),
 147   fDeltaP4(336),
 148   fDeltaPhi(336),
 149   fDeltaPhiP(336),
 150   fSignals(336),
 151   //
 152   fHisLaser(0),      //  N dim histogram of laser
 153   fHisNclIn(0),      //->Number of clusters inner
 154   fHisNclOut(0),     //->Number of clusters outer
 155   fHisNclIO(0),      //->Number of cluster inner outer
 156   fHisLclIn(0),      //->Level arm inner
 157   fHisLclOut(0),     //->Level arm outer
 158   fHisLclIO(0),      //->Number of cluster inner outer
 159   fHisdEdx(0),     //->dEdx histo
 160   fHisdZfit(0),     //->distance to the mirror after linear fit
 161   //
 162   //
 163   fHisChi2YIn1(0),      //->chi2 y inner - line
 164   fHisChi2YOut1(0),     //->chi2 y inner - line
 165   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 166   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 167   fHisChi2YIO1(0),      //->chi2 y IO    - common
 168   fHisChi2ZIn1(0),      //->chi2 z inner - line
 169   fHisChi2ZOut1(0),     //->chi2 z inner - line
 170   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 171   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 172   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 173   //
 174   //
 175   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 176   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 177   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 178   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 179   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 180   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 181   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 182   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 183   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 184   //
 185   //
 186   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 187   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 188   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 189   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 190   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 191   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 192   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 193   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 194   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 195   //
 196   fDeltaYres(336),   //->2D histo of residuals
 197   fDeltaZres(336),   //->2D histo fo residuals
 198   fDeltaYres2(336),   //->2D histo of residuals
 199   fDeltaZres2(336),   //->2D histo fo residuals
 200   fDeltaYresAbs(336), //->2D histo of residuals
 201   fHisYAbsErrors(0),  //-> total errors per beam in the abs res y analysis
 202   fDeltaZresAbs(336), //->2D histo of residuals
 203   fHisZAbsErrors(0),  //-> total errors per beam in the abs res z analysis
 204   //fDeltaYres3(336),   //->2D histo of residuals
 205   //fDeltaZres3(336),   //->2D histo fo residuals
 206   fFitAside(new TVectorD(5)),
 207   fFitCside(new TVectorD(5)),
 208   fFitACside(new TVectorD(6)),
 209   fEdgeXcuts(3),
 210   fEdgeYcuts(3),
 211   fNClCuts(5),
 212   fNcuts(0),
 213   fBeamSectorOuter(336),
 214   fBeamSectorInner(336),
 215   fBeamOffsetYOuter(336),
 216   fBeamSlopeYOuter(336),
 217   fBeamOffsetYInner(336),
 218   fBeamSlopeYInner(336),
 219   fBeamOffsetZOuter(336),
 220   fBeamSlopeZOuter(336),
 221   fBeamOffsetZInner(336),
 222   fBeamSlopeZInner(336),
 223   fInverseSlopeZ(kTRUE),
 224   fUseFixedDriftV(0),
 225   fFixedFitAside0(0.0),
 226   fFixedFitAside1(1.0),
 227   fFixedFitCside0(0.0),
 228   fFixedFitCside1(1.0)
 229 {
 230   //
 231   // Constructor
 232   //
 233   fTracksEsdParam.SetOwner(kTRUE);
 234 }
 235
 236 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
 237   AliTPCcalibBase(),
 238   fESD(0),
 239   fESDfriend(0),
 240   fNtracks(0),
 241   fTracksMirror(336),
 242   fTracksEsd(336),
 243   fTracksEsdParam(336),
 244   fTracksTPC(336),
 245   fFullCalib(full),
 246   //
 247   fDeltaZ(336),          // array of histograms of delta z for each track
 248   fDeltaP3(336),          // array of histograms of delta z for each track
 249   fDeltaP4(336),          // array of histograms of P3 for each track
 250   fDeltaPhi(336),          // array of histograms of P4 for each track
 251   fDeltaPhiP(336),          // array of histograms of delta z for each track
 252   fSignals(336),           // array of dedx signals
 253   //
 254   //
 255   fHisLaser(0),      //  N dim histogram of laser
 256   fHisNclIn(0),      //->Number of clusters inner
 257   fHisNclOut(0),     //->Number of clusters outer
 258   fHisNclIO(0),      //->Number of cluster inner outer
 259   fHisLclIn(0),      //->Level arm inner
 260   fHisLclOut(0),     //->Level arm outer
 261   fHisLclIO(0),      //->Number of cluster inner outer
 262   fHisdEdx(0), //->dEdx histo
 263   fHisdZfit(0),     //->distance to the mirror after linear fit
 264   //
 265   //
 266   fHisChi2YIn1(0),      //->chi2 y inner - line
 267   fHisChi2YOut1(0),     //->chi2 y inner - line
 268   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 269   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 270   fHisChi2YIO1(0),      //->chi2 y IO    - common
 271   fHisChi2ZIn1(0),      //->chi2 z inner - line
 272   fHisChi2ZOut1(0),     //->chi2 z inner - line
 273   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 274   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 275   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 276   //
 277   //
 278   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 279   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 280   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 281   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 282   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 283   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 284   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 285   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 286   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 287   //
 288   //
 289   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 290   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 291   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 292   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 293   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 294   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 295   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 296   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 297   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 298   //
 299   //
 300   //
 301   fDeltaYres(336),
 302   fDeltaZres(336),
 303   fDeltaYres2(336),
 304   fDeltaZres2(336),
 305   fDeltaYresAbs(336),
 306   fHisYAbsErrors(0),
 307   fDeltaZresAbs(336),
 308   fHisZAbsErrors(0),
 309   //  fDeltaYres3(336),
 310   //fDeltaZres3(336),
 311   fFitAside(new TVectorD(5)),        // drift fit - A side
 312   fFitCside(new TVectorD(5)),        // drift fit - C- side
 313   fFitACside(new TVectorD(6)),        // drift fit - AC- side
 314   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 315   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 316   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 317   fNcuts(0),           // number of cuts
 318   fBeamSectorOuter(336),
 319   fBeamSectorInner(336),
 320   fBeamOffsetYOuter(336),
 321   fBeamSlopeYOuter(336),
 322   fBeamOffsetYInner(336),
 323   fBeamSlopeYInner(336),
 324   fBeamOffsetZOuter(336),
 325   fBeamSlopeZOuter(336),
 326   fBeamOffsetZInner(336),
 327   fBeamSlopeZInner(336),
 328   fInverseSlopeZ(kTRUE),
 329   fUseFixedDriftV(0),
 330   fFixedFitAside0(0.0),
 331   fFixedFitAside1(1.0),
 332   fFixedFitCside0(0.0),
 333   fFixedFitCside1(1.0)
 334 {
 335   SetName(name);
 336   SetTitle(title);
 337   //
 338   // Constructor
 339   //
 340   fTracksEsdParam.SetOwner(kTRUE);
 341 }
 342
 343 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
 344   AliTPCcalibBase(calibLaser),
 345   fESD(0),
 346   fESDfriend(0),
 347   fNtracks(0),
 348   fTracksMirror(336),
 349   fTracksEsd(336),
 350   fTracksEsdParam(336),
 351   fTracksTPC(336),
 352   fFullCalib(calibLaser.fFullCalib),
 353   //
 354   fDeltaZ(calibLaser.fDeltaZ),          // array of histograms of delta z for each track
 355   fDeltaP3(((calibLaser.fDeltaP3))),          // array of histograms of delta z for each track
 356   fDeltaP4(((calibLaser.fDeltaP4))),          // array of histograms of P3 for each track
 357   fDeltaPhi(((calibLaser.fDeltaPhi))),          // array of histograms of P4 for each track
 358   fDeltaPhiP(((calibLaser.fDeltaPhiP))),          // array of histograms of delta z for each track
 359   fSignals(((calibLaser.fSignals))),           // array of dedx signals
 360   //
 361   //
 362   fHisLaser(0),      //  N dim histogram of laser
 363   fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))),      //->Number of clusters inner
 364   fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))),     //->Number of clusters outer
 365   fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))),      //->Number of cluster inner outer
 366   fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))),      //->Level arm inner
 367   fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))),     //->Level arm outer
 368   fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))),      //->Number of cluster inner outer
 369   fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))),      //
 370   fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))),     //->distance to the mirror after linear fit
 371   //
 372   //
 373   fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))),      //->chi2 y inner - line
 374   fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))),     //->chi2 y inner - line
 375   fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))),      //->chi2 y inner - parabola
 376   fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))),     //->chi2 y inner - parabola
 377   fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))),      //->chi2 y IO    - common
 378   fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))),      //->chi2 z inner - line
 379   fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))),     //->chi2 z inner - line
 380   fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))),      //->chi2 z inner - parabola
 381   fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))),     //->chi2 z inner - parabola
 382   fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))),      //->chi2 z IO    - common
 383   //
 384   //
 385   fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))),     //-> delta y   P0outer-P0inner - line
 386   fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))),     //-> delta y   P0outer-P0inner - parabola
 387   fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))),     //-> delta y   P0outer-P0inner - common parabola
 388   fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))),     //-> delta ky  P1outer-P1inner - line
 389   fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))),     //-> delta ky  P1outer-P1inner - parabola
 390   fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))),     //-> delta ky  P1outer-P1inner - common parabola
 391   fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))),   //-> Curv  P2inner - parabola
 392   fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))),  //-> Curv  P2outer - parabola
 393   fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))),   //-> Curv  P2outerinner - common parabola
 394   //
 395   //
 396   fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))),     //-> delta z   P0outer-P0inner - line
 397   fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))),     //-> delta z   P0outer-P0inner - parabola
 398   fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))),     //-> delta z   P0outer-P0inner - common parabola
 399   fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))),     //-> delta kz  P1outer-P1inner - line
 400   fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))),     //-> delta kz  P1outer-P1inner - parabola
 401   fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))),     //-> delta kz  P1outer-P1inner - common parabola
 402   fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))),   //-> Curv  P2inner - parabola
 403   fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))),  //-> Curv  P2outer - parabola
 404   fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))),   //-> Curv  P2outerinner - common parabola
 405   //
 406   //
 407   fDeltaYres(((calibLaser.fDeltaYres))),
 408   fDeltaZres(((calibLaser.fDeltaZres))),
 409   fDeltaYres2(((calibLaser.fDeltaYres))),
 410   fDeltaZres2(((calibLaser.fDeltaZres))),
 411   fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
 412   fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
 413   fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
 414   fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
 415   //  fDeltaYres3(((calibLaser.fDeltaYres))),
 416   //fDeltaZres3(((calibLaser.fDeltaZres))),
 417   fFitAside(new TVectorD(5)),        // drift fit - A side
 418   fFitCside(new TVectorD(5)),        // drift fit - C- side
 419   fFitACside(new TVectorD(6)),        // drift fit - C- side
 420   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 421   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 422   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 423   fNcuts(0),           // number of cuts
 424   fBeamSectorOuter(336),
 425   fBeamSectorInner(336),
 426   fBeamOffsetYOuter(336),
 427   fBeamSlopeYOuter(336),
 428   fBeamOffsetYInner(336),
 429   fBeamSlopeYInner(336),
 430   fBeamOffsetZOuter(336),
 431   fBeamSlopeZOuter(336),
 432   fBeamOffsetZInner(336),
 433   fBeamSlopeZInner(336),
 434   fInverseSlopeZ(calibLaser.fInverseSlopeZ),
 435   fUseFixedDriftV(calibLaser.fUseFixedDriftV),
 436   fFixedFitAside0(calibLaser.fFixedFitAside0),
 437   fFixedFitAside1(calibLaser.fFixedFitAside1),
 438   fFixedFitCside0(calibLaser.fFixedFitCside0),
 439   fFixedFitCside1(calibLaser.fFixedFitCside1)
 440 {
 441   //
 442   // copy constructor
 443   //
 444 }
 445
 446
 447
 448 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
 449   //
 450   // assgnment operator
 451   //
 452   if (this != &calibLaser) {
 453     new (this) AliTPCcalibLaser(calibLaser);
 454   }
 455   return *this;
 456
 457 }
 458
 459
 460
 461
 462 AliTPCcalibLaser::~AliTPCcalibLaser() {
 463   //
 464   // destructor
 465   //
 466   if ( fHisNclIn){
 467     delete fHisLaser;      //->
 468     delete fHisNclIn;      //->Number of clusters inner
 469     delete fHisNclOut;     //->Number of clusters outer
 470     delete fHisNclIO;      //->Number of cluster inner outer
 471     delete fHisLclIn;      //->Level arm inner
 472     delete fHisLclOut;     //->Level arm outer
 473     delete fHisLclIO;      //->Number of cluster inner outer
 474     delete fHisdEdx;
 475     delete fHisdZfit;
 476     //
 477     //
 478     delete fHisChi2YIn1;      //->chi2 y inner - line
 479     delete fHisChi2YOut1;     //->chi2 y inner - line
 480     delete fHisChi2YIn2;      //->chi2 y inner - parabola
 481     delete fHisChi2YOut2;     //->chi2 y inner - parabola
 482     delete fHisChi2YIO1;      //->chi2 y IO    - common
 483     delete fHisChi2ZIn1;      //->chi2 z inner - line
 484     delete fHisChi2ZOut1;     //->chi2 z inner - line
 485     delete fHisChi2ZIn2;      //->chi2 z inner - parabola
 486     delete fHisChi2ZOut2;     //->chi2 z inner - parabola
 487     delete fHisChi2ZIO1;      //->chi2 z IO    - common
 488     //
 489     //
 490     delete fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
 491     delete fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
 492     delete fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
 493     delete fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
 494     delete fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
 495     delete fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
 496     delete fHisPy2vP2In;   //-> Curv  P2inner - parabola
 497     delete fHisPy2vP2Out;  //-> Curv  P2outer - parabola
 498     delete fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
 499     //
 500     //
 501     delete fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
 502     delete fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
 503     delete fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
 504     delete fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
 505     delete fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
 506     delete fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
 507     delete fHisPz2vP2In;   //-> Curv  P2inner - parabola
 508     delete fHisPz2vP2Out;  //-> Curv  P2outer - parabola
 509     delete fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
 510   }
 511   //
 512   //
 513   //
 514   fDeltaZ.SetOwner();          //-> array of histograms of delta z for each track
 515   fDeltaP3.SetOwner();         //-> array of histograms of P3      for each track
 516   fDeltaP4.SetOwner();         //-> array of histograms of P4      for each track
 517   fDeltaPhi.SetOwner();        //-> array of histograms of delta z for each track
 518   fDeltaPhiP.SetOwner();       //-> array of histograms of delta z for each track
 519   fSignals.SetOwner();         //->Array of dedx signals
 520
 521   fDeltaZ.Delete();          //-> array of histograms of delta z for each track
 522   fDeltaP3.Delete();         //-> array of histograms of P3      for each track
 523   fDeltaP4.Delete();         //-> array of histograms of P4      for each track
 524   fDeltaPhi.Delete();        //-> array of histograms of delta z for each track
 525   fDeltaPhiP.Delete();       //-> array of histograms of delta z for each track
 526   fSignals.Delete();         //->Array of dedx signals
 527
 528   fDeltaYres.SetOwner();
 529   fDeltaYres.Delete();
 530   delete fHisYAbsErrors;
 531   fDeltaZres.SetOwner();
 532   fDeltaZres.Delete();
 533   delete fHisZAbsErrors;
 534   fDeltaYres2.SetOwner();
 535   fDeltaYres2.Delete();
 536   fDeltaZres2.SetOwner();
 537   fDeltaZres2.Delete();
 538
 539   fDeltaYresAbs.SetOwner();
 540   fDeltaYresAbs.Delete();
 541   fDeltaZresAbs.SetOwner();
 542   fDeltaZresAbs.Delete();
 543 }
 544
 545
 546
 547 void AliTPCcalibLaser::Process(AliESDEvent * event) {
 548   //
 549   //
 550   // Loop over tracks and call  Process function
 551   //
 552   const Int_t  kMinTracks=20;
 553   const Int_t  kMinClusters=40;
 554
 555   fESD = event;
 556   if (!fESD) {
 557     return;
 558   }
 559   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 560   if (!fESDfriend) {
 561     return;
 562   }
 563   if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
 564   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
 565   //
 566   // find CE background if present
 567   //
 568   if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
 569   TH1D hisCE("hhisCE","hhisCE",100,-100,100);
 570   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 571     AliESDtrack *track=fESD->GetTrack(i);
 572     if (!track) continue;
 573     hisCE.Fill(track->GetZ());
 574     hisCE.Fill(track->GetZ()+2);
 575     hisCE.Fill(track->GetZ()-2);
 576   }
 577   //
 578   //
 579
 580
 581   fTracksTPC.Clear();
 582   fTracksEsd.Clear();
 583   fTracksEsdParam.Delete();
 584   for (Int_t id=0; id<336;id++) {
 585     fCounter[id]=0;
 586     fClusterCounter[id]=0;
 587     fClusterSatur[id]=0;
 588   }
 589   //
 590   Int_t n=fESD->GetNumberOfTracks();
 591   Int_t counter=0;
 592   for (Int_t i=0;i<n;++i) {
 593     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
 594     if (!friendTrack) continue;
 595     AliESDtrack *track=fESD->GetTrack(i);
 596     Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
 597     if (binC>336) continue; //remove CE background
 598     TObject *calibObject=0;
 599     AliTPCseed *seed=0;
 600     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
 601       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
 602         break;
 603     if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
 604       //filter CE tracks
 605       Int_t id = FindMirror(track,seed);
 606       if (id>=0) counter++;
 607     }
 608     //
 609   }
 610   fNtracks=counter;
 611   if (counter<kMinTracks) return;
 612
 613   //FitDriftV();
 614   FitDriftV(0.3);
 615   if (!fFullCalib) return;
 616   static Bool_t init=kFALSE;
 617   if (!init){
 618     init = kTRUE;  // way around for PROOF - to be investigated
 619     MakeFitHistos();
 620   }
 621   //
 622   for (Int_t id=0; id<336; id++){
 623     //
 624     if (!fTracksEsdParam.At(id)) continue;
 625     if (fClusterSatur[id]>0.3)   continue;  // tracks in saturation
 626     DumpLaser(id);
 627     if ( AcceptLaser(id) && fFitZ[id]<0.5){
 628       RefitLaserJW(id);
 629       MakeDistHisto(id);
 630     }
 631   }
 632
 633 }
 634
 635 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
 636   //
 637   //
 638   //
 639   //  for (Int_t id=0; id<336; id++){
 640     //
 641     //
 642     if (!fTracksEsdParam.At(id)) return;
 643     if (!AcceptLaser(id)) return;
 644     Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
 645     //
 646     //
 647     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
 648     if (!hisdz) MakeFitHistos();
 649     hisdz = (TH1F*)fDeltaZ.At(id);
 650     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
 651     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
 652     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
 653     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
 654     TH1F * hisSignal = (TH1F*)fSignals.At(id);
 655     //
 656
 657     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 658     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 659     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 660     if (!param) return;
 661     if (!ltrp) return;
 662     if (!track) return;
 663     Double_t xyz[3];
 664     Double_t pxyz[3];
 665     Double_t lxyz[3];
 666     Double_t lpxyz[3];
 667     param->GetXYZ(xyz);
 668     param->GetPxPyPz(pxyz);
 669     ltrp->GetXYZ(lxyz);
 670     ltrp->GetPxPyPz(lpxyz);
 671     //
 672     Float_t dz   = param->GetZ()-ltrp->GetZ();
 673     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
 674     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
 675     if (hisdz) hisdz->Fill(dz);
 676     if (hisP3) hisP3->Fill(param->GetParameter()[3]);
 677     if (hisP4) hisP4->Fill(param->GetParameter()[4]);
 678     if (hisdphi) hisdphi->Fill(dphi);
 679     if (hisdphiP) hisdphiP->Fill(dphiP);
 680     if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
 681     // fill HisLaser
 682     xhis[0] = ltrp->GetId();
 683     xhis[1] = ltrp->GetSide();
 684     xhis[2] = ltrp->GetRod();
 685     xhis[3] = ltrp->GetBundle();
 686     xhis[4] = ltrp->GetBeam();
 687     xhis[5] = dphi;
 688     xhis[6] = fFitZ[id];
 689     xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
 690     xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
 691     xhis[9] = param->GetParameter()[4];
 692     xhis[10]= track->GetTPCNcls();
 693     xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
 694     // }
 695     fHisLaser->Fill(xhis);
 696 }
 697
 698 void AliTPCcalibLaser::FitDriftV(){
 699   //
 700   // Fit corrections to the drift velocity - linear approximation in the z and global y
 701   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
 702   //
 703   /*
 704     Formulas:
 705
 706     z  = s* (z0 - vd*(t-t0))
 707
 708     s  - side -1 and +1
 709     t0 - time 0
 710     vd - nominal drift velocity
 711     zs - miscalibrated position
 712
 713     zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
 714     vr  - relative change of the drift velocity
 715     dzt - vd*dt
 716     dr  = zz0-s*z
 717     ..
 718     ==>
 719     zs ~ z - s*vr*(z0-s*z)+s*dzt
 720     --------------------------------
 721     1. Correction function vr constant:
 722
 723
 724     dz = zs-z = -s*vr *(z0-s*z)+s*dzt
 725     dzs/dl = dz/dl +s*s*vr*dz/dl
 726     d(dz/dl) = vr*dz/dl
 727   */
 728
 729   const Float_t kZCut        = 240;    // remove the closest laser beam
 730   const Float_t kSaturCut    = 0.05;   // remove saturated lasers - cut on fraction of saturated
 731   const Float_t kDistCut     = 3;    // distance sigma cut
 732   const Float_t kDistCutAbs  = 0.25;
 733   const Float_t kMinClusters = 60;  // minimal amount of the clusters
 734   const Float_t kMinSignal   = 16;  // minimal mean height of the signal
 735   const Float_t kChi2Cut     = 0.1; // chi2 cut to accept drift fit
 736   static TLinearFitter fdriftA(3,"hyp2");
 737   static TLinearFitter fdriftC(3,"hyp2");
 738   static TLinearFitter fdriftAC(4,"hyp3");
 739   TVectorD fitA(3),fitC(3),fitAC(4);
 740
 741   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
 742   AliTPCParam  * tpcparam    = calib->GetParameters();
 743   //
 744   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
 745
 746   Float_t chi2A = 10;
 747   Float_t chi2C = 10;
 748   Float_t chi2AC = 10;
 749   Int_t npointsA=0;
 750   Int_t npointsC=0;
 751   Int_t npointsAC=0;
 752
 753
 754   for (Int_t iter=0; iter<3; iter++){
 755     fdriftA.ClearPoints();
 756     fdriftC.ClearPoints();
 757     fdriftAC.ClearPoints();
 758     //
 759     for (Int_t id=0; id<336; id++){
 760       if (!fTracksEsdParam.At(id)) continue;
 761       if (!AcceptLaser(id)) continue;
 762       if ( fClusterSatur[id]>kSaturCut)  continue;
 763       if ( fClusterCounter[id]<kMinClusters)  continue;
 764       AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 765       if (track->GetTPCsignal()<kMinSignal) continue;
 766       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 767       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 768
 769       Double_t xyz[3];
 770       Double_t pxyz[3];
 771       Double_t lxyz[3];
 772       Double_t lpxyz[3];
 773       param->GetXYZ(xyz);
 774       param->GetPxPyPz(pxyz);
 775       ltrp->GetXYZ(lxyz);
 776       ltrp->GetPxPyPz(lpxyz);
 777       if (TMath::Abs(lxyz[2])>kZCut) continue;
 778       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
 779       if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
 780       if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
 781       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
 782
 783       // drift distance
 784       Double_t zlength =  tpcparam->GetZLength(0);
 785       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
 786       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
 787       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
 788       if (ltrp->GetSide()==0){
 789         fdriftA.AddPoint(xxx,mdrift,1);
 790       }else{
 791         fdriftC.AddPoint(xxx,mdrift,1);
 792       }
 793       Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
 794       fdriftAC.AddPoint(xxx2,mdrift,1);
 795     }
 796     //
 797     if (fdriftA.GetNpoints()>10){
 798       //
 799       fdriftA.Eval();
 800       if (iter==0) fdriftA.EvalRobust(0.7);
 801       else fdriftA.EvalRobust(0.8);
 802       fdriftA.GetParameters(fitA);
 803       npointsA= fdriftA.GetNpoints();
 804       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
 805       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
 806         if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
 807         (*fFitAside)[0] = fitA[0];
 808         (*fFitAside)[1] = fitA[1];
 809         (*fFitAside)[2] = fitA[2];
 810         (*fFitAside)[3] = fdriftA.GetNpoints();
 811         (*fFitAside)[4] = chi2A;
 812       }
 813     }
 814     if (fdriftC.GetNpoints()>10){
 815       fdriftC.Eval();
 816       if (iter==0) fdriftC.EvalRobust(0.7);
 817       else fdriftC.EvalRobust(0.8);
 818       //
 819       fdriftC.GetParameters(fitC);
 820       npointsC= fdriftC.GetNpoints();
 821       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
 822       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
 823         if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
 824         (*fFitCside)[0] = fitC[0];
 825         (*fFitCside)[1] = fitC[1];
 826         (*fFitCside)[2] = fitC[2];
 827         (*fFitCside)[3] = fdriftC.GetNpoints();
 828         (*fFitCside)[4] = chi2C;
 829       }
 830     }
 831
 832     if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
 833       fdriftAC.Eval();
 834       if (iter==0) fdriftAC.EvalRobust(0.7);
 835       else fdriftAC.EvalRobust(0.8);
 836       //
 837       fdriftAC.GetParameters(fitAC);
 838       npointsAC= fdriftAC.GetNpoints();
 839       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
 840       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
 841     }
 842
 843     for (Int_t id=0; id<336; id++){
 844       if (!fTracksEsdParam.At(id)) continue;
 845       //
 846       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 847       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 848       Double_t xyz[3];
 849       Double_t pxyz[3];
 850       Double_t lxyz[3];
 851       Double_t lpxyz[3];
 852       param->GetXYZ(xyz);
 853       param->GetPxPyPz(pxyz);
 854       ltrp->GetXYZ(lxyz);
 855       ltrp->GetPxPyPz(lpxyz);
 856       Double_t zlength =  tpcparam->GetZLength(0);
 857       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
 858       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
 859
 860       Float_t fz =0;
 861       if (ltrp->GetSide()==0){
 862         fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
 863       }else{
 864         fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
 865       }
 866       if (npointsAC>10){
 867         fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
 868       }
 869       fFitZ[id]=mdrift-fz;
 870     }
 871     if (fStreamLevel>0){
 872       TTreeSRedirector *cstream = GetDebugStreamer();
 873       TTimeStamp tstamp(fTime);
 874       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 875       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 876       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 877       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 878       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 879       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 880       TVectorD vecGoofie(20);
 881       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 882       if (goofieArray)
 883         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 884           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 885           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 886         }
 887
 888       if (cstream){
 889         (*cstream)<<"driftv"<<
 890           "run="<<fRun<<              //  run number
 891           "event="<<fEvent<<          //  event number
 892           "time="<<fTime<<            //  time stamp of event
 893           "trigger="<<fTrigger<<      //  trigger
 894           "mag="<<fMagF<<             //  magnetic field
 895           // Environment values
 896           "press0="<<valuePressure0<<
 897           "press1="<<valuePressure1<<
 898           "pt0="<<ptrelative0<<
 899           "pt1="<<ptrelative1<<
 900           "temp0="<<temp0<<
 901           "temp1="<<temp1<<
 902           "vecGoofie.="<<&vecGoofie<<
 903           //
 904           //
 905           "iter="<<iter<<
 906           "driftA.="<<fFitAside<<
 907           "driftC.="<<fFitCside<<
 908           "driftAC.="<<fFitACside<<
 909           "chi2A="<<chi2A<<
 910           "chi2C="<<chi2C<<
 911           "chi2AC="<<chi2AC<<
 912           "nA="<<npointsA<<
 913           "nC="<<npointsC<<
 914           "nAC="<<npointsAC<<
 915           "\n";
 916       }
 917     }
 918   }
 919 }
 920 Bool_t  AliTPCcalibLaser::FitDriftV(Float_t minFraction){
 921   //
 922   // Fit corrections to the drift velocity - linear approximation in the z and global y
 923   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
 924   //
 925   // Source of outlyers :
 926   // 0. Track in the saturation - postpeak
 927   // 1. gating grid close the part of the signal for first bundle
 928
 929   // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
 930   // 1. Robust fit is used in the itteration number 0
 931   // only fraction of laser uted
 932   // 2. Only the tracks close to the fit used in the second itteration
 933   /*
 934     Formulas:
 935
 936     z  = s* (z0 - vd*(t-t0))
 937
 938     s  - side -1 and +1
 939     t0 - time 0
 940     vd - nominal drift velocity
 941     zs - miscalibrated position
 942
 943     zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
 944     vr  - relative change of the drift velocity
 945     dzt - vd*dt
 946     dr  = zz0-s*z
 947     ..
 948     ==>
 949     zs ~ z - s*vr*(z0-s*z)+s*dzt
 950     --------------------------------
 951     1. Correction function vr constant:
 952
 953
 954     dz = zs-z = -s*vr *(z0-s*z)+s*dzt
 955     dzs/dl = dz/dl +s*s*vr*dz/dl
 956     d(dz/dl) = vr*dz/dl
 957   */
 958   const Int_t   knLaser      = 336;     //n laser tracks
 959   const Float_t kFraction[2] = {0.70,0.95};    // robust fit fraction
 960
 961   const Float_t kSaturCut    = 0.05;    // remove saturated lasers - cut on fraction of saturated
 962   const Float_t kDistCut     = 3.;      // distance sigma cut - 3 sigma
 963   const Float_t kDistCutAbs  = 1.;      // absolute cut 1 cm
 964   const Float_t kMinClusters = 60.;      // minimal amount of the clusters
 965   const Float_t kMinSignal   = 10.;      // minimal mean height of the signal
 966   const Float_t kChi2Cut     = 1.0;     // chi2 cut to accept drift fit
 967   //
 968   static TLinearFitter fdriftA(3,"hyp2");
 969   static TLinearFitter fdriftC(3,"hyp2");
 970   static TLinearFitter fdriftAC(4,"hyp3");
 971   TVectorD fitA(3),fitC(3),fitAC(4);
 972
 973   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
 974   AliTPCParam  * tpcparam    = calib->GetParameters();
 975   //
 976   // reset old data
 977   //
 978   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
 979   if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
 980   if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
 981   for (Int_t i=0;i<5; i++){
 982     (*fFitCside)[i]=0;
 983     (*fFitAside)[i]=0;
 984   }
 985   //
 986   //
 987   Float_t chi2A = 10;
 988   Float_t chi2C = 10;
 989   Float_t chi2AC = 10;
 990   Int_t npointsA=0;
 991   Int_t npointsC=0;
 992   Int_t npointsAC=0;
 993   Int_t nbA[4]={0,0,0,0};
 994   Int_t nbC[4]={0,0,0,0};
 995   TVectorD vecZM(336);      // measured z potion of laser
 996   TVectorD vecA(336);       // accepted laser
 997   TVectorD vecZF(336);      // fitted position
 998   TVectorD vecDz(336);      // deltaZ
 999   TVectorD vecZS(336);      // surveyed position of laser
1000   // additional variable to cut
1001   TVectorD vecdEdx(336);    // dEdx
1002   TVectorD vecSy(336);      // shape y
1003   TVectorD vecSz(336);      // shape z
1004   //
1005   //
1006   for (Int_t id=0; id<336; id++){
1007     Int_t reject=0;
1008     AliTPCLaserTrack *ltrp =
1009       (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1010     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1011     AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
1012     vecZM(id)= (param==0) ? 0:param->GetZ();
1013     vecZS(id)= ltrp->GetZ();
1014     vecDz(id)= 0;
1015     vecA(id)=1;
1016     vecdEdx(id)=(seed)?seed->GetdEdx():0;
1017     vecSy(id) =(seed)?seed->CookShape(1):0;
1018     vecSz(id) =(seed)?seed->CookShape(2):0;
1019     //
1020     fFitZ[id]=0;
1021     if (!fTracksEsdParam.At(id)) reject|=1;
1022     if (!param) continue;
1023     if (!AcceptLaser(id))        reject|=2;
1024     if ( fClusterSatur[id]>kSaturCut)  reject|=4;
1025     if ( fClusterCounter[id]<kMinClusters)  reject|=8;
1026     vecA(id)=reject;
1027     if (reject>0) continue;
1028     if (ltrp->GetSide()==0){
1029       npointsA++;
1030       nbA[ltrp->GetBundle()]++;
1031     }
1032     if (ltrp->GetSide()==1){
1033       npointsC++;
1034       nbC[ltrp->GetBundle()]++;
1035     }
1036   }
1037   //
1038   // reject  "bad events"
1039   //
1040   Bool_t isOK=kTRUE;
1041   Int_t  naA=0;
1042   Int_t  naC=0;
1043   if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
1044     isOK=kFALSE;
1045   for (Int_t i=0;i<4;i++){
1046     //count accepted for all layers
1047     if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
1048     if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
1049   }
1050   if (naA<3 &&naC<3) isOK=kFALSE;
1051   if (isOK==kFALSE) return kFALSE;
1052
1053   //
1054   //
1055   //
1056   for (Int_t iter=0; iter<2; iter++){
1057     fdriftA.ClearPoints();
1058     fdriftC.ClearPoints();
1059     fdriftAC.ClearPoints();
1060     npointsA=0;  npointsC=0;  npointsAC=0;
1061     //
1062     for (Int_t id=0; id<336; id++){
1063       if (!fTracksEsdParam.At(id)) continue;
1064       if (!AcceptLaser(id)) continue;
1065       if ( fClusterSatur[id]>kSaturCut)  continue;
1066       if ( fClusterCounter[id]<kMinClusters)  continue;
1067       AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1068       if (track->GetTPCsignal()<kMinSignal) continue;
1069       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1070       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1071       Double_t xyz[3];
1072       Double_t pxyz[3];
1073       Double_t lxyz[3];
1074       Double_t lpxyz[3];
1075       param->GetXYZ(xyz);
1076       param->GetPxPyPz(pxyz);
1077       ltrp->GetXYZ(lxyz);
1078       ltrp->GetPxPyPz(lpxyz);
1079       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
1080       if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
1081       if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
1082       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
1083
1084 //       // drift distance
1085 //       Double_t zlength =  tpcparam->GetZLength(0);
1086 //       Double_t ldrift  = zlength-TMath::Abs(lxyz[2]);
1087 //       Double_t mdrift  = zlength-TMath::Abs(xyz[2]);
1088       //
1089       Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1090       Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1091       Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1092
1093
1094       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
1095       if (iter==0 &&ltrp->GetBundle()==0) continue;
1096       // skip bundle 0 - can be wrong in the 0 iteration
1097       if (ltrp->GetSide()==0){
1098         fdriftA.AddPoint(xxx,mdrift,1);
1099       }else{
1100         fdriftC.AddPoint(xxx,mdrift,1);
1101       }
1102       Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
1103       fdriftAC.AddPoint(xxx2,mdrift,1);
1104     }
1105     //
1106     if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
1107       //
1108       fdriftA.Eval();
1109       npointsA= fdriftA.GetNpoints();
1110       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
1111       fdriftA.EvalRobust(kFraction[iter]);
1112       fdriftA.GetParameters(fitA);
1113       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
1114         if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
1115         (*fFitAside)[0] = fitA[0];
1116         (*fFitAside)[1] = fitA[1];
1117         (*fFitAside)[2] = fitA[2];
1118         (*fFitAside)[3] = fdriftA.GetNpoints();
1119         (*fFitAside)[4] = chi2A;
1120       }
1121     }
1122     if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
1123       fdriftC.Eval();
1124       npointsC= fdriftC.GetNpoints();
1125       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
1126       fdriftC.EvalRobust(kFraction[iter]);
1127       fdriftC.GetParameters(fitC);
1128       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
1129         if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
1130         (*fFitCside)[0] = fitC[0];
1131         (*fFitCside)[1] = fitC[1];
1132         (*fFitCside)[2] = fitC[2];
1133         (*fFitCside)[3] = fdriftC.GetNpoints();
1134         (*fFitCside)[4] = chi2C;
1135       }
1136     }
1137
1138     if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
1139       fdriftAC.Eval();
1140       npointsAC= fdriftAC.GetNpoints();
1141       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
1142       fdriftAC.EvalRobust(kFraction[iter]);
1143       fdriftAC.GetParameters(fitAC);
1144       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
1145       (*fFitACside)[0] = fitAC[0];
1146       (*fFitACside)[1] = fitAC[1];
1147       (*fFitACside)[2] = fitAC[2];
1148       (*fFitACside)[3] = fdriftAC.GetNpoints();
1149       (*fFitACside)[4] = chi2AC;
1150     }
1151
1152     for (Int_t id=0; id<336; id++){
1153       if (!fTracksEsdParam.At(id)) continue;
1154       //
1155       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1156       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1157       Double_t xyz[3];
1158       Double_t pxyz[3];
1159       Double_t lxyz[3];
1160       Double_t lpxyz[3];
1161       param->GetXYZ(xyz);
1162       param->GetPxPyPz(pxyz);
1163       ltrp->GetXYZ(lxyz);
1164       ltrp->GetPxPyPz(lpxyz);
1165       //Double_t zlength =  tpcparam->GetZLength(0);
1166       //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
1167       //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
1168       Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
1169       Double_t ldrift  = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
1170       Double_t mdrift  = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
1171
1172
1173       Float_t fz =0;
1174       if (ltrp->GetSide()==0){
1175         fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
1176       }else{
1177         fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
1178       }
1179       if (npointsAC>10){
1180         //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
1181       }
1182       fFitZ[id]=mdrift-fz;
1183       vecZF[id]=fz;
1184       vecDz[id]=mdrift-fz;
1185     }
1186     if (fStreamLevel>0){
1187       TTreeSRedirector *cstream = GetDebugStreamer();
1188       TTimeStamp tstamp(fTime);
1189       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
1190       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
1191       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
1192       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
1193       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
1194       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
1195       TVectorD vecGoofie(20);
1196       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
1197       if (goofieArray)
1198         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
1199           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
1200           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
1201         }
1202
1203       if (cstream){
1204         (*cstream)<<"driftvN"<<
1205           "run="<<fRun<<              //  run number
1206           "event="<<fEvent<<          //  event number
1207           "time="<<fTime<<            //  time stamp of event
1208           "trigger="<<fTrigger<<      //  trigger
1209           "mag="<<fMagF<<             //  magnetic field
1210           // Environment values
1211           "press0="<<valuePressure0<<
1212           "press1="<<valuePressure1<<
1213           "pt0="<<ptrelative0<<
1214           "pt1="<<ptrelative1<<
1215           "temp0="<<temp0<<
1216           "temp1="<<temp1<<
1217           "vecGoofie.="<<&vecGoofie<<
1218           //
1219           //
1220           "vecZM.="<<&vecZM<<   // measured z position
1221           "vecZS.="<<&vecZS<<   // surveyed z position
1222           "vecZF.="<<&vecZF<<   // fitted   z position
1223           "vecDz.="<<&vecDz<<   // fitted   z position
1224           "vecA.="<<&vecA<<     // accept laser flag
1225           "vecdEdx.="<<&vecdEdx<<   // dEdx  - to cut on
1226           "vecSy.="<<&vecSy<<       // shape y - to cut on
1227           "vecSz.="<<&vecSz<<       // shape z - to cut on
1228           //
1229           "iter="<<iter<<
1230           "driftA.="<<fFitAside<<
1231           "driftC.="<<fFitCside<<
1232           "driftAC.="<<fFitACside<<
1233           "chi2A="<<chi2A<<
1234           "chi2C="<<chi2C<<
1235           "chi2AC="<<chi2AC<<
1236           "nA="<<npointsA<<
1237           "nC="<<npointsC<<
1238           "nAC="<<npointsAC<<
1239           "\n";
1240         /*
1241           //
1242           variables to check in debug mode:
1243           //
1244           chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
1245           chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
1246           chainDriftN->SetAlias("driftF","vecZF.fElements");
1247           chainDriftN->SetAlias("deltaZ","driftF-driftM");  //deltaZ
1248           TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
1249           TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
1250
1251          */
1252       }
1253     }
1254   }
1255   return kTRUE;
1256 }
1257
1258 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
1259   //
1260   // get distance between mirror and laser track
1261   //
1262   //
1263   Double_t xyz[3];
1264   Double_t lxyz[3];
1265   param->GetXYZ(xyz);
1266   ltrp->GetXYZ(lxyz);
1267   //
1268   //
1269   Double_t dist = 0;
1270   //radial distance
1271   dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
1272   //
1273   // z distance
1274   // apply drift correction if already exist
1275   //
1276   Float_t dz = 0;
1277   if (ltrp->GetSide()==0){
1278     if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
1279   }else{
1280     if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
1281   }
1282   if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
1283   dist+=TMath::Abs(dz);
1284   //
1285   // phi dist - divergence on 50 cm
1286   //
1287   dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
1288   return dist;
1289 }
1290
1291
1292 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
1293   //
1294   //
1295   //
1296   /*
1297   TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
1298   TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
1299   TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
1300   TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
1301
1302   TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
1303   */
1304   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
1305   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1306   //AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1307   Double_t xyz[3];
1308   Double_t lxyz[3];
1309   param->GetXYZ(xyz);
1310   ltrp->GetXYZ(lxyz);
1311   if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
1312   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE;    // cutZ -P1
1313   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE;  // cut -P2
1314   if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE;   // cut Tl -P3
1315   //
1316   //
1317
1318   return kTRUE;
1319 }
1320
1321 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
1322   //
1323   // Find corresponding mirror
1324   // add the corresponding tracks
1325
1326   if (!track->GetOuterParam()) return -1;
1327
1328   Float_t kRadius0  = 252;
1329   Float_t kRadius   = 254.2;
1330   Int_t countercl=0;
1331   Float_t counterSatur=0;
1332   Int_t csideA =0;
1333   Int_t csideC =0;
1334   for (Int_t irow=158;irow>-1;--irow) {
1335     AliTPCclusterMI *c=seed->GetClusterPointer(irow);
1336     if (!c) continue;
1337     Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
1338     Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1339     if (pedgeY<3) continue;
1340     if (pedgeX<3) continue;
1341     countercl++;
1342     if (c->GetDetector()%36<18) csideA++;
1343     if (c->GetDetector()%36>=18) csideC++;
1344     if (c->GetMax()>900) counterSatur++;
1345   }
1346   counterSatur/=(countercl+1);
1347   //
1348   //
1349   //
1350   if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0;  // cross laser track can not happen
1351
1352   Int_t side= 0;
1353   if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
1354
1355
1356   AliExternalTrackParam param(*(track->GetOuterParam()));
1357   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
1358   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
1359   AliTPCLaserTrack ltr;
1360   AliTPCLaserTrack *ltrp=0x0;
1361   //  AliTPCLaserTrack *ltrpjw=0x0;
1362   //
1363   Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param,side);
1364  // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
1365   //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
1366
1367   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
1368     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1369   else
1370     ltrp=&ltr;
1371
1372   if (id<0) return -1;
1373   if (ltrp->GetSide()!=side) return -1;
1374   fCounter[id]++;
1375   //
1376   //
1377   //
1378   //
1379   if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
1380   //
1381   //
1382   Float_t radius=TMath::Abs(ltrp->GetX());
1383   param.Rotate(ltrp->GetAlpha());
1384   AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kFALSE);
1385   //
1386   if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
1387   Bool_t accept=kTRUE;
1388   //
1389   // choose closer track
1390   //
1391   AliExternalTrackParam * param0  = (AliExternalTrackParam *)fTracksEsdParam.At(id);
1392   if (param0){
1393     Float_t dist0=GetDistance(param0,ltrp);
1394     Float_t dist1=GetDistance(&param,ltrp);
1395     if (dist0<dist1)    accept=kFALSE;
1396   }
1397
1398   if (accept){
1399     fClusterCounter[id]=countercl;
1400     fTracksEsdParam.AddAt(param.Clone(),id);
1401     fTracksEsd.AddAt(track,id);
1402     fTracksTPC.AddAt(seed,id);
1403   }
1404   return id;
1405 }
1406
1407
1408
1409 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1410   //
1411   //  Dump Laser info to the tree
1412   //
1413   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1414   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1415   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1416   //
1417   // Fast laser ID
1418   //
1419   Double_t xyz[3];
1420   Double_t pxyz[3];
1421   Double_t lxyz[3];
1422   Double_t lpxyz[3];
1423   param->GetXYZ(xyz);
1424   param->GetPxPyPz(pxyz);
1425   ltrp->GetXYZ(lxyz);
1426   ltrp->GetPxPyPz(lpxyz);
1427   Float_t dist3D   = GetDistance(param,ltrp);
1428   Float_t dist0    = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1429   Float_t distphi  = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1430
1431
1432   if (fStreamLevel>0){
1433     TTreeSRedirector *cstream = GetDebugStreamer();
1434     Int_t time = fESD->GetTimeStamp();
1435     Bool_t accept = AcceptLaser(id);
1436     if (cstream){
1437       (*cstream)<<"Track"<<
1438         //
1439         "run="<<fRun<<              //  run number
1440         "event="<<fEvent<<          //  event number
1441         "time="<<fTime<<            //  time stamp of event
1442         "trigger="<<fTrigger<<      //  trigger
1443         "mag="<<fMagF<<             //  magnetic field
1444
1445         "id="<<id<<
1446         "accept="<<accept<<
1447         "driftA.="<<fFitAside<<
1448         "driftC.="<<fFitCside<<
1449         "time="<<time<<
1450         "dist3D="<<dist3D<<
1451         "dist0="<<dist0<<
1452         "distphi="<<distphi<<
1453         //
1454         //
1455         "counter="<<fCounter[id]<<
1456         "clcounter="<<fClusterCounter[id]<<
1457         "satur="<<fClusterSatur[id]<<
1458         "fz="<<fFitZ[id]<<
1459         //
1460         "LTr.="<<ltrp<<
1461         "Esd.="<<track<<
1462         "Tr.="<<param<<
1463         "x0="<<xyz[0]<<
1464         "x1="<<xyz[1]<<
1465         "x2="<<xyz[2]<<
1466         "px0="<<pxyz[0]<<
1467         "px1="<<pxyz[1]<<
1468         "px2="<<pxyz[2]<<
1469         //
1470         "lx0="<<lxyz[0]<<
1471         "lx1="<<lxyz[1]<<
1472         "lx2="<<lxyz[2]<<
1473         "lpx0="<<lpxyz[0]<<
1474         "lpx1="<<lpxyz[1]<<
1475         "lpx2="<<lpxyz[2]<<
1476         "\n";
1477     }
1478   }
1479 }
1480
1481 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1482   //
1483   // Refit the track with different tracklet models:
1484   // 1. Per ROC using the kalman filter, different edge cuts
1485   // 2. Per ROC linear in y and z
1486   // 3. Per ROC quadratic in y and z
1487   // 4. Per track offset for each sector, linear for each sector, common quadratic
1488   // store x, y, z information for all models and the cluster to calculate the residuals
1489   //
1490
1491   // The clusters which do not fulfill given  criteria are skipped
1492   //
1493   // Cluters removed from fit
1494   //  1. Extended shape     > kShapeCut
1495   //  2. In saturation  Max > kMax
1496   //  3. Distance to edge   < cutEdge
1497   //
1498   // Clusters not used for the calibration:
1499   //
1500   //  1. Extended shape     > kShapeCut
1501   //  2. In saturation  Max > kMax
1502
1503
1504   AliTPCseed *track      = (AliTPCseed*)fTracksTPC.At(id);
1505   AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1506   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1507
1508   AliTPCclusterMI dummyCl;
1509
1510   //two tracklets
1511   Int_t kMaxTracklets=2;
1512   Float_t kShapeCut = 1.3;
1513   Float_t kRatioCut = 2.;
1514
1515   Float_t kMax      = 900.;
1516
1517
1518   //=============================================//
1519   // Linear Fitters for the Different Approaches //
1520   //=============================================//
1521   //linear fit model in y and z; inner - outer sector, combined with offset
1522   static TLinearFitter fy1I(2,"hyp1");
1523   static TLinearFitter fy1O(2,"hyp1");
1524   static TLinearFitter fz1I(2,"hyp1");
1525   static TLinearFitter fz1O(2,"hyp1");
1526   static TLinearFitter fy1IO(3,"hyp2");
1527   static TLinearFitter fz1IO(3,"hyp2");
1528   //quadratic fit model in y and z; inner - sector
1529   static TLinearFitter fy2I(3,"hyp2");
1530   static TLinearFitter fy2O(3,"hyp2");
1531   static TLinearFitter fz2I(3,"hyp2");
1532   static TLinearFitter fz2O(3,"hyp2");
1533   //common quadratic fit for IROC and OROC in y and z
1534   static TLinearFitter fy4(5,"hyp4");
1535   static TLinearFitter fz4(5,"hyp4");
1536
1537
1538   //set standard cuts
1539   if ( fNcuts==0 ){
1540       fNcuts=1;
1541       fEdgeXcuts[0]=4;
1542       fEdgeYcuts[0]=3;
1543       fNClCuts[0]=20;
1544   }
1545   //=============================//
1546   // Loop over all Tracklet Cuts //
1547   //=============================//
1548   for (Int_t icut=0; icut<fNcuts; icut++){
1549     Float_t xinMin = 2500, xinMax=-90;
1550     Float_t xoutMin=2500, xoutMax=-90;
1551     Float_t msigmaYIn=0;
1552     Float_t msigmaYOut=0;
1553     Float_t mqratioIn=0;
1554     Float_t mqratioOut=0;
1555
1556
1557       AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1558       //cut parameters
1559       Double_t edgeCutX = fEdgeXcuts[icut];
1560       Double_t edgeCutY = fEdgeYcuts[icut];
1561       Int_t    nclCut   = (Int_t)fNClCuts[icut];
1562       //=========================//
1563       // Parameters to calculate //
1564       //=========================//
1565       //fit parameter inner, outer tracklet and combined fit
1566       TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1567       TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1568       //
1569       TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1570       TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1571       TVectorD vecy4res(5),vecz4res(5);
1572       TVectorD vecy1resIO(3),vecz1resIO(3);
1573       // cluster and track positions for each row - used for residuals
1574       TVectorD vecgX(159);        // global X
1575       TVectorD vecgY(159);        // global Y
1576       TVectorD vecgZ(159);        // global Z
1577
1578       TVectorD vecX(159);        // x is the same for all (row center)
1579       TVectorD vecYkalman(159);  // y from kalman fit
1580       TVectorD vecZkalman(159);  // z from kalman fit
1581       TVectorD vecY1(159);       // y from pol1 fit per ROC
1582       TVectorD vecZ1(159);       // z from pol1 fit per ROC
1583       TVectorD vecY1IO(159);     // y from pol1 fit per ROC
1584       TVectorD vecZ1IO(159);     // z from pol1 fit per ROC
1585       TVectorD vecY2(159);       // y from pol2 fit per ROC
1586       TVectorD vecZ2(159);       // z from pol2 fit per ROC
1587       TVectorD vecY4(159);       // y from sector fit
1588       TVectorD vecZ4(159);       // z from sector fit
1589       TVectorD vecClY(159);      // y cluster position
1590       TVectorD vecClZ(159);      // z cluster position
1591       TVectorD vecSec(159);      // sector for each row
1592       TVectorD  isReject(159);    // flag - cluster to be rejected
1593       //chi2 of fits
1594       Double_t chi2I1z=0;       // chi2 of pol1 fit in z (inner)
1595       Double_t chi2I1y=0;       // chi2 of pol1 fit in y (inner)
1596       Double_t chi2O1z=0;       // chi2 of pol1 fit in z (outer)
1597       Double_t chi2O1y=0;       // chi2 of pol1 fit in y (outer)
1598       Double_t chi2IO1z=0;       // chi2 of pol1 fit in z (outer)
1599       Double_t chi2IO1y=0;       // chi2 of pol1 fit in y (outer)
1600       Double_t chi2I2z=0;       // chi2 of pol2 fit in z (inner)
1601       Double_t chi2I2y=0;       // chi2 of pol2 fit in y (inner)
1602       Double_t chi2O2z=0;       // chi2 of pol2 fit in z (outer)
1603       Double_t chi2O2y=0;       // chi2 of pol2 fit in y (outer)
1604       Double_t chi2IOz=0;       // chi2 of hyp4 fit in z (inner+outer)
1605       Double_t chi2IOy=0;       // chi2 of hyp4 fit in y (inner+outer)
1606       //more
1607       Int_t innerSector = -1;    // number of inner sector
1608       Int_t outerSector = -1;    // number of outer sector
1609       Int_t nclI=0;              // number of clusters (inner)
1610       Int_t nclO=0;              // number of clusters (outer)
1611       //=================================================//
1612       // Perform the Kalman Fit using the Tracklet Class //
1613       //=================================================//
1614       AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1615       TObjArray tracklets=
1616           AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1617                                           kFALSE,nclCut,kMaxTracklets);
1618       // tracklet pointers
1619       AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1620       AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1621       AliTPCTracklet *tr=0x0;
1622       AliTPCTracklet dummy;
1623       //continue if we didn't find a tracklet
1624       if ( !trInner && !trOuter ) continue;
1625       //================================================================================//
1626       // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1627       //================================================================================//
1628       if ( trInner && trOuter ){
1629           if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1630           if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1631               tr = trInner;
1632               trInner=trOuter;
1633               trOuter=tr;
1634               AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1635           }
1636       } else {
1637           if ( trInner ){
1638               if ( !trInner->GetInner() ) continue;
1639               trOuter=&dummy;
1640               if ( trInner->GetSector()>35 ){
1641                   trOuter=trInner;
1642                   trInner=&dummy;
1643               }
1644           } else { //trOuter
1645               if ( !trOuter->GetInner() ) continue;
1646               trInner=&dummy;
1647               if ( trOuter->GetSector()<36 ){
1648                   trInner=trOuter;
1649                   trOuter=&dummy;
1650               }
1651           }
1652       }
1653       innerSector = trInner->GetSector();
1654       if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1655       outerSector = trOuter->GetSector();
1656       if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1657
1658       // array of clusters
1659       TClonesArray arrCl("AliTPCclusterMI",159);
1660       arrCl.ExpandCreateFast(159);
1661       //=======================================//
1662       // fill fitters with cluster information //
1663       //=======================================//
1664       AliDebug(3,"Filing Cluster Information");
1665
1666       //
1667       //
1668
1669       for (Int_t irow=158;irow>-1;--irow) {
1670           AliTPCclusterMI *c=track->GetClusterPointer(irow);
1671           AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1672           cl=dummyCl;
1673           vecX[irow]   = 0;
1674           vecClY[irow] = 0;
1675           vecClZ[irow] = 0;
1676           Float_t meanY=0, sumY=0;
1677           for (Int_t drow=-1;drow<=1;drow++) {
1678             if (irow+drow<0) continue;
1679             if (irow+drow>158) continue;
1680             AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1681             if (!ccurrent) continue;
1682             Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1683             if ( roc!=innerSector && roc!=outerSector ) continue;
1684             if (ccurrent->GetX()<10) continue;
1685             if (ccurrent->GetY()==0) continue;
1686             meanY+=ccurrent->GetY();
1687             sumY++;
1688           }
1689           if (sumY>0)  meanY/=sumY;
1690
1691           //
1692           vecSec[irow]=-1;
1693           if (!c) continue;
1694           Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1695           Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1696
1697           //
1698           Int_t roc = static_cast<Int_t>(c->GetDetector());
1699           if ( roc!=innerSector && roc!=outerSector ) continue;
1700           vecSec[irow]=roc;
1701           //store clusters in clones array
1702           cl=*c;
1703           //
1704           if (c->GetMax()<4)   continue;  // noise cluster?
1705           if (TMath::Abs(c->GetY())<0.0001)   continue;  // noise cluster?
1706           //cluster position
1707           vecX[irow]   = c->GetX();
1708           vecClY[irow] = c->GetY();
1709           vecClZ[irow] = c->GetZ();
1710           //
1711 //        Float_t gxyz[3];
1712 //        c->GetGlobalXYZ(gxyz);
1713 //        vecgX[irow]   = gxyz[0];
1714 //        vecgY[irow]   = gxyz[1];
1715 //        vecgZ[irow]   = gxyz[2];
1716           //
1717           Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1718           Double_t y = vecClY[irow];
1719           Double_t z = vecClZ[irow];
1720           //
1721           Double_t x2[2]={x,x*x};   //linear and parabolic parameters
1722           Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1723           Double_t xIO[2]={0,x};    //common linear + offset IROC-OROC
1724           if ( roc == innerSector ) {
1725               x4[0]=1; //offset inner - outer sector
1726               x4[1]=x; //slope parameter inner sector
1727               xIO[0]=1;
1728           } else {
1729               x4[2]=x; //slope parameter outer sector
1730           }
1731           x4[3]=x*x;   //common parabolic shape
1732           if (pedgeX < fEdgeXcuts[icut]) continue;
1733           if (pedgeY < fEdgeYcuts[icut]) continue;
1734           if (c->GetMax()>900) continue;  // cluster in saturation
1735           //
1736           if ( roc==innerSector ){
1737               fy1I.AddPoint(x2,y);
1738               fz1I.AddPoint(x2,z);
1739               fy2I.AddPoint(x2,y);
1740               fz2I.AddPoint(x2,z);
1741               ++nclI;
1742               if (c->GetX()<xinMin) xinMin=c->GetX();
1743               if (c->GetX()>xinMax) xinMax=c->GetX();
1744
1745               msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1746               mqratioIn +=c->GetMax()/c->GetQ();
1747
1748           }
1749           if ( roc==outerSector ){
1750               fy1O.AddPoint(x2,y);
1751               fz1O.AddPoint(x2,z);
1752               fy2O.AddPoint(x2,y);
1753               fz2O.AddPoint(x2,z);
1754               ++nclO;
1755               if (c->GetX()<xoutMin) xoutMin=c->GetX();
1756               if (c->GetX()>xoutMax) xoutMax=c->GetX();
1757               msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1758               mqratioOut +=c->GetMax()/c->GetQ();
1759
1760           }
1761           fy4.AddPoint(x4,y);
1762           fz4.AddPoint(x4,z);
1763           fy1IO.AddPoint(xIO,y);
1764           fz1IO.AddPoint(xIO,z);
1765       }
1766       if (nclI>0)  {
1767         msigmaYIn/=nclI;
1768         mqratioIn/=nclI;
1769       }
1770       if (nclO>0)  {
1771         msigmaYOut/=nclO;
1772         mqratioOut/=nclO;
1773       }
1774       //======================================//
1775       // Evaluate and retrieve fit parameters //
1776       //======================================//
1777       AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1778       //inner sector
1779       if (  (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1780           fy1I.Eval();
1781           fz1I.Eval();
1782           fy2I.Eval();
1783           fz2I.Eval();
1784           fy1I.GetParameters(vecy1resInner);
1785           fz1I.GetParameters(vecz1resInner);
1786           fy2I.GetParameters(vecy2resInner);
1787           fz2I.GetParameters(vecz2resInner);
1788           chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1789           chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1790           chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1791           chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1792       }
1793       //outer sector
1794       if (  (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1795           fy1O.Eval();
1796           fz1O.Eval();
1797           fy2O.Eval();
1798           fz2O.Eval();
1799           fy1O.GetParameters(vecy1resOuter);
1800           fz1O.GetParameters(vecz1resOuter);
1801           fy2O.GetParameters(vecy2resOuter);
1802           fz2O.GetParameters(vecz2resOuter);
1803           chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1804           chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1805           chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1806           chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1807       }
1808       //combined hyp4 fit
1809       if ( innerSector>0 && outerSector>0 ){
1810           if (fy4.GetNpoints()>0) {
1811               fy4.Eval();
1812               fy4.GetParameters(vecy4res);
1813               chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1814           }
1815           if (fz4.GetNpoints()>0) {
1816               fz4.Eval();
1817               fz4.GetParameters(vecz4res);
1818               chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1819           }
1820           if (fy1IO.GetNpoints()>0) {
1821             fy1IO.Eval();
1822             fy1IO.GetParameters(vecy1resIO);
1823             chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1824           }
1825           if (fz1IO.GetNpoints()>0) {
1826             fz1IO.Eval();
1827             fz1IO.GetParameters(vecz1resIO);
1828             chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1829           }
1830       }
1831       //clear points
1832       fy4.ClearPoints();  fz4.ClearPoints();
1833       fy1I.ClearPoints(); fy1O.ClearPoints();
1834       fz1I.ClearPoints(); fz1O.ClearPoints();
1835       fy2I.ClearPoints(); fy2O.ClearPoints();
1836       fz2I.ClearPoints(); fz2O.ClearPoints();
1837       fy1IO.ClearPoints(); fz1IO.ClearPoints();
1838       //==============================//
1839       // calculate tracklet positions //
1840       //==============================//
1841       AliDebug(4,"Calculate tracklet positions");
1842       for (Int_t irow=158;irow>-1;--irow) {
1843         isReject[irow]=0;
1844          AliTPCclusterMI *c=track->GetClusterPointer(irow);
1845          if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1846            isReject[irow]+=1;
1847          }
1848
1849          if (!c) { //no cluster
1850            isReject[irow]+=2;
1851          }else{
1852            if (c->GetMax()>kMax){   //saturation
1853              isReject[irow]+=4;
1854            }
1855            if ( vecSec[irow] == outerSector ) {   //extended shape
1856              if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut)                isReject[irow]+=8;
1857              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut)    isReject[irow]+=16;
1858            }else{
1859              if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut)                isReject[irow]+=8;
1860              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut)    isReject[irow]+=16;
1861            }
1862          }
1863
1864
1865
1866           if ( vecSec[irow]==-1 ) continue;  //no cluster info
1867           if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1868           tr=&dummy;
1869           Double_t x=vecX[irow];
1870           Double_t xref=x-133.4;
1871           //
1872           Double_t yoffInner=0;
1873           Double_t zoffInner=0;
1874           Double_t yoffInner1=0;
1875           Double_t zoffInner1=0;
1876           Double_t yslopeInner=0;
1877           Double_t yslopeOuter=0;
1878           Double_t zslopeInner=0;
1879           Double_t zslopeOuter=0;
1880           //positions of hyperplane fits
1881           if ( vecSec[irow] == outerSector ) {
1882               tr=trOuter;
1883               vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1884               vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1885               vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1886               vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1887               yslopeOuter=vecy4res[3];
1888               zslopeOuter=vecz4res[3];
1889           } else {
1890               tr=trInner;
1891               vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1892               vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1893               vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1894               vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1895               yoffInner=vecy4res[1];
1896               zoffInner=vecz4res[1];
1897               yoffInner1=vecy1resIO[1];
1898               zoffInner1=vecz1resIO[1];
1899               yslopeInner=vecy4res[2];
1900               zslopeInner=vecz4res[2];
1901           }
1902           vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1903           vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1904           vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1905           vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1906           //positions of kalman fits
1907           Double_t gxyz[3],xyz[3];
1908           AliExternalTrackParam *param = 0x0;
1909           //
1910           param=tr->GetInner();
1911           if (param){
1912               param->GetXYZ(gxyz);
1913               Float_t bz = AliTracker::GetBz(gxyz);
1914               param->GetYAt(x, bz, xyz[1]);
1915               param->GetZAt(x, bz, xyz[2]);
1916               vecYkalman[irow]=xyz[1];
1917               vecZkalman[irow]=xyz[2];
1918           }
1919           //
1920           //
1921           //
1922
1923       }
1924       //=====================================================================//
1925       // write results from the different tracklet fits with debug streamers //
1926       //=====================================================================//
1927       if (fStreamLevel>4){
1928           TTreeSRedirector *cstream = GetDebugStreamer();
1929           if (cstream){
1930             Float_t dedx = track->GetdEdx();
1931             (*cstream)<<"FitModels"<<
1932               "run="<<fRun<<              //  run number
1933               "event="<<fEvent<<          //  event number
1934               "time="<<fTime<<            //  time stamp of event
1935               "trigger="<<fTrigger<<      //  trigger
1936               "mag="<<fMagF<<             //  magnetic field
1937               //
1938               "cutNr="      << icut <<
1939               "edgeCutX="   << edgeCutX <<
1940               "edgeCutY="   << edgeCutY <<
1941               "nclCut="     << nclCut <<
1942               "innerSector="<< innerSector <<
1943               "outerSector="<< outerSector <<
1944               "dEdx="       << dedx <<
1945               "LTr.="       << ltrp <<
1946               "Tr.="        << extparam <<
1947               "yPol1In.="   << &vecy1resInner <<
1948               "zPol1In.="   << &vecz1resInner <<
1949               "yPol1InOut.="<< &vecy1resIO <<
1950               "zPol1InOut.="<< &vecz1resIO <<
1951               "yPol2In.="   << &vecy2resInner <<
1952               "zPol2In.="   << &vecz2resInner <<
1953               "yPol1Out.="  << &vecy1resOuter <<
1954               "zPol1Out.="  << &vecz1resOuter <<
1955               "yPol2Out.="  << &vecy2resOuter <<
1956               "zPol2Out.="  << &vecz2resOuter <<
1957               "yInOut.="    << &vecy4res <<
1958               "zInOut.="    << &vecz4res <<
1959               "chi2y1In="   << chi2I1y <<
1960               "chi2z1In="   << chi2I1z <<
1961               "chi2y1InOut="<< chi2IO1y <<
1962               "chi2z1InOut="<< chi2IO1z <<
1963               "chi2y1Out="  << chi2O1y <<
1964               "chi2z1Out="  << chi2O1z <<
1965               "chi2y2In="   << chi2I2y <<
1966               "chi2z2In="   << chi2I2z <<
1967               "chi2y2Out="  << chi2O2y <<
1968               "chi2z2Out="  << chi2O2z <<
1969               "chi2yInOut=" << chi2IOy <<
1970               "chi2zInOut=" << chi2IOz <<
1971               "trletIn.="   << trInner <<
1972               "trletOut.="  << trOuter <<
1973               "nclI="       << nclI <<
1974               "nclO="       << nclO <<
1975               "xinMin="     << xinMin<<
1976               "xinMax="     << xinMax<<
1977               "xoutMin="    << xoutMin<<
1978               "xoutMax="     << xoutMax<<
1979               "msigmaYIn=" <<msigmaYIn<<
1980               "msigmaYOut=" <<msigmaYOut<<
1981               "mqratioIn=" <<mqratioIn<<
1982               "mqratioOut=" <<  mqratioOut      <<
1983               "\n";
1984           }
1985       }
1986
1987       // wirte residuals information
1988       if (fStreamLevel>5){
1989           TTreeSRedirector *cstream = GetDebugStreamer();
1990           if (cstream){
1991             Float_t dedx = track->GetdEdx();
1992             (*cstream)<<"Residuals"<<
1993               "run="<<fRun<<              //  run number
1994               "event="<<fEvent<<          //  event number
1995               "time="<<fTime<<            //  time stamp of event
1996               "trigger="<<fTrigger<<      //  trigger
1997               "mag="<<fMagF<<             //  magnetic field
1998               //
1999               "cutNr="      << icut <<
2000               "edgeCutX="   << edgeCutX <<
2001               "edgeCutY="   << edgeCutY <<
2002               "nclCut="     << nclCut   <<
2003               "LTr.="       << ltrp <<
2004               "Tr.="        << extparam<<
2005               "dEdx="       << dedx <<
2006               "Cl.="        << &arrCl <<
2007               "vX.="        << &vecgX<<   // global x
2008               "vY.="        << &vecgY<<   // global y
2009               "vZ.="        << &vecgZ<<   // global z
2010               "TrX.="       << &vecX <<
2011               "TrYpol1.="   << &vecY1 <<
2012               "TrZpol1.="   << &vecZ1 <<
2013               "TrYpol2.="   << &vecY2 <<
2014               "TrZpol2.="   << &vecZ2 <<
2015               "TrYpol1InOut.="<< &vecY1IO <<
2016               "TrZpol1InOut.="<< &vecZ1IO <<
2017               "TrYInOut.="  << &vecY4 <<
2018               "TrZInOut.="  << &vecZ4 <<
2019               "ClY.="       << &vecClY <<
2020               "ClZ.="       << &vecClZ <<
2021               "isReject.="  << &isReject<<
2022               "sec.="       << &vecSec <<
2023               "nclI="       << nclI <<
2024               "nclO="       << nclO <<
2025               "xinMin="     << xinMin<<
2026               "xinMax="     << xinMax<<
2027               "xoutMin="    << xoutMin<<
2028               "xoutMax="     << xoutMax<<
2029               "msigmaYIn=" <<msigmaYIn<<
2030               "msigmaYOut=" <<msigmaYOut<<
2031               "mqratioIn=" <<mqratioIn<<
2032               "mqratioOut=" <<  mqratioOut      <<
2033               "yInOut.="    << &vecy4res <<
2034               "zInOut.="    << &vecz4res <<
2035               //chi2s
2036               "chi2y1In="   << chi2I1y <<    //
2037               "chi2z1In="   << chi2I1z <<
2038               "chi2y1Out="  << chi2O1y <<
2039               "chi2z1Out="  << chi2O1z <<
2040               "chi2y1InOut="<< chi2IO1y <<
2041               "chi2z1InOut="<< chi2IO1z <<
2042               "chi2y2In="   << chi2I2y <<
2043               "chi2z2In="   << chi2I2z <<
2044               "chi2y2Out="  << chi2O2y <<
2045               "chi2z2Out="  << chi2O2z <<
2046               "chi2yInOut=" << chi2IOy <<
2047               "chi2zInOut=" << chi2IOz <<
2048               // fit parameters
2049               "yPol1In.="   << &vecy1resInner <<
2050               "zPol1In.="   << &vecz1resInner <<
2051               "yPol2In.="   << &vecy2resInner <<
2052               "zPol2In.="   << &vecz2resInner <<
2053               "yPol1Out.="  << &vecy1resOuter <<
2054               "zPol1Out.="  << &vecz1resOuter <<
2055               "yPol1InOut.="<< &vecy1resIO <<
2056               "zPol1InOut.="<< &vecz1resIO <<
2057               "yPol2Out.="  << &vecy2resOuter <<
2058               "zPol2Out.="  << &vecz2resOuter <<
2059
2060               "\n";
2061
2062           }
2063       }
2064       //==========================//
2065       // Fill Residual Histograms //
2066       //==========================//
2067       if (!fHisNclIn) MakeFitHistos();
2068
2069       TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2070       TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2071       TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2072       TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2073       TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2074       TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
2075       //      TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2076       //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2077       //
2078       for (Int_t irow=158;irow>-1;--irow) {
2079         if (vecSec[irow]==-1)continue; // no cluster info
2080         if (isReject[irow]>0.5) continue;  //
2081         Double_t ycl  = vecClY[irow];
2082         Double_t yfit = vecY1[irow];
2083         Double_t yfit2 = vecY2[irow];
2084         Double_t x    = vecX[irow];
2085         Double_t yabsbeam = -1000;
2086         if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2087           yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
2088         else if(innerSector==fBeamSectorInner[id])
2089           yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
2090
2091         //      Double_t yfit3 = vecY2[irow];
2092         Double_t zcl  = vecClZ[irow];
2093         Double_t zfit = vecZ1[irow];
2094         Double_t zfit2 = vecZ2[irow];
2095         //Double_t zfit3 = vecZ2[irow];
2096
2097         // dz abs
2098         // The expressions for zcorrected has been obtained by
2099         // inverting the fits in the FitDriftV() method (ignoring the
2100         // global y dependence for now):
2101         // A side:
2102         // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
2103         // =>
2104         // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
2105         //
2106         // C side:
2107         // 250 + zmeasured  = [0] + [1]*(250+zreal) + .... (yglobal)
2108         // =>
2109         // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
2110
2111         Double_t dzabs = -1000;
2112         Double_t zcorrected = -1000;
2113         if (ltrp->GetSide()==0){
2114           if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
2115             // ignore global y dependence for now
2116             zcorrected = 0;
2117             if(!fUseFixedDriftV)
2118               zcorrected = (zcl + (*fFitAside)[0] -
2119                             (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
2120             else
2121               zcorrected = (zcl + fFixedFitAside0 -
2122                             (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
2123             //      zcorrected = zcl;
2124             if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2125               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2126             else if(innerSector==fBeamSectorInner[id])
2127               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2128           }
2129         } else {
2130           if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
2131
2132             if(!fUseFixedDriftV)
2133               zcorrected = (zcl - (*fFitCside)[0] +
2134                             (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
2135             else
2136               zcorrected = (zcl - fFixedFitCside0 +
2137                             (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
2138
2139             //      zcorrected = zcl;
2140             if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
2141               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
2142             else if(innerSector==fBeamSectorInner[id])
2143               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
2144           }
2145         }
2146
2147         if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
2148           if (profy){
2149               profy->Fill(irow,ycl-yfit);
2150               profy2->Fill(irow,ycl-yfit2);
2151               if(yabsbeam<-100) {
2152                 fHisYAbsErrors->Fill(id);
2153                 //              profyabs->Fill(irow,-0.99);
2154               } else
2155                 profyabs->Fill(irow,ycl-yabsbeam);
2156
2157               //              profy3->Fill(irow,ycl-yfit3);
2158           }
2159           if (profz) {
2160               profz->Fill(irow,zcl-zfit);
2161               profz2->Fill(irow,zcl-zfit2);
2162               //profz3->Fill(irow,zcl-zfit3);
2163               if(dzabs<-100) {
2164
2165                 fHisZAbsErrors->Fill(id);
2166               }else
2167                 profzabs->Fill(irow,dzabs);
2168           }
2169         }
2170       }
2171       //
2172       //
2173       // Fill laser fit histograms
2174       //
2175       Float_t dedx = track->GetdEdx();
2176       if (nclI>20&&nclO>20){
2177         fHisNclIn->Fill(id,nclI);      //->Number of clusters inner
2178         fHisNclOut->Fill(id,nclO);     //->Number of clusters outer
2179         fHisNclIO->Fill(id,nclI+nclO);      //->Number of cluster inner outer
2180         //
2181         fHisLclIn->Fill(id,xinMax-xinMin);        //->Level arm inner
2182         fHisLclOut->Fill(id,xoutMax-xoutMin);     //->Level arm outer
2183         fHisLclIO->Fill(id,xoutMax-xinMin);       //->Number of cluster inner outer
2184         //
2185         fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
2186         fHisdZfit->Fill(id,fFitZ[id]);
2187         //
2188         //
2189         fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y));      //->chi2 y inner - line
2190         fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y));     //->chi2 y inner - line
2191         fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y));      //->chi2 y inner - parabola
2192         fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y));     //->chi2 y inner - parabola
2193         fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy));      //->chi2 y IO    - common
2194
2195
2196         fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z));      //->chi2 z inner - line
2197         fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z));     //->chi2 z inner - line
2198         fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y));      //->chi2 z inner - parabola
2199         fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z));     //->chi2 z inner - parabola
2200         fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz));      //->chi2 z IO    - common
2201         //
2202         //
2203         fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]);     //-> delta y   P0outer-P0inner - line
2204         fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]);     //-> delta y   P0outer-P0inner - parabola
2205         fHisPy3vP0->Fill(id,vecy4res[1]);                           //-> delta y   P0outer-P0inner - common parabola
2206         //
2207         fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]);     //-> delta ky  P1outer-P1inner - line
2208         fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]);     //-> delta ky  P1outer-P1inner - parabola
2209         fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]);               //-> delta ky  P1outer-P1inner - common parabola
2210         //
2211         fHisPy3vP2IO->Fill(id,vecy4res[4]);        //-> Curv  P2outerinner - common parabola
2212         fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]);     //-> delta z   P0outer-P0inner - line
2213         fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]);     //-> delta z   P0outer-P0inner - parabola
2214         fHisPz3vP0->Fill(id,vecz4res[1]);                           //-> delta z   P0outer-P0inner - common parabola
2215         //
2216         fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]);     //-> delta kz  P1outer-P1inner - line
2217         fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]);     //-> delta kz  P1outer-P1inner - parabola
2218         fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]);               //-> delta kz  P1outer-P1inner - common parabola
2219         fHisPz3vP2IO->Fill(id,vecz4res[4]);        //-> Curv  P2outerinner - common parabola
2220       }
2221       if(nclI>20){
2222         fHisPy2vP2In->Fill(id,vecy2resInner[2]);   //-> Curv  P2inner - parabola
2223         fHisPz2vP2In->Fill(id,vecz2resInner[2]);   //-> Curv  P2inner - parabola
2224       }
2225       //
2226       if (nclO>20){
2227         fHisPz2vP2Out->Fill(id,vecz2resOuter[2]);  //-> Curv  P2outer - parabola
2228         fHisPy2vP2Out->Fill(id,vecy2resOuter[2]);  //-> Curv  P2outer - parabola
2229       }
2230
2231   }
2232 }
2233
2234
2235
2236 void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
2237   //
2238   //  Dump information about laser beams
2239   //  isOK variable indicates usability of the beam
2240   //  Beam is not usable if:
2241   //     a.  No entries in range (krmsCut0)
2242   //     b.  Big sperad          (krmscut1)
2243   //     c.  RMSto fit sigma bigger then (kmultiCut)
2244   //     d.  Too big angular spread
2245   //
2246
2247   const Float_t krmsCut0=0.001;
2248   const Float_t krmsCut1=0.16;
2249   const Float_t kmultiCut=2;
2250   const Float_t kcutP0=0.002;
2251   //
2252   AliTPCcalibLaser *laser = this;
2253   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
2254   TF1 fg("fg","gaus");
2255   AliTPCParam  * tpcparam    = 0;
2256   // start set up for absolute residuals analysis
2257   //
2258   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
2259   tpcparam    = calib->GetParameters();
2260   if (!tpcparam) tpcparam    = new AliTPCParamSR;
2261   tpcparam->Update();
2262   AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
2263   Float_t current=0;
2264   Float_t bfield      = 0, bz=0;
2265   if (grp){
2266     current = grp->GetL3Current((AliGRPObject::Stats)0);
2267     bfield = 5*current/30000.;
2268     bz = 5*current/30000.;
2269     printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
2270   }
2271
2272   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2273   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2274   TLinearFitter lfabsyInner(2);
2275   lfabsyInner.SetFormula("1 ++ x");
2276   TLinearFitter lfabszInner(2);
2277   lfabszInner.SetFormula("1 ++ x");
2278
2279   TLinearFitter lfabsyOuter(2);
2280   lfabsyOuter.SetFormula("1 ++ x");
2281   TLinearFitter lfabszOuter(2);
2282   lfabszOuter.SetFormula("1 ++ x");
2283   // end set up for absolute residuals analysis
2284
2285   //
2286   //
2287   for (Int_t id=0; id<336; id++){
2288     Bool_t isOK=kTRUE;
2289     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
2290     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
2291     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
2292     TH1F * hisP3    = (TH1F*)laser->fDeltaP3.At(id);
2293     TH1F * hisP4    = (TH1F*)laser->fDeltaP4.At(id);
2294     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
2295     //if (!hisphi) continue;
2296     Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
2297     //if (entries<minEntries) continue;
2298     //
2299     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
2300     if (!ltrp) {
2301      AliTPCLaserTrack::LoadTracks();
2302       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2303     }
2304     pcstream->GetFile()->cd();
2305     if (hisphi)  hisphi->Write();
2306     if (hisphiP) hisphiP->Write();
2307     if (hisZ)    hisZ->Write();
2308     if (hisP3)    hisP3->Write();
2309     if (hisP4)    hisP4->Write();
2310
2311     Float_t meanphi = hisphi->GetMean();
2312     Float_t rmsphi = hisphi->GetRMS();
2313     //
2314     Float_t meanphiP = hisphiP->GetMean();
2315     Float_t rmsphiP = hisphiP->GetRMS();
2316     Float_t meanZ = hisZ->GetMean();
2317     Float_t rmsZ = hisZ->GetRMS();
2318     if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
2319       hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
2320     Double_t gphi1 = fg.GetParameter(1);
2321     Double_t gphi2 = fg.GetParameter(2);
2322     if (hisphiP->GetRMS()>0)
2323       hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
2324     Double_t gphiP1 = fg.GetParameter(1);
2325     Double_t gphiP2 = fg.GetParameter(2);
2326     //
2327     if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
2328       hisZ->Fit(&fg,"","");
2329     Double_t gz1 = fg.GetParameter(1);
2330     Double_t gz2 = fg.GetParameter(2);
2331     //
2332     if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
2333       hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
2334     Double_t gp31 = fg.GetParameter(1);
2335     Double_t gp32 = fg.GetParameter(2);
2336     //
2337     if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2338       hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2339     Double_t gp41 = fg.GetParameter(1);
2340     Double_t gp42 = fg.GetParameter(2);
2341     //
2342     Float_t meanS=hisS->GetMean();
2343     //
2344     Double_t lxyz[3];
2345     Double_t lpxyz[3];
2346     ltrp->GetXYZ(lxyz);
2347     ltrp->GetPxPyPz(lpxyz);
2348
2349     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2350     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2351     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
2352     if (gphiP2>kcutP0) isOK=kFALSE;
2353     //
2354     //
2355     //
2356     TH1 *his =0;
2357     //
2358     his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2359     Float_t mnclIn =  his->GetMean();
2360     delete his;
2361     his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2362     Float_t mnclOut =  his->GetMean();
2363     delete his;
2364     his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2365     Float_t mnclIO =  his->GetMean();
2366     delete his;
2367     his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2368     Float_t mLclIn =  his->GetMean();
2369     delete his;
2370     his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2371     Float_t mLclOut =  his->GetMean();
2372     delete his;
2373     his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2374     Float_t mLclIO =  his->GetMean();
2375     delete his;
2376     //
2377     his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2378     Float_t mdEdx =  his->GetMean();
2379     delete his;
2380     //
2381     //
2382     //
2383     //
2384     his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - line
2385     Float_t  mChi2YIn1=  his->GetMean();
2386     delete his;
2387     his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - line
2388     Float_t mChi2YOut1 =  his->GetMean();
2389     delete his;
2390     his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - parabola
2391     Float_t mChi2YIn2 =  his->GetMean();
2392     delete his;
2393     his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - parabola
2394     Float_t mChi2YOut2 =  his->GetMean();
2395     delete his;
2396     his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 y IO    - common
2397     Float_t mChi2YIO1 =  his->GetMean();
2398     delete his;
2399     his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - line
2400     Float_t mChi2ZIn1 =  his->GetMean();
2401     delete his;
2402     his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - line
2403     Float_t mChi2ZOut1 =  his->GetMean();
2404     delete his;
2405     his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - parabola
2406     Float_t mChi2ZIn2 =  his->GetMean();
2407     delete his;
2408     his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - parabola
2409     Float_t mChi2ZOut2 =  his->GetMean();
2410     delete his;
2411     his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 z IO    - common
2412     Float_t mChi2ZIO1 =  his->GetMean();
2413     delete his;
2414     //
2415     // fit res. histos
2416     //
2417     his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2418     Float_t  edZfit  = his->GetEntries();
2419     Float_t  mdZfit =  his->GetMean();
2420     Float_t  rdZfit  = his->GetRMS();
2421     delete his;
2422
2423     his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - line
2424     Float_t  ePy1vP0  = his->GetEntries();
2425     Float_t  mPy1vP0  = his->GetMean();
2426     Float_t  rPy1vP0  = his->GetRMS();
2427     delete his;
2428
2429     his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - parabola
2430     Float_t  ePy2vP0  = his->GetEntries();
2431     Float_t  mPy2vP0  = his->GetMean();
2432     Float_t  rPy2vP0  = his->GetRMS();
2433     delete his;
2434
2435     his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - common parabola
2436     Float_t  ePy3vP0  = his->GetEntries();
2437     Float_t  mPy3vP0  = his->GetMean();
2438     Float_t  rPy3vP0  = his->GetRMS();
2439     delete his;
2440
2441     his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - line
2442     Float_t  ePy1vP1  = his->GetEntries();
2443     Float_t  mPy1vP1  = his->GetMean();
2444     Float_t  rPy1vP1  = his->GetRMS();
2445     delete his;
2446
2447     his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - parabola
2448     Float_t  ePy2vP1  = his->GetEntries();
2449     Float_t  mPy2vP1  = his->GetMean();
2450     Float_t  rPy2vP1  = his->GetRMS();
2451     delete his;
2452
2453     his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - common parabola
2454     Float_t  ePy3vP1  = his->GetEntries();
2455     Float_t  mPy3vP1  = his->GetMean();
2456     Float_t  rPy3vP1  = his->GetRMS();
2457     delete his;
2458
2459     his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2460     Float_t  ePy2vP2In  = his->GetEntries();
2461     Float_t  mPy2vP2In  = his->GetMean();
2462     Float_t  rPy2vP2In  = his->GetRMS();
2463     delete his;
2464
2465     his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2466     Float_t  ePy2vP2Out  = his->GetEntries();
2467     Float_t  mPy2vP2Out  = his->GetMean();
2468     Float_t  rPy2vP2Out  = his->GetRMS();
2469     delete his;
2470
2471     his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2472     Float_t  ePy3vP2IO  = his->GetEntries();
2473     Float_t  mPy3vP2IO  = his->GetMean();
2474     Float_t  rPy3vP2IO  = his->GetRMS();
2475     delete his;
2476
2477     //
2478     //
2479     his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - line
2480     Float_t  ePz1vP0  = his->GetEntries();
2481     Float_t  mPz1vP0  = his->GetMean();
2482     Float_t  rPz1vP0  = his->GetRMS();
2483     delete his;
2484
2485     his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - parabola
2486     Float_t  ePz2vP0  = his->GetEntries();
2487     Float_t  mPz2vP0  = his->GetMean();
2488     Float_t  rPz2vP0  = his->GetRMS();
2489     delete his;
2490
2491     his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - common parabola
2492     Float_t  ePz3vP0  = his->GetEntries();
2493     Float_t  mPz3vP0  = his->GetMean();
2494     Float_t  rPz3vP0  = his->GetRMS();
2495     delete his;
2496
2497     his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - line
2498     Float_t  ePz1vP1  = his->GetEntries();
2499     Float_t  mPz1vP1  = his->GetMean();
2500     Float_t  rPz1vP1  = his->GetRMS();
2501     delete his;
2502
2503     his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - parabola
2504     Float_t  ePz2vP1  = his->GetEntries();
2505     Float_t  mPz2vP1  = his->GetMean();
2506     Float_t  rPz2vP1  = his->GetRMS();
2507     delete his;
2508
2509     his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - common parabola
2510     Float_t  ePz3vP1  = his->GetEntries();
2511     Float_t  mPz3vP1  = his->GetMean();
2512     Float_t  rPz3vP1  = his->GetRMS();
2513     delete his;
2514
2515     his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2516     Float_t  ePz2vP2In  = his->GetEntries();
2517     Float_t  mPz2vP2In  = his->GetMean();
2518     Float_t  rPz2vP2In  = his->GetRMS();
2519     delete his;
2520
2521     his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2522     Float_t  ePz2vP2Out  = his->GetEntries();
2523     Float_t  mPz2vP2Out  = his->GetMean();
2524     Float_t  rPz2vP2Out  = his->GetRMS();
2525     delete his;
2526
2527     his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2528     Float_t  ePz3vP2IO  = his->GetEntries();
2529     Float_t  mPz3vP2IO  = his->GetMean();
2530     Float_t  rPz3vP2IO  = his->GetRMS();
2531     delete his;
2532
2533     // Fit absolute laser residuals
2534     TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2535     TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2536
2537     Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2538     Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2539
2540     TVectorD vecX(159);        // X
2541     TVectorD vecY(159);        // Y
2542     TVectorD vecR(159);        // R
2543     TVectorD vecDY(159);       // absolute residuals in Y
2544     TVectorD vecDZ(159);       // absolute residuals in Z
2545     TVectorD vecN(159);        // number of clusters
2546     TVectorD vecEy(159);       //error y
2547     TVectorD vecEz(159);       //error z
2548     TVectorD vecPhi(159);      // local tangent
2549     TVectorD vecPhiR(159);     // local tangent
2550
2551     lfabsyInner.ClearPoints();
2552     lfabszInner.ClearPoints();
2553     lfabsyOuter.ClearPoints();
2554     lfabszOuter.ClearPoints();
2555     // dummy fit values
2556     Int_t    nClInY         = 0;
2557     Float_t  yAbsInOffset  = -100;
2558     Float_t  yAbsInSlope   = -100;
2559     Float_t  yAbsInDeltay  = -100;
2560     Int_t    nClInZ         = 0;
2561     Float_t  zAbsInOffset  = -100;
2562     Float_t  zAbsInSlope   = -100;
2563     Float_t  zAbsInDeltaz  = -100;
2564     Int_t    nClOutY         = 0;
2565     Float_t  yAbsOutOffset  = -100;
2566     Float_t  yAbsOutSlope   = -100;
2567     Float_t  yAbsOutDeltay  = -100;
2568     Int_t    nClOutZ         = 0;
2569     Float_t  zAbsOutOffset  = -100;
2570     Float_t  zAbsOutSlope   = -100;
2571     Float_t  zAbsOutDeltaz  = -100;
2572
2573     Float_t  lasTanPhiLocIn = -100;
2574     Float_t  lasTanPhiLocOut = -100;
2575
2576     if(histAbsY->GetEntries()>0) {
2577
2578       Double_t rotAngOut = 10;
2579       Double_t rotAngIn = 10;
2580       if((secInner%36)!=(secOuter%36))
2581         rotAngIn += 20; // 30 degrees
2582
2583       // Calculate laser mirror X position in local frame
2584       Double_t laserposOut =
2585         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2586       Double_t laserposIn =
2587         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2588
2589       // Calculate laser tan phi in local frame
2590       lasTanPhiLocOut  = TMath::ASin(ltrp->GetSnp());
2591       if(lasTanPhiLocOut<0) {
2592         lasTanPhiLocIn  = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2593         lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2594       } else {
2595
2596         lasTanPhiLocIn  = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2597         lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2598       }
2599
2600       lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2601       lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2602
2603       TProfile* yprof  = histAbsY->ProfileX("yprof");
2604       TProfile* zprof  = histAbsZ->ProfileX("zprof");
2605
2606       for(Int_t bin = 1; bin<=159; bin++) {
2607
2608         if(yprof->GetBinEntries(bin)<10&&
2609            zprof->GetBinEntries(bin)<10) {
2610           continue;
2611         }
2612
2613         // There is a problem in defining inner and outer sectors for
2614         // the outer beams (0 and 6) where both sectors are OROCs. To
2615         // make sure there is no overlap row 94 to 99 are cutted.
2616         if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2617           continue;
2618
2619         Int_t row = (bin-1);
2620         if(row>62)
2621           row -= 63;
2622
2623         Bool_t isOuter = kTRUE;
2624         Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2625
2626         if(bin<=62 ||
2627            (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2628
2629           isOuter = kFALSE;
2630           sector = TMath::Nint(fBeamSectorInner[id]);
2631         }
2632
2633
2634         Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2635         vecN[bin-1] =yprof->GetBinEntries(bin);
2636         vecEy[bin-1]=yprof->GetBinError(bin);
2637         vecEz[bin-1]=zprof->GetBinError(bin);
2638         vecX[bin-1] = x;
2639         vecDY[bin-1] = yprof->GetBinContent(bin);
2640         vecDZ[bin-1] = zprof->GetBinContent(bin);
2641         if(!isOuter) { // inner
2642           vecPhi[bin-1]=lasTanPhiLocIn;
2643           // Calculate local y from residual and database
2644           Double_t y  = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2645             + vecDY[bin-1];
2646           vecY[bin-1] = y;
2647           Double_t r  = TMath::Sqrt(x*x + y*y);
2648           vecR[bin-1] = r;
2649           // Find angle between laser vector and R vector
2650           // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2651           Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2652           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2653           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2654           if(lasTanPhiLocIn<0)
2655             vecPhiR[bin-1]*=-1; // to have the same sign
2656
2657           if(yprof->GetBinEntries(bin)>=10) {
2658             lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2659                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2660           }
2661           if(zprof->GetBinEntries(bin)>=10) {
2662             lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2663                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2664           }
2665         } else { // outer
2666           vecPhi[bin-1]=lasTanPhiLocOut;
2667           // Calculate local y from residual and database
2668           Double_t y  = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2669             + vecDY[bin-1];
2670           vecY[bin-1] = y;
2671           Double_t r  = TMath::Sqrt(x*x + y*y);
2672           vecR[bin-1] = r;
2673
2674           Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2675           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2676           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2677           if(lasTanPhiLocOut<0)
2678             vecPhiR[bin-1]*=-1; // to have the same sign
2679
2680           if(yprof->GetBinEntries(bin)>=10) {
2681             lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2682                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2683           }
2684           if(zprof->GetBinEntries(bin)>=10) {
2685             lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2686                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2687           }
2688         }
2689
2690       }
2691
2692       delete yprof; delete zprof;
2693
2694
2695       // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2696       nClInY = lfabsyInner.GetNpoints();
2697       if(lfabsyInner.GetNpoints()>10) {
2698         lfabsyInner.EvalRobust(0.95);
2699
2700         TVectorD result(2);
2701         lfabsyInner.GetParameters(result);
2702         yAbsInOffset = result[0];
2703         yAbsInSlope  = result[1];
2704         yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2705       }
2706       nClInZ = lfabszInner.GetNpoints();
2707       if(lfabszInner.GetNpoints()>10) {
2708         lfabszInner.EvalRobust(0.95);
2709
2710         TVectorD result(2);
2711         lfabszInner.GetParameters(result);
2712         zAbsInOffset = result[0];
2713         zAbsInSlope  = result[1];
2714         zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2715       }
2716       nClOutY = lfabsyOuter.GetNpoints();
2717       if(lfabsyOuter.GetNpoints()>10) {
2718         lfabsyOuter.EvalRobust(0.95);
2719
2720         TVectorD result(2);
2721         lfabsyOuter.GetParameters(result);
2722         yAbsOutOffset = result[0];
2723         yAbsOutSlope  = result[1];
2724         yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2725       }
2726       nClOutZ = lfabszOuter.GetNpoints();
2727       if(lfabszOuter.GetNpoints()>10) {
2728         lfabszOuter.EvalRobust(0.95);
2729
2730         TVectorD result(2);
2731         lfabszOuter.GetParameters(result);
2732         zAbsOutOffset = result[0];
2733         zAbsOutSlope  = result[1];
2734         zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2735       }
2736     }
2737
2738
2739     Int_t itime=-1;
2740     Float_t  coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2741     Float_t  coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2742     Float_t  coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2743     Float_t  coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2744     //
2745     Float_t  skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2746     Float_t  skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2747     //
2748     Float_t  ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2749     Float_t  ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2750
2751     //
2752     (*pcstream)<<"Mean"<<
2753       "run="<<run<<               //
2754       //voltages
2755       "UcIA="<<coverIA<<
2756       "UcIC="<<coverIC<<
2757       "UcOA="<<coverOA<<
2758       "UcOC="<<coverOC<<
2759       "UsA="<<skirtA<<
2760       "UsC="<<skirtC<<
2761       "UggA="<<ggOffA<<
2762       "UggC="<<ggOffC<<
2763       //
2764       "isOK="<<isOK<<             //
2765       "id="<<id<<                 // track id
2766       "entries="<<entries<<       // number of entries
2767       "bz="<<bfield<<             // bfield
2768       "LTr.="<<ltrp<<             // refernece track
2769       //
2770       "mnclIn="<<mnclIn<<         // mean number of clusters in inner
2771       "mnclOut="<<mnclOut<<       // mean number of clusters in outer
2772       "mnclIO="<<mnclIO<<         // mean number of clusters in inner+outer
2773       "mLclIn="<<mLclIn<<         // mean number of clusters in inner
2774       "mLclOut="<<mLclOut<<       // mean number of clusters in outer
2775       "mLclIO="<<mLclIO<<         // mean number of clusters in inner+outer
2776       "mdEdx="<<mdEdx<<           // mean dEdx
2777       "edZfit="<<edZfit<<           // entries z fit
2778       "mdZfit="<<mdZfit<<           // mean z fit
2779       "rdZfit="<<rdZfit<<           // RMS z fit
2780       //
2781       //
2782       "mChi2YIn1="<<mChi2YIn1<<       //->chi2 y inner - line
2783       "mChi2YOut1="<<mChi2YOut1<<     //->chi2 y inner - line
2784       "mChi2YIn2="<<mChi2YIn2<<       //->chi2 y inner - parabola
2785       "mChi2YOut2="<<mChi2YOut2<<     //->chi2 y inner - parabola
2786       "mChi2YIO1="<<mChi2YIO1<<       //->chi2 y IO    - common
2787       "mChi2ZIn1="<<mChi2ZIn1<<       //->chi2 z inner - line
2788       "mChi2ZOut1="<<mChi2ZOut1<<     //->chi2 z inner - line
2789       "mChi2ZIn2="<<mChi2ZIn2<<       //->chi2 z inner - parabola
2790       "mChi2ZOut2="<<mChi2ZOut2<<     //->chi2 z inner - parabola
2791       "mChi2ZIO1="<<mChi2ZIO1<<       //->chi2 z IO    - common
2792       //
2793       //
2794       "ePy1vP0="<<ePy1vP0<<
2795       "mPy1vP0="<<mPy1vP0<<
2796       "rPy1vP0="<<rPy1vP0<<
2797       "ePy2vP0="<<ePy2vP0<<
2798       "mPy2vP0="<<mPy2vP0<<
2799       "rPy2vP0="<<rPy2vP0<<
2800       "ePy3vP0="<<ePy3vP0<<
2801       "mPy3vP0="<<mPy3vP0<<
2802       "rPy3vP0="<<rPy3vP0<<
2803       "ePy1vP1="<<ePy1vP1<<
2804       "mPy1vP1="<<mPy1vP1<<
2805       "rPy1vP1="<<rPy1vP1<<
2806       "ePy2vP1="<<ePy2vP1<<
2807       "mPy2vP1="<<mPy2vP1<<
2808       "rPy2vP1="<<rPy2vP1<<
2809       "ePy3vP1="<<ePy3vP1<<
2810       "mPy3vP1="<<mPy3vP1<<
2811       "rPy3vP1="<<rPy3vP1<<
2812       "ePy2vP2In="<<ePy2vP2In<<
2813       "mPy2vP2In="<<mPy2vP2In<<
2814       "rPy2vP2In="<<rPy2vP2In<<
2815       "ePy2vP2Out="<<ePy2vP2Out<<
2816       "mPy2vP2Out="<<mPy2vP2Out<<
2817       "rPy2vP2Out="<<rPy2vP2Out<<
2818       "ePy3vP2IO="<<ePy3vP2IO<<
2819       "mPy3vP2IO="<<mPy3vP2IO<<
2820       "rPy3vP2IO="<<rPy3vP2IO<<
2821       //
2822       //
2823       "ePz1vP0="<<ePz1vP0<<
2824       "mPz1vP0="<<mPz1vP0<<
2825       "rPz1vP0="<<rPz1vP0<<
2826       "ePz2vP0="<<ePz2vP0<<
2827       "mPz2vP0="<<mPz2vP0<<
2828       "rPz2vP0="<<rPz2vP0<<
2829       "ePz3vP0="<<ePz3vP0<<
2830       "mPz3vP0="<<mPz3vP0<<
2831       "rPz3vP0="<<rPz3vP0<<
2832       "ePz1vP1="<<ePz1vP1<<
2833       "mPz1vP1="<<mPz1vP1<<
2834       "rPz1vP1="<<rPz1vP1<<
2835       "ePz2vP1="<<ePz2vP1<<
2836       "mPz2vP1="<<mPz2vP1<<
2837       "rPz2vP1="<<rPz2vP1<<
2838       "ePz3vP1="<<ePz3vP1<<
2839       "mPz3vP1="<<mPz3vP1<<
2840       "rPz3vP1="<<rPz3vP1<<
2841       "ePz2vP2In="<<ePz2vP2In<<
2842       "mPz2vP2In="<<mPz2vP2In<<
2843       "rPz2vP2In="<<rPz2vP2In<<
2844       "ePz2vP2Out="<<ePz2vP2Out<<
2845       "mPz2vP2Out="<<mPz2vP2Out<<
2846       "rPz2vP2Out="<<rPz2vP2Out<<
2847       "ePz3vP2IO="<<ePz3vP2IO<<
2848       "mPz3vP2IO="<<mPz3vP2IO<<
2849       "rPz3vP2IO="<<rPz3vP2IO<<
2850       //
2851       //
2852       //
2853       "lx0="<<lxyz[0]<<            // reference x
2854       "lx1="<<lxyz[1]<<            // reference y
2855       "lx2="<<lxyz[2]<<            // refernece z
2856       "lpx0="<<lpxyz[0]<<           // reference x
2857       "lpx1="<<lpxyz[1]<<          // reference y
2858       "lpx2="<<lpxyz[2]<<          // refernece z
2859       //
2860       "msig="<<meanS<<
2861       //
2862       "mphi="<<meanphi<<         //
2863       "rmsphi="<<rmsphi<<        //
2864       "gphi1="<<gphi1<<
2865       "gphi2="<<gphi2<<
2866       //
2867       "mphiP="<<meanphiP<<       //
2868       "rmsphiP="<<rmsphiP<<      //
2869       "gphiP1="<<gphiP1<<
2870       "gphiP2="<<gphiP2<<
2871       //
2872       "meanZ="<<meanZ<<
2873       "rmsZ="<<rmsZ<<
2874       "gz1="<<gz1<<
2875       "gz2="<<gz2<<
2876       //
2877       "gp31="<<gp31<<            //gaus mean - tgl
2878       "gp32="<<gp32<<            //gaus rms  -tgl
2879       "gp41="<<gp41<<            //gaus mean - P4
2880       "gp42="<<gp42<<            //gaus rms  - P4
2881       // Parameters from abs res analysis
2882       "SecIn="<<secInner<<              // inner sector
2883       "SecOut="<<secOuter<<             // outer sector
2884       "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2885       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
2886       "X.="<<&vecX<<       // local x
2887       "Y.="<<&vecY<<       // local y
2888       "R.="<<&vecR<<       // radius
2889       "dY.="<<&vecDY<<     // abs y residuals
2890       "dZ.="<<&vecDZ<<     // abs z residuals
2891       "eY.="<<&vecEy<<     // error of y residuals
2892       "eZ.="<<&vecEz<<     // error of z residuals
2893       "kY.="<<&vecPhi<<    // local tangent y (fixed for sector)
2894       "kYR.="<<&vecPhiR<<  // tangent between laser and R vector (varies inside sector)
2895       "nCl.="<<&vecN<<     // number of clusters
2896       //
2897       "nClInY="<<nClInY<<               // Number of clusters for inner y
2898       "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2899       "yAbsInSlope="<<yAbsInSlope <<  // fitted slope absres (inner y)
2900       "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2901       "nClInZ="<<nClInZ<<               // Number of clusters for inner z
2902       "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2903       "zAbsInSlope="<<zAbsInSlope <<  // fitted slope absres (inner z)
2904       "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2905       //
2906       "nClOutY="<<nClOutY<<               // Number of clusters for outer y
2907       "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
2908       "yAbsOutSlope="<<yAbsOutSlope <<  // fitted slope absres (outer y)
2909       "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
2910       "nClOutZ="<<nClOutZ<<               // Number of clusters for outer z
2911       "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
2912       "zAbsOutSlope="<<zAbsOutSlope <<  // fitted slope absres (outer z)
2913       "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
2914       //
2915       "\n";
2916   }
2917   delete pcstream;
2918   /*
2919     Browse the content
2920     TFile fmean("laserMean.root");
2921
2922
2923    */
2924
2925
2926 }
2927
2928
2929
2930 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
2931   //
2932   //
2933   //
2934   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
2935   TFile * f = pcstream->GetFile();
2936   f->mkdir("dirphi");
2937   f->mkdir("dirphiP");
2938   f->mkdir("dirZ");
2939   TF1 fp("p1","pol1");
2940   //
2941   //
2942   char    cut[1000];
2943   char grname[1000];
2944   char grnamefull[1000];
2945
2946   Double_t mphi[100];
2947   Double_t mphiP[100];
2948   Double_t smphi[100];
2949   Double_t smphiP[100];
2950   Double_t mZ[100];
2951   Double_t smZ[100];
2952   Double_t bz[100];
2953   Double_t sbz[100];
2954   // fit parameters
2955   Double_t pphi[3];
2956   Double_t pphiP[3];
2957   Double_t pmZ[3];
2958
2959   //
2960   for (Int_t id=0; id<336; id++){
2961     // id =205;
2962     sprintf(cut,"fId==%d&&%s",id,cutUser);
2963     Int_t entries = chain->Draw("bz",cut,"goff");
2964     if (entries<3) continue;
2965     AliTPCLaserTrack *ltrp = 0;
2966     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
2967     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2968     Double_t lxyz[3];
2969     Double_t lpxyz[3];
2970     ltrp->GetXYZ(lxyz);
2971     ltrp->GetPxPyPz(lpxyz);
2972
2973     chain->Draw("bz",cut,"goff");
2974     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
2975     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
2976     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
2977     //
2978     chain->Draw("gphi1",cut,"goff");
2979     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
2980     chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
2981     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
2982     //
2983     chain->Draw("gphiP1",cut,"goff");
2984     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
2985     chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
2986     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
2987     //
2988     chain->Draw("gz1",cut,"goff");
2989     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
2990     chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
2991     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
2992     //
2993     //
2994     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
2995             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
2996     // store data
2997     // phi
2998     f->cd("dirphi");
2999     Float_t phiB0 =0;
3000     Float_t phiPB0=0;
3001     Float_t zB0=0;
3002     for (Int_t ientry=0; ientry<entries; ientry++){
3003       if (TMath::Abs(bz[ientry])<0.05){
3004         phiB0  = mphi[ientry];
3005         phiPB0 = mphiP[ientry];
3006         zB0    = mZ[ientry];
3007       }
3008     }
3009     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3010     grphi->Draw("a*");
3011     grphi->Fit(&fp);
3012     pphi[0] = fp.GetParameter(0);                          // offset
3013     pphi[1] = fp.GetParameter(1);                          // slope
3014     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3015     sprintf(grname,"phi_id%d",id);
3016     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
3017     grphi->GetXaxis()->SetTitle("b_{z} (T)");
3018     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3019     grphi->SetMaximum(1.2);
3020     grphi->SetMinimum(-1.2);
3021     grphi->Draw("a*");
3022
3023     grphi->Write();
3024     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3025     // phiP
3026     f->cd("dirphiP)");
3027     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3028     grphiP->Draw("a*");
3029     grphiP->Fit(&fp);
3030     pphiP[0] = fp.GetParameter(0);                          // offset
3031     pphiP[1] = fp.GetParameter(1);                          // slope
3032     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3033     sprintf(grname,"phiP_id%d",id);
3034     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
3035     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3036     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3037     grphiP->SetMaximum(pphiP[0]+0.005);
3038     grphiP->SetMinimum(pphiP[0]-0.005);
3039
3040     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3041     grphiP->Write();
3042     //
3043     //Z
3044     f->cd("dirZ");
3045     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3046     grmZ->Draw("a*");
3047     grmZ->Fit(&fp);
3048     pmZ[0] = fp.GetParameter(0);                          // offset
3049     pmZ[1] = fp.GetParameter(1);                          // slope
3050     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3051     sprintf(grname,"mZ_id%d",id);
3052     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
3053     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3054     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3055
3056     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3057     grmZ->Write();
3058     //
3059     // P4
3060     //
3061
3062     for (Int_t ientry=0; ientry<entries; ientry++){
3063       (*pcstream)<<"Mean"<<
3064         "id="<<id<<
3065         "LTr.="<<ltrp<<
3066         "entries="<<entries<<
3067         "bz="<<bz[ientry]<<
3068         "lx0="<<lxyz[0]<<          // reference x
3069         "lx1="<<lxyz[1]<<          // reference y
3070         "lx2="<<lxyz[2]<<          // refernece z
3071         "lpx0="<<lpxyz[0]<<          // reference x
3072         "lpx1="<<lpxyz[1]<<          // reference y
3073         "lpx2="<<lpxyz[2]<<          // refernece z
3074         //values
3075         "phiB0="<<phiB0<<          // position shift at 0 field
3076         "phiPB0="<<phiPB0<<        // angular  shift at 0 field
3077         "zB0="<<zB0<<              // z shift for 0 field
3078         //
3079         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3080         "pphi0="<<pphi[0]<<   // offset
3081         "pphi1="<<pphi[1]<<   // slope
3082         "pphi2="<<pphi[2]<<   // norm chi2
3083         //
3084         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3085         "pphiP0="<<pphiP[0]<< // offset
3086         "pphiP1="<<pphiP[1]<< // slope
3087         "pphiP2="<<pphiP[2]<< // norm chi2
3088         //
3089         "gz1="<<mZ[ientry]<<
3090         "pmZ0="<<pmZ[0]<<     // offset
3091         "pmZ1="<<pmZ[1]<<     // slope
3092         "pmZ2="<<pmZ[2]<<     // norm chi2
3093         "\n";
3094     }
3095   }
3096
3097   delete pcstream;
3098
3099 }
3100
3101
3102 void AliTPCcalibLaser::Analyze(){
3103   //
3104   //
3105   //
3106 }
3107
3108
3109 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3110
3111   TIterator* iter = li->MakeIterator();
3112   AliTPCcalibLaser* cal = 0;
3113   static Int_t counter0=0;
3114   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3115     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3116       //      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3117       return -1;
3118     }
3119     printf("Marging number %d\n", counter0);
3120     counter0++;
3121     //
3122     MergeFitHistos(cal);
3123     TH1F *h=0x0;
3124     TH1F *hm=0x0;
3125     TH2F *h2=0x0;
3126     TH2F *h2m=0x0;
3127     //    TProfile *hp=0x0;
3128     //TProfile *hpm=0x0;
3129
3130     for (Int_t id=0; id<336; id++){
3131       // merge fDeltaZ histograms
3132       hm = (TH1F*)cal->fDeltaZ.At(id);
3133       h  = (TH1F*)fDeltaZ.At(id);
3134       if (!h) {
3135         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3136         h->SetDirectory(0);
3137         fDeltaZ.AddAt(h,id);
3138       }
3139       if (hm) h->Add(hm);
3140       // merge fP3 histograms
3141       hm = (TH1F*)cal->fDeltaP3.At(id);
3142       h  = (TH1F*)fDeltaP3.At(id);
3143       if (!h) {
3144         h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3145         h->SetDirectory(0);
3146         fDeltaP3.AddAt(h,id);
3147       }
3148       if (hm) h->Add(hm);
3149       // merge fP4 histograms
3150       hm = (TH1F*)cal->fDeltaP4.At(id);
3151       h  = (TH1F*)fDeltaP4.At(id);
3152       if (!h) {
3153         h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3154         h->SetDirectory(0);
3155         fDeltaP4.AddAt(h,id);
3156       }
3157       if (hm) h->Add(hm);
3158
3159       //
3160       // merge fDeltaPhi histograms
3161       hm = (TH1F*)cal->fDeltaPhi.At(id);
3162       h  = (TH1F*)fDeltaPhi.At(id);
3163       if (!h) {
3164         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3165         h->SetDirectory(0);
3166         fDeltaPhi.AddAt(h,id);
3167       }
3168       if (hm) h->Add(hm);
3169       // merge fDeltaPhiP histograms
3170       hm = (TH1F*)cal->fDeltaPhiP.At(id);
3171       h  = (TH1F*)fDeltaPhiP.At(id);
3172       if (!h) {
3173         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3174         h->SetDirectory(0);
3175         fDeltaPhiP.AddAt(h,id);
3176       }
3177       if (hm) h->Add(hm);
3178       // merge fSignals histograms
3179       hm = (TH1F*)cal->fSignals.At(id);
3180       h  = (TH1F*)fSignals.At(id);
3181       if (!h) {
3182         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3183         h->SetDirectory(0);
3184         fSignals.AddAt(h,id);
3185       }
3186       if (hm) h->Add(hm);
3187       //
3188       //
3189       // merge ProfileY histograms -0
3190       h2m = (TH2F*)cal->fDeltaYres.At(id);
3191       h2  = (TH2F*)fDeltaYres.At(id);
3192       if (h2m&&h2) h2->Add(h2m);
3193       //
3194       h2m = (TH2F*)cal->fDeltaZres.At(id);
3195       h2  = (TH2F*)fDeltaZres.At(id);
3196       if (h2m&&h2) h2->Add(h2m);
3197       // merge ProfileY histograms - 2
3198       h2m = (TH2F*)cal->fDeltaYres2.At(id);
3199       h2  = (TH2F*)fDeltaYres2.At(id);
3200       if (h2m&&h2) h2->Add(h2m);
3201       //
3202       h2m = (TH2F*)cal->fDeltaZres2.At(id);
3203       h2  = (TH2F*)fDeltaZres2.At(id);
3204       if (h2m&&h2) h2->Add(h2m);
3205
3206       // merge ProfileY histograms - abs
3207       h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3208       h2  = (TH2F*)fDeltaYresAbs.At(id);
3209       if (h2m&&h2) h2->Add(h2m);
3210
3211       h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3212       h2  = (TH2F*)fDeltaZresAbs.At(id);
3213       if (h2m&&h2) h2->Add(h2m);
3214       // merge ProfileY histograms - 3
3215       //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3216       //h2  = (TH2F*)fDeltaYres3.At(id);
3217       //if (h2m) h2->Add(h2m);
3218       //
3219       //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3220       //h2  = (TH2F*)fDeltaZres3.At(id);
3221       //if (h2m) h->Add(h2m);
3222       //
3223       //
3224     }
3225   }
3226   return 0;
3227 }
3228
3229 void   AliTPCcalibLaser::MakeFitHistos(){
3230   //
3231   // Make a fit histograms
3232   //
3233   // Number of clusters
3234   //
3235   //TH2F           *fHisNclIn;      //->Number of clusters inner
3236   //TH2F           *fHisNclOut;     //->Number of clusters outer
3237   //TH2F           *fHisNclIO;      //->Number of cluster inner outer
3238   // TH2F           *fHisdEdx;       //->dEdx histo
3239   fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3240   fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3241   fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3242   //
3243   fHisLclIn  = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3244   fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3245   fHisLclIO  = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3246   //
3247   fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3248   fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3249
3250   //
3251   //  Chi2
3252   //
3253   //   TH2F           *fHisChi2YIn1;      //->chi2 y inner - line
3254   //   TH2F           *fHisChi2YOut1;     //->chi2 y inner - line
3255   //   TH2F           *fHisChi2YIn2;      //->chi2 y inner - parabola
3256   //   TH2F           *fHisChi2YOut2;     //->chi2 y inner - parabola
3257   //   TH2F           *fHisChi2YIO1;      //->chi2 y IO    - common
3258   fHisChi2YIn1   = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3259   fHisChi2YOut1  = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3260   fHisChi2YIn2   = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3261   fHisChi2YOut2  = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3262   fHisChi2YIO1   = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3263   //   TH2F           *fHisChi2ZIn1;      //->chi2 z inner - line
3264   //   TH2F           *fHisChi2ZOut1;     //->chi2 z inner - line
3265   //   TH2F           *fHisChi2ZIn2;      //->chi2 z inner - parabola
3266   //   TH2F           *fHisChi2ZOut2;     //->chi2 z inner - parabola
3267   //   TH2F           *fHisChi2ZIO1;      //->chi2 z IO    - common
3268   fHisChi2ZIn1   = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3269   fHisChi2ZOut1  = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3270   fHisChi2ZIn2   = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3271   fHisChi2ZOut2  = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3272   fHisChi2ZIO1   = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3273   //
3274   // Fit
3275   //
3276   //
3277   //   TH2F           *fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
3278   //   TH2F           *fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
3279   //   TH2F           *fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
3280   //   TH2F           *fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
3281   //   TH2F           *fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
3282   //   TH2F           *fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
3283   //   TH2F           *fHisPy2vP2In;   //-> Curv  P2inner - parabola
3284   //   TH2F           *fHisPy2vP2Out;  //-> Curv  P2outer - parabola
3285   //   TH2F           *fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
3286   fHisPy1vP0   = new TH2F("HisPy1vP0",   "HisPy1vP0",336,0,336,500,-0.3,0.3);
3287   fHisPy2vP0   = new TH2F("HisPy2vP0",   "HisPy2vP0",336,0,336,500,-0.3,0.3);
3288   fHisPy3vP0   = new TH2F("HisPy3vP0",   "HisPy3vP0",336,0,336,500,-0.3,0.3);
3289   fHisPy1vP1   = new TH2F("HisPy1vP1",   "HisPy1vP1",336,0,336,500,-0.01,0.01);
3290   fHisPy2vP1   = new TH2F("HisPy2vP1",   "HisPy2vP1",336,0,336,500,-0.01,0.01);
3291   fHisPy3vP1   = new TH2F("HisPy3vP1",   "HisPy3vP1",336,0,336,500,-0.01,0.01);
3292   fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3293   fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3294   fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3295   //
3296   //
3297   //   TH2F           *fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
3298   //   TH2F           *fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
3299   //   TH2F           *fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
3300   //   TH2F           *fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
3301   //   TH2F           *fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
3302   //   TH2F           *fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
3303   //   TH2F           *fHisPz2vP2In;   //-> Curv  P2inner - parabola
3304   //   TH2F           *fHisPz2vP2Out;  //-> Curv  P2outer - parabola
3305   //   TH2F           *fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
3306   fHisPz1vP0   = new TH2F("HisPz1vP0",   "HisPz1vP0",336,0,336,500,-0.3,0.3);
3307   fHisPz2vP0   = new TH2F("HisPz2vP0",   "HisPz2vP0",336,0,336,500,-0.3,0.3);
3308   fHisPz3vP0   = new TH2F("HisPz3vP0",   "HisPz3vP0",336,0,336,500,-0.3,0.3);
3309   fHisPz1vP1   = new TH2F("HisPz1vP1",   "HisPz1vP1",336,0,336,500,-0.01,0.01);
3310   fHisPz2vP1   = new TH2F("HisPz2vP1",   "HisPz2vP1",336,0,336,500,-0.01,0.01);
3311   fHisPz3vP1   = new TH2F("HisPz3vP1",   "HisPz3vP1",336,0,336,500,-0.01,0.01);
3312   fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3313   fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3314   fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3315
3316   fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3317   fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3318
3319   fHisNclIn->SetDirectory(0);      //->Number of clusters inner
3320   fHisNclOut->SetDirectory(0);     //->Number of clusters outer
3321   fHisNclIO->SetDirectory(0);      //->Number of cluster inner outer
3322   fHisLclIn->SetDirectory(0);      //->Level arm inner
3323   fHisLclOut->SetDirectory(0);     //->Level arm outer
3324   fHisLclIO->SetDirectory(0);      //->Number of cluster inner outer
3325   fHisdEdx->SetDirectory(0);      //->Number of cluster inner outer
3326   fHisdZfit->SetDirectory(0);      //->Number of cluster inner outer
3327   //
3328   //
3329   fHisChi2YIn1->SetDirectory(0);      //->chi2 y inner - line
3330   fHisChi2YOut1->SetDirectory(0);     //->chi2 y inner - line
3331   fHisChi2YIn2->SetDirectory(0);      //->chi2 y inner - parabola
3332   fHisChi2YOut2->SetDirectory(0);     //->chi2 y inner - parabola
3333   fHisChi2YIO1->SetDirectory(0);      //->chi2 y IO    - common
3334   fHisChi2ZIn1->SetDirectory(0);      //->chi2 z inner - line
3335   fHisChi2ZOut1->SetDirectory(0);     //->chi2 z inner - line
3336   fHisChi2ZIn2->SetDirectory(0);      //->chi2 z inner - parabola
3337   fHisChi2ZOut2->SetDirectory(0);     //->chi2 z inner - parabola
3338   fHisChi2ZIO1->SetDirectory(0);      //->chi2 z IO    - common
3339   //
3340   //
3341   fHisPy1vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - line
3342   fHisPy2vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - parabola
3343   fHisPy3vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - common parabola
3344   fHisPy1vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - line
3345   fHisPy2vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - parabola
3346   fHisPy3vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - common parabola
3347   fHisPy2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
3348   fHisPy2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
3349   fHisPy3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
3350   //
3351   //
3352   fHisPz1vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - line
3353   fHisPz2vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - parabola
3354   fHisPz3vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - common parabola
3355   fHisPz1vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - line
3356   fHisPz2vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - parabola
3357   fHisPz3vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - common parabola
3358   fHisPz2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
3359   fHisPz2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
3360   fHisPz3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
3361
3362   fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3363   fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3364
3365
3366
3367   //
3368   //
3369   //
3370   for (Int_t id=0; id<336;id++){
3371     TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3372     TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3373     TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3374     TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
3375     TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3376     TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3377     //    TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3378     //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3379     if (!profy){
3380       profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3381       profy->SetDirectory(0);
3382       fDeltaYres.AddAt(profy,id);
3383       profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3384       profy2->SetDirectory(0);
3385       fDeltaYres2.AddAt(profy2,id);
3386       if(!fUseFixedDriftV)
3387         profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3388       else
3389         profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
3390       profyabs->SetDirectory(0);
3391       fDeltaYresAbs.AddAt(profyabs,id);
3392       //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3393       //profy3->SetDirectory(0);
3394       //fDeltaYres3.AddAt(profy3,id);
3395     }
3396     if (!profz){
3397       profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3398       profz->SetDirectory(0);
3399       fDeltaZres.AddAt(profz,id);
3400       profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3401       profz2->SetDirectory(0);
3402       fDeltaZres2.AddAt(profz2,id);
3403       if(!fUseFixedDriftV)
3404         profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
3405       else
3406         profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
3407       profzabs->SetDirectory(0);
3408       fDeltaZresAbs.AddAt(profzabs,id);
3409       //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3410       //profz3->SetDirectory(0);
3411       //fDeltaZres3.AddAt(profz3,id);
3412     }
3413   }
3414   //
3415   //
3416   for (Int_t id=0; id<336;id++){
3417     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3418     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3419     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3420
3421     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3422     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3423     TH1F * hisSignal = (TH1F*)fSignals.At(id);
3424
3425     if (!hisdz){
3426       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3427       hisdz->SetDirectory(0);
3428       fDeltaZ.AddAt(hisdz,id);
3429
3430       hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3431       hisP3->SetDirectory(0);
3432       fDeltaP3.AddAt(hisP3,id);
3433       //
3434       hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3435       hisP4->SetDirectory(0);
3436       fDeltaP4.AddAt(hisP4,id);
3437
3438       //
3439       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3440       hisdphi->SetDirectory(0);
3441       fDeltaPhi.AddAt(hisdphi,id);
3442       //
3443       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3444       hisdphiP->SetDirectory(0);
3445       fDeltaPhiP.AddAt(hisdphiP,id);
3446       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3447       hisSignal->SetDirectory(0);
3448       fSignals.AddAt(hisSignal,id);
3449     }
3450   }
3451
3452   SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3453   SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3454   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3455   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3456   //
3457   // Make THnSparse
3458   //
3459   //                   id   side   rod  bundle  beam  dP0  dP1  dP2  dP3  dP4 ncl dEdx
3460   Int_t binsLaser[12]= {336,  //id
3461                         2,    //side
3462                         6,    //rod
3463                         4,    //bundle
3464                         7,    //beam
3465                         300,  //dP0
3466                         300,  //dP1
3467                         300,  //dP2
3468                         300,  //dP3
3469                         300,  //dP4
3470                         80,   //ncl
3471                         50};   //dEdx
3472   Double_t xminLaser[12]= {0,  //id
3473                         0,     //side
3474                         0,     //rod
3475                         0,     //bundle
3476                         0,     //beam
3477                         -1,    //dP0
3478                         -1,    //dP1
3479                         -0.01,   //dP2
3480                         -0.01,  //dP3
3481                         -0.1,  //dP4
3482                         0,   //ncl
3483                         0};   //sqrt dEdx
3484   Double_t xmaxLaser[12]= {336,  //id
3485                         2,    //side
3486                         6,    //rod
3487                         4,    //bundle
3488                         7,    //beam
3489                         1,  //dP0
3490                         1,  //dP1
3491                         0.01,  //dP2
3492                         0.01,  //dP3
3493                         0.1,  //dP4
3494                         160,   //ncl
3495                         40};   //sqrt dEdx
3496
3497   TString nameLaser[12]= {"id",
3498                           "side",
3499                           "rod",
3500                           "bundle",
3501                           "beam",
3502                           "dP0",
3503                           "dP1",
3504                           "dP2",
3505                           "dP3",
3506                           "dP4",
3507                           "ncl",
3508                           "sqrt dEdx"};
3509   TString titleLaser[12]= {"id",
3510                           "side",
3511                           "rod",
3512                           "bundle",
3513                           "beam",
3514                           "#Delta_{P0}",
3515                           "#Delta_{P1}",
3516                           "#Delta_{P2}",
3517                           "#Delta_{P3}",
3518                           "#Delta_{P4}",
3519                           "N_{cl}",
3520                           "#sqrt{dEdx}"};
3521   fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3522   for (Int_t iaxis=1; iaxis<12; iaxis++){
3523     fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3524     fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3525   }
3526 }
3527
3528 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3529   //
3530   // Merge content of histograms
3531   //
3532   // Only first histogram is checked - all other should be the same
3533   if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3534
3535   if (!laser->fHisNclIn) return;  // empty histograms
3536   if (!fHisNclIn) MakeFitHistos();
3537   //
3538   fHisNclIn->Add(laser->fHisNclIn  );      //->Number of clusters inner
3539   fHisNclOut->Add(laser->fHisNclOut  );     //->Number of clusters outer
3540   fHisNclIO->Add(laser->fHisNclIO  );      //->Number of cluster inner outer
3541   fHisLclIn->Add(laser->fHisLclIn  );      //->Level arm inner
3542   fHisLclOut->Add(laser->fHisLclOut  );     //->Level arm outer
3543   fHisLclIO->Add(laser->fHisLclIO  );      //->Number of cluster inner outer
3544   fHisdEdx->Add(laser->fHisdEdx  );      //->dedx
3545   fHisdZfit->Add(laser->fHisdZfit  );      //->dedx
3546   //
3547   //
3548   fHisChi2YIn1->Add(laser->fHisChi2YIn1  );      //->chi2 y inner - line
3549   fHisChi2YOut1->Add(laser->fHisChi2YOut1  );     //->chi2 y inner - line
3550   fHisChi2YIn2->Add(laser->fHisChi2YIn2  );      //->chi2 y inner - parabola
3551   fHisChi2YOut2->Add(laser->fHisChi2YOut2  );     //->chi2 y inner - parabola
3552   fHisChi2YIO1->Add(laser->fHisChi2YIO1  );      //->chi2 y IO    - common
3553   fHisChi2ZIn1->Add(laser->fHisChi2ZIn1  );      //->chi2 z inner - line
3554   fHisChi2ZOut1->Add(laser->fHisChi2ZOut1  );     //->chi2 z inner - line
3555   fHisChi2ZIn2->Add(laser->fHisChi2ZIn2  );      //->chi2 z inner - parabola
3556   fHisChi2ZOut2->Add(laser->fHisChi2ZOut2  );     //->chi2 z inner - parabola
3557   fHisChi2ZIO1->Add(laser->fHisChi2ZIO1  );      //->chi2 z IO    - common
3558   //
3559   //
3560   fHisPy1vP0->Add(laser->fHisPy1vP0  );     //-> delta y   P0outer-P0inner - line
3561   fHisPy2vP0->Add(laser->fHisPy2vP0  );     //-> delta y   P0outer-P0inner - parabola
3562   fHisPy3vP0->Add(laser->fHisPy3vP0  );     //-> delta y   P0outer-P0inner - common parabola
3563   fHisPy1vP1->Add(laser->fHisPy1vP1  );     //-> delta ky  P1outer-P1inner - line
3564   fHisPy2vP1->Add(laser->fHisPy2vP1  );     //-> delta ky  P1outer-P1inner - parabola
3565   fHisPy3vP1->Add(laser->fHisPy3vP1  );     //-> delta ky  P1outer-P1inner - common parabola
3566   fHisPy2vP2In->Add(laser-> fHisPy2vP2In  );   //-> Curv  P2inner - parabola
3567   fHisPy2vP2Out->Add(laser->fHisPy2vP2Out  );  //-> Curv  P2outer - parabola
3568   fHisPy3vP2IO->Add(laser->fHisPy3vP2IO  );   //-> Curv  P2outerinner - common parabola
3569   //
3570   //
3571   fHisPz1vP0->Add(laser->fHisPz1vP0  );     //-> delta z   P0outer-P0inner - line
3572   fHisPz2vP0->Add(laser->fHisPz2vP0  );     //-> delta z   P0outer-P0inner - parabola
3573   fHisPz3vP0->Add(laser->fHisPz3vP0  );     //-> delta z   P0outer-P0inner - common parabola
3574   fHisPz1vP1->Add(laser->fHisPz1vP1  );     //-> delta kz  P1outer-P1inner - line
3575   fHisPz2vP1->Add(laser->fHisPz2vP1  );     //-> delta kz  P1outer-P1inner - parabola
3576   fHisPz3vP1->Add(laser->fHisPz3vP1  );     //-> delta kz  P1outer-P1inner - common parabola
3577   fHisPz2vP2In->Add(laser->fHisPz2vP2In  );   //-> Curv  P2inner - parabola
3578   fHisPz2vP2Out->Add(laser->fHisPz2vP2Out  );  //-> Curv  P2outer - parabola
3579   fHisPz3vP2IO->Add(laser->fHisPz3vP2IO  );   //-> Curv  P2outerinner - common parabola
3580   fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3581   fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3582   //
3583   //
3584   //
3585
3586
3587
3588
3589 }
3590
3591
3592
3593
3594 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3595   //
3596   // Dump fit information - collect information from the streamers
3597   //
3598   /*
3599     TChain * chainFit=0;
3600     TChain * chainTrack=0;
3601     TChain * chain=0;
3602     //
3603     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3604     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3605     AliXRDPROOFtoolkit tool;
3606     chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3607     chainTrack->Lookup();
3608     chainTrack->SetProof(kTRUE);
3609     chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3610     chainDrift->Lookup();
3611     chainDrift->SetProof(kTRUE);
3612
3613     chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3614     chain->Lookup();
3615     chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3616     chainFit->Lookup();
3617     chainFit->SetProof(kTRUE);
3618     chain->SetProof(kTRUE);
3619     AliTPCLaserTrack::LoadTracks();
3620     //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3621
3622   */
3623   //
3624   // Fit cuts
3625   //
3626   TCut cutP3("abs(Tr.fP[3])<0.1");
3627   TCut cutP4("abs(Tr.fP[4])<0.5");
3628   TCut cutPx = cutP3+cutP4;
3629   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3630   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3631   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3632   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3633   //
3634   TCut cutdEdx("sqrt(dEdx)>3");
3635   TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3636   TCut cutN("nclO>20&&nclI>20");
3637   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3638   //
3639   // Cluster cuts
3640   //
3641   TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3642   TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3643   TCut cutClX("abs(Cl[].fX)>10");
3644   TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3645   TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3646   TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3647   TCut cutQ("sqrt(Cl[].fMax)>4");
3648   TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3649
3650
3651   TH1F * phisAl     = 0;
3652   TH1F * phisAccept = 0;
3653   TH1F * phisOut    = 0;
3654   TProfile * pdEdx  = 0;
3655
3656   TProfile * pP0    = 0;
3657   TProfile * pP1    = 0;
3658   TProfile * pP2    = 0;
3659   TProfile * pP3    = 0;
3660   TProfile * pP4    = 0;
3661   //
3662   TProfile * pNclI  = 0;
3663   TProfile * pNclO  = 0;
3664   //
3665   TProfile * pchi2YIn   =0;
3666   TProfile * pchi2ZIn   =0;
3667   TProfile * pchi2YOut  =0;
3668   TProfile * pchi2ZOut  =0;
3669   TProfile * pchi2YInOut =0;
3670   TProfile * pchi2ZInOut =0;;
3671   // laser counters
3672   chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3673   phisAl = (TH1F*)gROOT->FindObject("hisAl");
3674   chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3675   phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3676   chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3677   phisOut = (TH1F*)gROOT->FindObject("hisOut");
3678   //
3679   chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3680   pdEdx   = (TProfile*)gROOT->FindObject("hdEdx");
3681   // track param
3682   //
3683   chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3684   pP0   = (TProfile*)gROOT->FindObject("hP0");
3685   chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3686   pP1   = (TProfile*)gROOT->FindObject("hP1");
3687   chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3688   pP2   = (TProfile*)gROOT->FindObject("hP2");
3689   chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3690   pP3   = (TProfile*)gROOT->FindObject("hP3");
3691   chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3692   pP4   = (TProfile*)gROOT->FindObject("hP4");
3693   //
3694   chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3695   pNclI   = (TProfile*)gROOT->FindObject("hNclI");
3696   chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3697   pNclO   = (TProfile*)gROOT->FindObject("hNclO");
3698   //
3699   //
3700   chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3701   pchi2YIn   = (TProfile*)gROOT->FindObject("hChi2YIn");
3702   chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3703   pchi2YOut   = (TProfile*)gROOT->FindObject("hChi2YOut");
3704   chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3705   pchi2YInOut   = (TProfile*)gROOT->FindObject("hChi2YInOut");
3706   chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3707   pchi2ZIn   = (TProfile*)gROOT->FindObject("hChi2ZIn");
3708   chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3709   pchi2ZOut   = (TProfile*)gROOT->FindObject("hChi2ZOut");
3710   chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3711   pchi2ZInOut   = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3712   //
3713   // second derivatives
3714   //
3715   TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3716   chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3717   TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3718   chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3719   TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3720   chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3721   //
3722   phisPy2In->FitSlicesY();
3723   TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3724   TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3725   TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3726   //
3727   phisPy2Out->FitSlicesY();
3728   TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3729   TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3730   TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3731   //
3732   phisPy2InOut->FitSlicesY();
3733   TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3734   TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3735   TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3736   //
3737   TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3738   chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3739   TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3740   chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3741   TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3742   chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3743   //
3744   phisPz2In->FitSlicesY();
3745   TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3746   TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3747   TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3748   //
3749   phisPz2Out->FitSlicesY();
3750   TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3751   TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3752   TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3753   //
3754   phisPz2InOut->FitSlicesY();
3755   TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3756   TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3757   TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3758   //
3759   //
3760   //
3761
3762
3763   {
3764     TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3765     for (Int_t ilaser=0; ilaser<336; ilaser++){
3766       Float_t all=phisAl->GetBinContent(ilaser+1);
3767       Float_t accept=phisAccept->GetBinContent(ilaser+1);
3768       Float_t out=phisOut->GetBinContent(ilaser+1);
3769       Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3770       Float_t mP0 = pP0->GetBinContent(ilaser+1);
3771       Float_t mP1 = pP1->GetBinContent(ilaser+1);
3772       Float_t mP2 = pP2->GetBinContent(ilaser+1);
3773       Float_t mP3 = pP3->GetBinContent(ilaser+1);
3774       Float_t mP4 = pP4->GetBinContent(ilaser+1);
3775
3776
3777       Float_t nclI  = pNclI->GetBinContent(ilaser+1);
3778       Float_t nclO  = pNclO->GetBinContent(ilaser+1);
3779       //
3780       Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3781       Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3782       Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3783       Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3784       Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3785       Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3786       //
3787       Float_t entriesPy2In  = phisPy2InEntries->GetBinContent(ilaser+1);
3788       Float_t meanPy2In     = phisPy2InMean->GetBinContent(ilaser+1);
3789       Float_t sigmaPy2In    = phisPy2InSigma->GetBinContent(ilaser+1);
3790       //
3791       Float_t entriesPy2Out  = phisPy2OutEntries->GetBinContent(ilaser+1);
3792       Float_t meanPy2Out     = phisPy2OutMean->GetBinContent(ilaser+1);
3793       Float_t sigmaPy2Out    = phisPy2OutSigma->GetBinContent(ilaser+1);
3794       //
3795       Float_t entriesPy2InOut  = phisPy2InOutEntries->GetBinContent(ilaser+1);
3796       Float_t meanPy2InOut     = phisPy2InOutMean->GetBinContent(ilaser+1);
3797       Float_t sigmaPy2InOut    = phisPy2InOutSigma->GetBinContent(ilaser+1);
3798       //
3799       Float_t entriesPz2In  = phisPz2InEntries->GetBinContent(ilaser+1);
3800       Float_t meanPz2In     = phisPz2InMean->GetBinContent(ilaser+1);
3801       Float_t sigmaPz2In    = phisPz2InSigma->GetBinContent(ilaser+1);
3802       //
3803       Float_t entriesPz2Out  = phisPz2OutEntries->GetBinContent(ilaser+1);
3804       Float_t meanPz2Out     = phisPz2OutMean->GetBinContent(ilaser+1);
3805       Float_t sigmaPz2Out    = phisPz2OutSigma->GetBinContent(ilaser+1);
3806       //
3807       Float_t entriesPz2InOut  = phisPz2InOutEntries->GetBinContent(ilaser+1);
3808       Float_t meanPz2InOut     = phisPz2InOutMean->GetBinContent(ilaser+1);
3809       Float_t sigmaPz2InOut    = phisPz2InOutSigma->GetBinContent(ilaser+1);
3810
3811       AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3812       (*pcstream)<<"Scan"<<
3813         "Run="<<id<<
3814         "LTr.="<<ltrp<<
3815         "all="<<all<<
3816         "accept="<<accept<<
3817         "out="<<out<<
3818         "sdedx="<<sdedx<<
3819         "mP0="<<mP0<<
3820         "mP1="<<mP1<<
3821         "mP2="<<mP2<<
3822         "mP3="<<mP3<<
3823         "mP4="<<mP4<<
3824         "nclI="<<nclI<<
3825         "nclO="<<nclO<<
3826         "chi2YIn="<<chi2YIn<<
3827         "chi2YOut="<<chi2YOut<<
3828         "chi2YInOut="<<chi2YInOut<<
3829         "chi2ZIn="<<chi2ZIn<<
3830         "chi2ZOut="<<chi2ZOut<<
3831         "chi2ZInOut="<<chi2ZInOut<<
3832         //
3833         "nPy2In="<<entriesPy2In<<
3834         "mPy2In="<<meanPy2In<<
3835         "sPy2In="<<sigmaPy2In<<
3836         //
3837         "nPy2Out="<<entriesPy2Out<<
3838         "mPy2Out="<<meanPy2Out<<
3839         "sPy2Out="<<sigmaPy2Out<<
3840         //
3841         "nPy2InOut="<<entriesPy2InOut<<
3842         "mPy2InOut="<<meanPy2InOut<<
3843         "sPy2InOut="<<sigmaPy2InOut<<
3844         //
3845         "nPz2In="<<entriesPz2In<<
3846         "mPz2In="<<meanPz2In<<
3847         "sPz2In="<<sigmaPz2In<<
3848         //
3849         "nPz2Out="<<entriesPz2Out<<
3850         "mPz2Out="<<meanPz2Out<<
3851         "sPz2Out="<<sigmaPz2Out<<
3852         //
3853         "nPz2InOut="<<entriesPz2InOut<<
3854         "mPz2InOut="<<meanPz2InOut<<
3855         "sPz2InOut="<<sigmaPz2InOut<<
3856         "\n";
3857     }
3858
3859     delete pcstream;
3860   }
3861 }
3862
3863 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3864                                          TVectorD& meanSlope,
3865                                          TVectorD& sectorArray,
3866                                          Int_t option)
3867 {
3868   // This method should ideally go in AliTPCLaser
3869   // option == 0 (pads outer - closest to beam)
3870   // option == 1 (pads inner)
3871   // option == 2 (time outer)
3872   // option == 3 (time inner)
3873   Int_t nFailures = 0;
3874
3875   for(Int_t id = 0; id < 336; id++) {
3876
3877     if (!AliTPCLaserTrack::GetTracks())
3878       AliTPCLaserTrack::LoadTracks();
3879     AliTPCLaserTrack *ltrp =
3880       (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3881
3882     AliExternalTrackParam trackParam(*ltrp);
3883
3884     Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3885     if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3886       deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3887
3888     Double_t angle = trackParam.GetAlpha();
3889     if(angle<0)
3890       angle += 2*TMath::Pi();
3891     if(trackParam.GetSnp()>0) // track points to sector "before"
3892       angle -= deltaangle*TMath::DegToRad();
3893     else // track points to sector "after"
3894       angle += deltaangle*TMath::DegToRad();
3895
3896     Bool_t success = trackParam.Rotate(angle);
3897
3898     if(!success) {
3899       //      cout << "WARNING: Rotate failed for ID: " << id << endl;
3900       nFailures++;
3901     }
3902
3903     angle *= TMath::RadToDeg();
3904
3905     Int_t sector = TMath::Nint((angle-10.0)/20.0);
3906     if(sector<0)
3907       sector += 18;
3908     else if(sector>=18)
3909       sector -= 18;
3910     if(ltrp->GetSide()==1) // C side
3911       sector += 18;
3912     if(option==0 || option==2)
3913       sector += 36;
3914     if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
3915       sector += 36;
3916
3917     sectorArray[id] = sector;
3918
3919     const Double_t x0 = 0;
3920
3921     Double_t slopey  = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3922     Double_t slopez  = trackParam.GetTgl();
3923     //  One needs a factor sqrt(1+slopey**2) to take into account the
3924     //  longer path length
3925     slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3926     if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3927       slopez *= -1;
3928     //    Double_t offsetz = trackParam.GetZ();
3929     Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
3930     Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
3931     if(option==2 || option==3) {
3932       meanOffset[id] = offsetz; meanSlope[id] = slopez;
3933     } else {
3934       meanOffset[id] = offsety; meanSlope[id] = slopey;
3935     }
3936   }
3937
3938   if(nFailures>0)
3939     AliWarning(Form("Rotate method failed %d times", nFailures));
3940 }
3941
3942
3943   /*
3944     TFile f("vscan.root");
3945    */
3946
3947   /*
3948     pad binning effect
3949    chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
3950    //
3951    chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
3952    //
3953
3954 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);
3955
3956
3957 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)
3958
3959   */
3960
3961
3962
3963
3964
3965   /*
3966   // check edge effects
3967   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
3968   //
3969   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
3970
3971   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)
3972
3973
3974
3975   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
3976   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
3977
3978 */
3979
3980
3981
3982
3983
3984   /*
3985     Edge y effect
3986
3987     dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
3988
3989
3990     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)
3991
3992     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)
3993
3994
3995
3996
3997
3998     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)
3999
4000     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)
4001
4002
4003
4004     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)
4005
4006     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)
4007
4008
4009
4010 */
4011
4012
4013 /*
4014
4015 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4016
4017 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4018
4019
4020
4021 chainFit->Draw("LTr.fId","nclI>10",100000)
4022
4023 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4024
4025 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4026
4027 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4028
4029 */
4030
4031
4032
4033
4034
4035
4036 /*
4037  gSystem->Load("libSTAT.so")
4038  TStatToolkit toolkit;
4039  Double_t chi2;
4040  TVectorD fitParam;
4041  TMatrixD covMatrix;
4042  Int_t npoints;
4043
4044  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
4045
4046
4047 TString fstring="";
4048 //
4049 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
4050 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
4051 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
4052 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4
4053 //
4054 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
4055 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
4056 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
4057 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
4058 //
4059 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
4060 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
4061 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
4062 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
4063
4064
4065
4066
4067  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4068
4069  treeT->SetAlias("fit",strq0->Data());
4070
4071
4072  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4073
4074  treeT->SetAlias("fitP",strqP->Data());
4075
4076
4077  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4078  treeT->SetAlias("fitD",strqDrift->Data());
4079
4080
4081 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4082 {
4083 for (Int_t i=0; i<6;i++){
4084 treeT->SetLineColor(i+2);
4085 treeT->SetMarkerSize(1);
4086 treeT->SetMarkerStyle(22+i);
4087 treeT->SetMarkerColor(i+2);
4088
4089 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4090 }
4091 }
4092  */
4093
4094
4095
4096 /*
4097   TTree * tree = (TTree*)f.Get("FitModels");
4098
4099   TEventList listLFit0("listLFit0","listLFit0");
4100   TEventList listLFit1("listLFit1","listLFit1");
4101   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4102   tree->SetEventList(&listLFit0);
4103
4104
4105
4106
4107   gSystem->Load("libSTAT.so")
4108   TStatToolkit toolkit;
4109   Double_t chi2;
4110   TVectorD fitParam;
4111   TMatrixD covMatrix;
4112   Int_t npoints;
4113
4114   chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4115   chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4116
4117
4118   TString fstring="";
4119   fstring+="cos(dp)++";
4120   fstring+="sin(dp)++";
4121   fstring+="cos(dt)++";
4122   fstring+="sin(dt)++";
4123
4124   TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4125
4126
4127
4128 */
4129
4130
4131
4132 /*
4133   Edge effects
4134   //
4135   //
4136   //
4137   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4138   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4139   AliXRDPROOFtoolkit tool;
4140   TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4141   chainTrack->Lookup();
4142   chainTrack->SetProof(kTRUE);
4143
4144   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4145   chain->Lookup();
4146   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4147   chainFit->Lookup();
4148   chainFit->SetProof(kTRUE);
4149   chain->SetProof(kTRUE);
4150   //
4151   // Fit cuts
4152   //
4153   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4154   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4155   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4156   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4157   //
4158   TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4159   TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4160   TCut cutN("nclO>20&&nclI>20");
4161   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4162   //
4163   // Cluster cuts
4164   //
4165   TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4166   TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4167   TCut cutClX("abs(Cl.fX)>10");
4168   TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4169   TCut cutCl=cutClY+cutClZ+cutClX;
4170
4171
4172   // check edge effects
4173   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4174   //
4175   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4176
4177   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)
4178
4179
4180
4181   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4182   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
4183
4184 */
4185