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