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