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