TPC/AliTPCcalibLaser.cxx

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