TPC/AliTPCcalibLaser.cxx

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