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