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