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     Double_t meanp3 = hisP3->GetMean();
2345     Double_t rmsp3  = hisP3->GetRMS();
2346     //
2347     if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
2348       hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
2349     Double_t gp41 = fg.GetParameter(1);
2350     Double_t gp42 = fg.GetParameter(2);
2351     Double_t meanp4 = hisP4->GetMean();
2352     Double_t rmsp4  = hisP4->GetRMS();
2353     //
2354     Float_t meanS=hisS->GetMean();
2355     //
2356     Double_t lxyz[3];
2357     Double_t lpxyz[3];
2358     ltrp->GetXYZ(lxyz);
2359     ltrp->GetPxPyPz(lpxyz);
2360
2361     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
2362     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
2363     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
2364     if (gphiP2>kcutP0) isOK=kFALSE;
2365     //
2366     //
2367     //
2368     TH1 *his =0;
2369     //
2370     his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
2371     Float_t mnclIn =  his->GetMean();
2372     delete his;
2373     his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
2374     Float_t mnclOut =  his->GetMean();
2375     delete his;
2376     his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
2377     Float_t mnclIO =  his->GetMean();
2378     delete his;
2379     his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
2380     Float_t mLclIn =  his->GetMean();
2381     delete his;
2382     his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
2383     Float_t mLclOut =  his->GetMean();
2384     delete his;
2385     his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
2386     Float_t mLclIO =  his->GetMean();
2387     delete his;
2388     //
2389     his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
2390     Float_t mdEdx =  his->GetMean();
2391     delete his;
2392     //
2393     //
2394     //
2395     //
2396     his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - line
2397     Float_t  mChi2YIn1=  his->GetMean();
2398     delete his;
2399     his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - line
2400     Float_t mChi2YOut1 =  his->GetMean();
2401     delete his;
2402     his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - parabola
2403     Float_t mChi2YIn2 =  his->GetMean();
2404     delete his;
2405     his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - parabola
2406     Float_t mChi2YOut2 =  his->GetMean();
2407     delete his;
2408     his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 y IO    - common
2409     Float_t mChi2YIO1 =  his->GetMean();
2410     delete his;
2411     his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - line
2412     Float_t mChi2ZIn1 =  his->GetMean();
2413     delete his;
2414     his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - line
2415     Float_t mChi2ZOut1 =  his->GetMean();
2416     delete his;
2417     his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - parabola
2418     Float_t mChi2ZIn2 =  his->GetMean();
2419     delete his;
2420     his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - parabola
2421     Float_t mChi2ZOut2 =  his->GetMean();
2422     delete his;
2423     his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 z IO    - common
2424     Float_t mChi2ZIO1 =  his->GetMean();
2425     delete his;
2426     //
2427     // fit res. histos
2428     //
2429     his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
2430     Float_t  edZfit  = his->GetEntries();
2431     Float_t  mdZfit =  his->GetMean();
2432     Float_t  rdZfit  = his->GetRMS();
2433     delete his;
2434
2435     his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - line
2436     Float_t  ePy1vP0  = his->GetEntries();
2437     Float_t  mPy1vP0  = his->GetMean();
2438     Float_t  rPy1vP0  = his->GetRMS();
2439     delete his;
2440
2441     his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - parabola
2442     Float_t  ePy2vP0  = his->GetEntries();
2443     Float_t  mPy2vP0  = his->GetMean();
2444     Float_t  rPy2vP0  = his->GetRMS();
2445     delete his;
2446
2447     his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - common parabola
2448     Float_t  ePy3vP0  = his->GetEntries();
2449     Float_t  mPy3vP0  = his->GetMean();
2450     Float_t  rPy3vP0  = his->GetRMS();
2451     delete his;
2452
2453     his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - line
2454     Float_t  ePy1vP1  = his->GetEntries();
2455     Float_t  mPy1vP1  = his->GetMean();
2456     Float_t  rPy1vP1  = his->GetRMS();
2457     delete his;
2458
2459     his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - parabola
2460     Float_t  ePy2vP1  = his->GetEntries();
2461     Float_t  mPy2vP1  = his->GetMean();
2462     Float_t  rPy2vP1  = his->GetRMS();
2463     delete his;
2464
2465     his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - common parabola
2466     Float_t  ePy3vP1  = his->GetEntries();
2467     Float_t  mPy3vP1  = his->GetMean();
2468     Float_t  rPy3vP1  = his->GetRMS();
2469     delete his;
2470
2471     his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2472     Float_t  ePy2vP2In  = his->GetEntries();
2473     Float_t  mPy2vP2In  = his->GetMean();
2474     Float_t  rPy2vP2In  = his->GetRMS();
2475     delete his;
2476
2477     his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2478     Float_t  ePy2vP2Out  = his->GetEntries();
2479     Float_t  mPy2vP2Out  = his->GetMean();
2480     Float_t  rPy2vP2Out  = his->GetRMS();
2481     delete his;
2482
2483     his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2484     Float_t  ePy3vP2IO  = his->GetEntries();
2485     Float_t  mPy3vP2IO  = his->GetMean();
2486     Float_t  rPy3vP2IO  = his->GetRMS();
2487     delete his;
2488
2489     //
2490     //
2491     his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - line
2492     Float_t  ePz1vP0  = his->GetEntries();
2493     Float_t  mPz1vP0  = his->GetMean();
2494     Float_t  rPz1vP0  = his->GetRMS();
2495     delete his;
2496
2497     his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - parabola
2498     Float_t  ePz2vP0  = his->GetEntries();
2499     Float_t  mPz2vP0  = his->GetMean();
2500     Float_t  rPz2vP0  = his->GetRMS();
2501     delete his;
2502
2503     his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - common parabola
2504     Float_t  ePz3vP0  = his->GetEntries();
2505     Float_t  mPz3vP0  = his->GetMean();
2506     Float_t  rPz3vP0  = his->GetRMS();
2507     delete his;
2508
2509     his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - line
2510     Float_t  ePz1vP1  = his->GetEntries();
2511     Float_t  mPz1vP1  = his->GetMean();
2512     Float_t  rPz1vP1  = his->GetRMS();
2513     delete his;
2514
2515     his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - parabola
2516     Float_t  ePz2vP1  = his->GetEntries();
2517     Float_t  mPz2vP1  = his->GetMean();
2518     Float_t  rPz2vP1  = his->GetRMS();
2519     delete his;
2520
2521     his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - common parabola
2522     Float_t  ePz3vP1  = his->GetEntries();
2523     Float_t  mPz3vP1  = his->GetMean();
2524     Float_t  rPz3vP1  = his->GetRMS();
2525     delete his;
2526
2527     his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2528     Float_t  ePz2vP2In  = his->GetEntries();
2529     Float_t  mPz2vP2In  = his->GetMean();
2530     Float_t  rPz2vP2In  = his->GetRMS();
2531     delete his;
2532
2533     his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2534     Float_t  ePz2vP2Out  = his->GetEntries();
2535     Float_t  mPz2vP2Out  = his->GetMean();
2536     Float_t  rPz2vP2Out  = his->GetRMS();
2537     delete his;
2538
2539     his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2540     Float_t  ePz3vP2IO  = his->GetEntries();
2541     Float_t  mPz3vP2IO  = his->GetMean();
2542     Float_t  rPz3vP2IO  = his->GetRMS();
2543     delete his;
2544
2545     // Fit absolute laser residuals
2546     TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
2547     TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
2548
2549     Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
2550     Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
2551
2552     TVectorD vecX(159);        // X
2553     TVectorD vecY(159);        // Y
2554     TVectorD vecR(159);        // R
2555     TVectorD vecDY(159);       // absolute residuals in Y
2556     TVectorD vecDZ(159);       // absolute residuals in Z
2557     TVectorD vecN(159);        // number of clusters
2558     TVectorD vecEy(159);       //error y
2559     TVectorD vecEz(159);       //error z
2560     TVectorD vecPhi(159);      // local tangent
2561     TVectorD vecPhiR(159);     // local tangent
2562     // magnetic field integrals
2563     TVectorD vecIBR(159);        // radial
2564     TVectorD vecIBRPhi(159);     // r-phi
2565     TVectorD vecIBLX(159);       // local x
2566     TVectorD vecIBLY(159);       // local y
2567     TVectorD vecIBGX(159);       // local x
2568     TVectorD vecIBGY(159);       // local y
2569     TVectorD vecIBZ(159);        // z
2570     //
2571     for (Int_t irow=0;irow<159;irow++){
2572       vecIBR[irow]=0;
2573       vecIBRPhi[irow]=0;
2574       vecIBLX[irow]=0;
2575       vecIBLY[irow]=0;
2576       vecIBGX[irow]=0;
2577       vecIBGY[irow]=0;
2578       vecIBZ[irow]=0;
2579       Double_t gx    =(*(ltrp->fVecGX))[irow];
2580       Double_t gy    =(*(ltrp->fVecGY))[irow];
2581       Int_t    lsec  =TMath::Nint((*(ltrp->fVecSec))[irow]);
2582       Double_t   ca  =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
2583       Double_t   sa  =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
2584       xyz[2]=(*(ltrp->fVecGZ))[irow];
2585       xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
2586       xyz[1]=TMath::ATan2(gy,gx);
2587       Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
2588       if (magF){
2589         magF->GetTPCIntCyl(xyz,bxyz);
2590         magF->GetTPCInt(gxyz,bgxyz);
2591         vecIBR[irow]=bxyz[0];
2592         vecIBRPhi[irow]=bxyz[1];
2593         //
2594         vecIBGX[irow]=bgxyz[0];
2595         vecIBGY[irow]=bgxyz[1];
2596         //
2597         vecIBLX[irow]=  bgxyz[0]*ca+bgxyz[1]*sa;
2598         vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
2599         //
2600
2601         vecIBZ[irow]=bxyz[2];
2602       }
2603     }
2604
2605
2606     lfabsyInner.ClearPoints();
2607     lfabszInner.ClearPoints();
2608     lfabsyOuter.ClearPoints();
2609     lfabszOuter.ClearPoints();
2610     // dummy fit values
2611     Int_t    nClInY         = 0;
2612     Float_t  yAbsInOffset  = -100;
2613     Float_t  yAbsInSlope   = -100;
2614     Float_t  yAbsInDeltay  = -100;
2615     Int_t    nClInZ         = 0;
2616     Float_t  zAbsInOffset  = -100;
2617     Float_t  zAbsInSlope   = -100;
2618     Float_t  zAbsInDeltaz  = -100;
2619     Int_t    nClOutY         = 0;
2620     Float_t  yAbsOutOffset  = -100;
2621     Float_t  yAbsOutSlope   = -100;
2622     Float_t  yAbsOutDeltay  = -100;
2623     Int_t    nClOutZ         = 0;
2624     Float_t  zAbsOutOffset  = -100;
2625     Float_t  zAbsOutSlope   = -100;
2626     Float_t  zAbsOutDeltaz  = -100;
2627
2628     Float_t  lasTanPhiLocIn = -100;
2629     Float_t  lasTanPhiLocOut = -100;
2630
2631     if(histAbsY->GetEntries()>0) {
2632
2633       Double_t rotAngOut = 10;
2634       Double_t rotAngIn = 10;
2635       if((secInner%36)!=(secOuter%36))
2636         rotAngIn += 20; // 30 degrees
2637
2638       // Calculate laser mirror X position in local frame
2639       Double_t laserposOut =
2640         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
2641       Double_t laserposIn =
2642         TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
2643
2644       // Calculate laser tan phi in local frame
2645       lasTanPhiLocOut  = TMath::ASin(ltrp->GetSnp());
2646       if(lasTanPhiLocOut<0) {
2647         lasTanPhiLocIn  = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
2648         lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
2649       } else {
2650
2651         lasTanPhiLocIn  = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
2652         lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
2653       }
2654
2655       lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
2656       lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
2657
2658       TProfile* yprof  = histAbsY->ProfileX("yprof");
2659       TProfile* zprof  = histAbsZ->ProfileX("zprof");
2660
2661       for(Int_t bin = 1; bin<=159; bin++) {
2662
2663         if(yprof->GetBinEntries(bin)<10&&
2664            zprof->GetBinEntries(bin)<10) {
2665           continue;
2666         }
2667
2668         // There is a problem in defining inner and outer sectors for
2669         // the outer beams (0 and 6) where both sectors are OROCs. To
2670         // make sure there is no overlap row 94 to 99 are cutted.
2671         if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
2672           continue;
2673
2674         Int_t row = (bin-1);
2675         if(row>62)
2676           row -= 63;
2677
2678         Bool_t isOuter = kTRUE;
2679         Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
2680
2681         if(bin<=62 ||
2682            (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
2683
2684           isOuter = kFALSE;
2685           sector = TMath::Nint(fBeamSectorInner[id]);
2686         }
2687
2688
2689         Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
2690         vecN[bin-1] =yprof->GetBinEntries(bin);
2691         vecEy[bin-1]=yprof->GetBinError(bin);
2692         vecEz[bin-1]=zprof->GetBinError(bin);
2693         vecX[bin-1] = x;
2694         vecDY[bin-1] = yprof->GetBinContent(bin);
2695         vecDZ[bin-1] = zprof->GetBinContent(bin);
2696         if(!isOuter) { // inner
2697           vecPhi[bin-1]=lasTanPhiLocIn;
2698           // Calculate local y from residual and database
2699           Double_t y  = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
2700             + vecDY[bin-1];
2701           vecY[bin-1] = y;
2702           Double_t r  = TMath::Sqrt(x*x + y*y);
2703           vecR[bin-1] = r;
2704           // Find angle between laser vector and R vector
2705           // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
2706           Double_t cosPhi = x + y*fBeamSlopeYInner[id];
2707           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
2708           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2709           if(lasTanPhiLocIn<0)
2710             vecPhiR[bin-1]*=-1; // to have the same sign
2711
2712           if(yprof->GetBinEntries(bin)>=10) {
2713             lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
2714                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2715           }
2716           if(zprof->GetBinEntries(bin)>=10) {
2717             lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
2718                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2719           }
2720         } else { // outer
2721           vecPhi[bin-1]=lasTanPhiLocOut;
2722           // Calculate local y from residual and database
2723           Double_t y  = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
2724             + vecDY[bin-1];
2725           vecY[bin-1] = y;
2726           Double_t r  = TMath::Sqrt(x*x + y*y);
2727           vecR[bin-1] = r;
2728
2729           Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
2730           cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
2731           vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
2732           if(lasTanPhiLocOut<0)
2733             vecPhiR[bin-1]*=-1; // to have the same sign
2734
2735           if(yprof->GetBinEntries(bin)>=10) {
2736             lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
2737                                   TMath::Max(yprof->GetBinError(bin), 0.001));
2738           }
2739           if(zprof->GetBinEntries(bin)>=10) {
2740             lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
2741                                   TMath::Max(zprof->GetBinError(bin), 0.001));
2742           }
2743         }
2744         // global position
2745
2746       }
2747
2748       delete yprof; delete zprof;
2749
2750
2751       // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
2752       nClInY = lfabsyInner.GetNpoints();
2753       if(lfabsyInner.GetNpoints()>10) {
2754         lfabsyInner.EvalRobust(0.95);
2755
2756         TVectorD result(2);
2757         lfabsyInner.GetParameters(result);
2758         yAbsInOffset = result[0];
2759         yAbsInSlope  = result[1];
2760         yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
2761       }
2762       nClInZ = lfabszInner.GetNpoints();
2763       if(lfabszInner.GetNpoints()>10) {
2764         lfabszInner.EvalRobust(0.95);
2765
2766         TVectorD result(2);
2767         lfabszInner.GetParameters(result);
2768         zAbsInOffset = result[0];
2769         zAbsInSlope  = result[1];
2770         zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
2771       }
2772       nClOutY = lfabsyOuter.GetNpoints();
2773       if(lfabsyOuter.GetNpoints()>10) {
2774         lfabsyOuter.EvalRobust(0.95);
2775
2776         TVectorD result(2);
2777         lfabsyOuter.GetParameters(result);
2778         yAbsOutOffset = result[0];
2779         yAbsOutSlope  = result[1];
2780         yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
2781       }
2782       nClOutZ = lfabszOuter.GetNpoints();
2783       if(lfabszOuter.GetNpoints()>10) {
2784         lfabszOuter.EvalRobust(0.95);
2785
2786         TVectorD result(2);
2787         lfabszOuter.GetParameters(result);
2788         zAbsOutOffset = result[0];
2789         zAbsOutSlope  = result[1];
2790         zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
2791       }
2792     }
2793
2794
2795     Int_t itime=-1;
2796     Float_t  coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
2797     Float_t  coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
2798     Float_t  coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
2799     Float_t  coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
2800     //
2801     Float_t  skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
2802     Float_t  skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
2803     //
2804     Float_t  ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
2805     Float_t  ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
2806
2807     //
2808     (*pcstream)<<"Mean"<<
2809       "run="<<run<<               //
2810       //voltages
2811       "UcIA="<<coverIA<<
2812       "UcIC="<<coverIC<<
2813       "UcOA="<<coverOA<<
2814       "UcOC="<<coverOC<<
2815       "UsA="<<skirtA<<
2816       "UsC="<<skirtC<<
2817       "UggA="<<ggOffA<<
2818       "UggC="<<ggOffC<<
2819       //
2820       "isOK="<<isOK<<             //
2821       "id="<<id<<                 // track id
2822       "entries="<<entries<<       // number of entries
2823       "bz="<<bfield<<             // bfield
2824       "LTr.="<<ltrp<<             // refernece track
2825       //
2826       "mnclIn="<<mnclIn<<         // mean number of clusters in inner
2827       "mnclOut="<<mnclOut<<       // mean number of clusters in outer
2828       "mnclIO="<<mnclIO<<         // mean number of clusters in inner+outer
2829       "mLclIn="<<mLclIn<<         // mean number of clusters in inner
2830       "mLclOut="<<mLclOut<<       // mean number of clusters in outer
2831       "mLclIO="<<mLclIO<<         // mean number of clusters in inner+outer
2832       "mdEdx="<<mdEdx<<           // mean dEdx
2833       "edZfit="<<edZfit<<           // entries z fit
2834       "mdZfit="<<mdZfit<<           // mean z fit
2835       "rdZfit="<<rdZfit<<           // RMS z fit
2836       //
2837       //
2838       "mChi2YIn1="<<mChi2YIn1<<       //->chi2 y inner - line
2839       "mChi2YOut1="<<mChi2YOut1<<     //->chi2 y inner - line
2840       "mChi2YIn2="<<mChi2YIn2<<       //->chi2 y inner - parabola
2841       "mChi2YOut2="<<mChi2YOut2<<     //->chi2 y inner - parabola
2842       "mChi2YIO1="<<mChi2YIO1<<       //->chi2 y IO    - common
2843       "mChi2ZIn1="<<mChi2ZIn1<<       //->chi2 z inner - line
2844       "mChi2ZOut1="<<mChi2ZOut1<<     //->chi2 z inner - line
2845       "mChi2ZIn2="<<mChi2ZIn2<<       //->chi2 z inner - parabola
2846       "mChi2ZOut2="<<mChi2ZOut2<<     //->chi2 z inner - parabola
2847       "mChi2ZIO1="<<mChi2ZIO1<<       //->chi2 z IO    - common
2848       //
2849       //
2850       "ePy1vP0="<<ePy1vP0<<
2851       "mPy1vP0="<<mPy1vP0<<
2852       "rPy1vP0="<<rPy1vP0<<
2853       "ePy2vP0="<<ePy2vP0<<
2854       "mPy2vP0="<<mPy2vP0<<
2855       "rPy2vP0="<<rPy2vP0<<
2856       "ePy3vP0="<<ePy3vP0<<
2857       "mPy3vP0="<<mPy3vP0<<
2858       "rPy3vP0="<<rPy3vP0<<
2859       "ePy1vP1="<<ePy1vP1<<
2860       "mPy1vP1="<<mPy1vP1<<
2861       "rPy1vP1="<<rPy1vP1<<
2862       "ePy2vP1="<<ePy2vP1<<
2863       "mPy2vP1="<<mPy2vP1<<
2864       "rPy2vP1="<<rPy2vP1<<
2865       "ePy3vP1="<<ePy3vP1<<
2866       "mPy3vP1="<<mPy3vP1<<
2867       "rPy3vP1="<<rPy3vP1<<
2868       "ePy2vP2In="<<ePy2vP2In<<
2869       "mPy2vP2In="<<mPy2vP2In<<
2870       "rPy2vP2In="<<rPy2vP2In<<
2871       "ePy2vP2Out="<<ePy2vP2Out<<
2872       "mPy2vP2Out="<<mPy2vP2Out<<
2873       "rPy2vP2Out="<<rPy2vP2Out<<
2874       "ePy3vP2IO="<<ePy3vP2IO<<
2875       "mPy3vP2IO="<<mPy3vP2IO<<
2876       "rPy3vP2IO="<<rPy3vP2IO<<
2877       //
2878       //
2879       "ePz1vP0="<<ePz1vP0<<
2880       "mPz1vP0="<<mPz1vP0<<
2881       "rPz1vP0="<<rPz1vP0<<
2882       "ePz2vP0="<<ePz2vP0<<
2883       "mPz2vP0="<<mPz2vP0<<
2884       "rPz2vP0="<<rPz2vP0<<
2885       "ePz3vP0="<<ePz3vP0<<
2886       "mPz3vP0="<<mPz3vP0<<
2887       "rPz3vP0="<<rPz3vP0<<
2888       "ePz1vP1="<<ePz1vP1<<
2889       "mPz1vP1="<<mPz1vP1<<
2890       "rPz1vP1="<<rPz1vP1<<
2891       "ePz2vP1="<<ePz2vP1<<
2892       "mPz2vP1="<<mPz2vP1<<
2893       "rPz2vP1="<<rPz2vP1<<
2894       "ePz3vP1="<<ePz3vP1<<
2895       "mPz3vP1="<<mPz3vP1<<
2896       "rPz3vP1="<<rPz3vP1<<
2897       "ePz2vP2In="<<ePz2vP2In<<
2898       "mPz2vP2In="<<mPz2vP2In<<
2899       "rPz2vP2In="<<rPz2vP2In<<
2900       "ePz2vP2Out="<<ePz2vP2Out<<
2901       "mPz2vP2Out="<<mPz2vP2Out<<
2902       "rPz2vP2Out="<<rPz2vP2Out<<
2903       "ePz3vP2IO="<<ePz3vP2IO<<
2904       "mPz3vP2IO="<<mPz3vP2IO<<
2905       "rPz3vP2IO="<<rPz3vP2IO<<
2906       //
2907       //
2908       //
2909       "lx0="<<lxyz[0]<<            // reference x
2910       "lx1="<<lxyz[1]<<            // reference y
2911       "lx2="<<lxyz[2]<<            // refernece z
2912       "lpx0="<<lpxyz[0]<<           // reference x
2913       "lpx1="<<lpxyz[1]<<          // reference y
2914       "lpx2="<<lpxyz[2]<<          // refernece z
2915       //
2916       "msig="<<meanS<<
2917       //
2918       "mphi="<<meanphi<<         //
2919       "rmsphi="<<rmsphi<<        //
2920       "gphi1="<<gphi1<<
2921       "gphi2="<<gphi2<<
2922       //
2923       "mphiP="<<meanphiP<<       //
2924       "rmsphiP="<<rmsphiP<<      //
2925       "gphiP1="<<gphiP1<<
2926       "gphiP2="<<gphiP2<<
2927       //
2928       "meanZ="<<meanZ<<
2929       "rmsZ="<<rmsZ<<
2930       "gz1="<<gz1<<
2931       "gz2="<<gz2<<
2932       //
2933       "gp31="<<gp31<<            //gaus mean - tgl
2934       "gp32="<<gp32<<            //gaus rms  -tgl
2935       "meanp3="<<meanp3<<            //mean - tgl
2936       "rmsp3="<<rmsp3<<            //rms  -tgl
2937       "gp41="<<gp41<<            //gaus mean - P4
2938       "gp42="<<gp42<<            //gaus rms  - P4
2939       "meanp4="<<meanp4<<            //mean - P4
2940       "rmsp4="<<rmsp4<<            //rms  - P4
2941       // Parameters from abs res analysis
2942       "SecIn="<<secInner<<              // inner sector
2943       "SecOut="<<secOuter<<             // outer sector
2944       "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
2945       "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
2946       "ibr.="<<&vecIBR<<   // radial filed integral
2947       "ibrphi.="<<&vecIBRPhi<<   // r=phifiled integral
2948       "ibr.="<<&vecIBR<<   // radial filed integral
2949       "ibz.="<<&vecIBZ<<   // z filed integral
2950       //
2951       "iblx.="<<&vecIBLX<<   // local bx  integral
2952       "ibly.="<<&vecIBLY<<   // local by integral
2953       "ibgx.="<<&vecIBGX<<   // global bx  integral
2954       "ibgy.="<<&vecIBGY<<   // global by integral
2955       //
2956       "X.="<<&vecX<<       // local x
2957       "Y.="<<&vecY<<       // local y
2958       "R.="<<&vecR<<       // radius
2959       "dY.="<<&vecDY<<     // abs y residuals
2960       "dZ.="<<&vecDZ<<     // abs z residuals
2961       "eY.="<<&vecEy<<     // error of y residuals
2962       "eZ.="<<&vecEz<<     // error of z residuals
2963       "kY.="<<&vecPhi<<    // local tangent y (fixed for sector)
2964       "kYR.="<<&vecPhiR<<  // tangent between laser and R vector (varies inside sector)
2965       "nCl.="<<&vecN<<     // number of clusters
2966       //
2967       "nClInY="<<nClInY<<               // Number of clusters for inner y
2968       "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
2969       "yAbsInSlope="<<yAbsInSlope <<  // fitted slope absres (inner y)
2970       "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
2971       "nClInZ="<<nClInZ<<               // Number of clusters for inner z
2972       "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
2973       "zAbsInSlope="<<zAbsInSlope <<  // fitted slope absres (inner z)
2974       "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
2975       //
2976       "nClOutY="<<nClOutY<<               // Number of clusters for outer y
2977       "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
2978       "yAbsOutSlope="<<yAbsOutSlope <<  // fitted slope absres (outer y)
2979       "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
2980       "nClOutZ="<<nClOutZ<<               // Number of clusters for outer z
2981       "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
2982       "zAbsOutSlope="<<zAbsOutSlope <<  // fitted slope absres (outer z)
2983       "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
2984       //
2985       "\n";
2986   }
2987   delete pcstream;
2988   /*
2989     Browse the content
2990     TFile fmean("laserMean.root");
2991
2992
2993    */
2994
2995
2996 }
2997
2998
2999
3000 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
3001   //
3002   //
3003   //
3004   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
3005   TFile * f = pcstream->GetFile();
3006   f->mkdir("dirphi");
3007   f->mkdir("dirphiP");
3008   f->mkdir("dirZ");
3009   TF1 fp("p1","pol1");
3010   //
3011   //
3012   char    cut[1000];
3013   char grname[1000];
3014   char grnamefull[1000];
3015
3016   Double_t mphi[100];
3017   Double_t mphiP[100];
3018   Double_t smphi[100];
3019   Double_t smphiP[100];
3020   Double_t mZ[100];
3021   Double_t smZ[100];
3022   Double_t bz[100];
3023   Double_t sbz[100];
3024   // fit parameters
3025   Double_t pphi[3];
3026   Double_t pphiP[3];
3027   Double_t pmZ[3];
3028
3029   //
3030   for (Int_t id=0; id<336; id++){
3031     // id =205;
3032     sprintf(cut,"fId==%d&&%s",id,cutUser);
3033     Int_t entries = chain->Draw("bz",cut,"goff");
3034     if (entries<3) continue;
3035     AliTPCLaserTrack *ltrp = 0;
3036     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
3037     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3038     Double_t lxyz[3];
3039     Double_t lpxyz[3];
3040     ltrp->GetXYZ(lxyz);
3041     ltrp->GetPxPyPz(lpxyz);
3042
3043     chain->Draw("bz",cut,"goff");
3044     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
3045     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
3046     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
3047     //
3048     chain->Draw("gphi1",cut,"goff");
3049     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
3050     chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
3051     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
3052     //
3053     chain->Draw("gphiP1",cut,"goff");
3054     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
3055     chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
3056     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
3057     //
3058     chain->Draw("gz1",cut,"goff");
3059     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
3060     chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
3061     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
3062     //
3063     //
3064     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
3065             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
3066     // store data
3067     // phi
3068     f->cd("dirphi");
3069     Float_t phiB0 =0;
3070     Float_t phiPB0=0;
3071     Float_t zB0=0;
3072     for (Int_t ientry=0; ientry<entries; ientry++){
3073       if (TMath::Abs(bz[ientry])<0.05){
3074         phiB0  = mphi[ientry];
3075         phiPB0 = mphiP[ientry];
3076         zB0    = mZ[ientry];
3077       }
3078     }
3079     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
3080     grphi->Draw("a*");
3081     grphi->Fit(&fp);
3082     pphi[0] = fp.GetParameter(0);                          // offset
3083     pphi[1] = fp.GetParameter(1);                          // slope
3084     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3085     sprintf(grname,"phi_id%d",id);
3086     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
3087     grphi->GetXaxis()->SetTitle("b_{z} (T)");
3088     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
3089     grphi->SetMaximum(1.2);
3090     grphi->SetMinimum(-1.2);
3091     grphi->Draw("a*");
3092
3093     grphi->Write();
3094     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
3095     // phiP
3096     f->cd("dirphiP)");
3097     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
3098     grphiP->Draw("a*");
3099     grphiP->Fit(&fp);
3100     pphiP[0] = fp.GetParameter(0);                          // offset
3101     pphiP[1] = fp.GetParameter(1);                          // slope
3102     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3103     sprintf(grname,"phiP_id%d",id);
3104     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
3105     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
3106     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
3107     grphiP->SetMaximum(pphiP[0]+0.005);
3108     grphiP->SetMinimum(pphiP[0]-0.005);
3109
3110     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
3111     grphiP->Write();
3112     //
3113     //Z
3114     f->cd("dirZ");
3115     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
3116     grmZ->Draw("a*");
3117     grmZ->Fit(&fp);
3118     pmZ[0] = fp.GetParameter(0);                          // offset
3119     pmZ[1] = fp.GetParameter(1);                          // slope
3120     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
3121     sprintf(grname,"mZ_id%d",id);
3122     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
3123     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
3124     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
3125
3126     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
3127     grmZ->Write();
3128     //
3129     // P4
3130     //
3131
3132     for (Int_t ientry=0; ientry<entries; ientry++){
3133       (*pcstream)<<"Mean"<<
3134         "id="<<id<<
3135         "LTr.="<<ltrp<<
3136         "entries="<<entries<<
3137         "bz="<<bz[ientry]<<
3138         "lx0="<<lxyz[0]<<          // reference x
3139         "lx1="<<lxyz[1]<<          // reference y
3140         "lx2="<<lxyz[2]<<          // refernece z
3141         "lpx0="<<lpxyz[0]<<          // reference x
3142         "lpx1="<<lpxyz[1]<<          // reference y
3143         "lpx2="<<lpxyz[2]<<          // refernece z
3144         //values
3145         "phiB0="<<phiB0<<          // position shift at 0 field
3146         "phiPB0="<<phiPB0<<        // angular  shift at 0 field
3147         "zB0="<<zB0<<              // z shift for 0 field
3148         //
3149         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
3150         "pphi0="<<pphi[0]<<   // offset
3151         "pphi1="<<pphi[1]<<   // slope
3152         "pphi2="<<pphi[2]<<   // norm chi2
3153         //
3154         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
3155         "pphiP0="<<pphiP[0]<< // offset
3156         "pphiP1="<<pphiP[1]<< // slope
3157         "pphiP2="<<pphiP[2]<< // norm chi2
3158         //
3159         "gz1="<<mZ[ientry]<<
3160         "pmZ0="<<pmZ[0]<<     // offset
3161         "pmZ1="<<pmZ[1]<<     // slope
3162         "pmZ2="<<pmZ[2]<<     // norm chi2
3163         "\n";
3164     }
3165   }
3166
3167   delete pcstream;
3168
3169 }
3170
3171
3172 void AliTPCcalibLaser::Analyze(){
3173   //
3174   //
3175   //
3176 }
3177
3178
3179 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
3180
3181   TIterator* iter = li->MakeIterator();
3182   AliTPCcalibLaser* cal = 0;
3183   static Int_t counter0=0;
3184   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
3185     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
3186       //      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
3187       return -1;
3188     }
3189     printf("Marging number %d\n", counter0);
3190     counter0++;
3191     //
3192     MergeFitHistos(cal);
3193     TH1F *h=0x0;
3194     TH1F *hm=0x0;
3195     TH2F *h2=0x0;
3196     TH2F *h2m=0x0;
3197     //    TProfile *hp=0x0;
3198     //TProfile *hpm=0x0;
3199
3200     for (Int_t id=0; id<336; id++){
3201       // merge fDeltaZ histograms
3202       hm = (TH1F*)cal->fDeltaZ.At(id);
3203       h  = (TH1F*)fDeltaZ.At(id);
3204       if (!h) {
3205         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3206         h->SetDirectory(0);
3207         fDeltaZ.AddAt(h,id);
3208       }
3209       if (hm) h->Add(hm);
3210       // merge fP3 histograms
3211       hm = (TH1F*)cal->fDeltaP3.At(id);
3212       h  = (TH1F*)fDeltaP3.At(id);
3213       if (!h) {
3214         h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3215         h->SetDirectory(0);
3216         fDeltaP3.AddAt(h,id);
3217       }
3218       if (hm) h->Add(hm);
3219       // merge fP4 histograms
3220       hm = (TH1F*)cal->fDeltaP4.At(id);
3221       h  = (TH1F*)fDeltaP4.At(id);
3222       if (!h) {
3223         h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3224         h->SetDirectory(0);
3225         fDeltaP4.AddAt(h,id);
3226       }
3227       if (hm) h->Add(hm);
3228
3229       //
3230       // merge fDeltaPhi histograms
3231       hm = (TH1F*)cal->fDeltaPhi.At(id);
3232       h  = (TH1F*)fDeltaPhi.At(id);
3233       if (!h) {
3234         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3235         h->SetDirectory(0);
3236         fDeltaPhi.AddAt(h,id);
3237       }
3238       if (hm) h->Add(hm);
3239       // merge fDeltaPhiP histograms
3240       hm = (TH1F*)cal->fDeltaPhiP.At(id);
3241       h  = (TH1F*)fDeltaPhiP.At(id);
3242       if (!h) {
3243         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3244         h->SetDirectory(0);
3245         fDeltaPhiP.AddAt(h,id);
3246       }
3247       if (hm) h->Add(hm);
3248       // merge fSignals histograms
3249       hm = (TH1F*)cal->fSignals.At(id);
3250       h  = (TH1F*)fSignals.At(id);
3251       if (!h) {
3252         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3253         h->SetDirectory(0);
3254         fSignals.AddAt(h,id);
3255       }
3256       if (hm) h->Add(hm);
3257       //
3258       //
3259       // merge ProfileY histograms -0
3260       h2m = (TH2F*)cal->fDeltaYres.At(id);
3261       h2  = (TH2F*)fDeltaYres.At(id);
3262       if (h2m&&h2) h2->Add(h2m);
3263       //
3264       h2m = (TH2F*)cal->fDeltaZres.At(id);
3265       h2  = (TH2F*)fDeltaZres.At(id);
3266       if (h2m&&h2) h2->Add(h2m);
3267       // merge ProfileY histograms - 2
3268       h2m = (TH2F*)cal->fDeltaYres2.At(id);
3269       h2  = (TH2F*)fDeltaYres2.At(id);
3270       if (h2m&&h2) h2->Add(h2m);
3271       //
3272       h2m = (TH2F*)cal->fDeltaZres2.At(id);
3273       h2  = (TH2F*)fDeltaZres2.At(id);
3274       if (h2m&&h2) h2->Add(h2m);
3275
3276       // merge ProfileY histograms - abs
3277       h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
3278       h2  = (TH2F*)fDeltaYresAbs.At(id);
3279       if (h2m&&h2) h2->Add(h2m);
3280
3281       h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
3282       h2  = (TH2F*)fDeltaZresAbs.At(id);
3283       if (h2m&&h2) h2->Add(h2m);
3284       // merge ProfileY histograms - 3
3285       //h2m = (TH2F*)cal->fDeltaYres3.At(id);
3286       //h2  = (TH2F*)fDeltaYres3.At(id);
3287       //if (h2m) h2->Add(h2m);
3288       //
3289       //h2m = (TH2F*)cal->fDeltaZres3.At(id);
3290       //h2  = (TH2F*)fDeltaZres3.At(id);
3291       //if (h2m) h->Add(h2m);
3292       //
3293       //
3294     }
3295   }
3296   return 0;
3297 }
3298
3299 void   AliTPCcalibLaser::MakeFitHistos(){
3300   //
3301   // Make a fit histograms
3302   //
3303   // Number of clusters
3304   //
3305   //TH2F           *fHisNclIn;      //->Number of clusters inner
3306   //TH2F           *fHisNclOut;     //->Number of clusters outer
3307   //TH2F           *fHisNclIO;      //->Number of cluster inner outer
3308   // TH2F           *fHisdEdx;       //->dEdx histo
3309   fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
3310   fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
3311   fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
3312   //
3313   fHisLclIn  = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
3314   fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
3315   fHisLclIO  = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
3316   //
3317   fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
3318   fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
3319
3320   //
3321   //  Chi2
3322   //
3323   //   TH2F           *fHisChi2YIn1;      //->chi2 y inner - line
3324   //   TH2F           *fHisChi2YOut1;     //->chi2 y inner - line
3325   //   TH2F           *fHisChi2YIn2;      //->chi2 y inner - parabola
3326   //   TH2F           *fHisChi2YOut2;     //->chi2 y inner - parabola
3327   //   TH2F           *fHisChi2YIO1;      //->chi2 y IO    - common
3328   fHisChi2YIn1   = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
3329   fHisChi2YOut1  = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
3330   fHisChi2YIn2   = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
3331   fHisChi2YOut2  = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
3332   fHisChi2YIO1   = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
3333   //   TH2F           *fHisChi2ZIn1;      //->chi2 z inner - line
3334   //   TH2F           *fHisChi2ZOut1;     //->chi2 z inner - line
3335   //   TH2F           *fHisChi2ZIn2;      //->chi2 z inner - parabola
3336   //   TH2F           *fHisChi2ZOut2;     //->chi2 z inner - parabola
3337   //   TH2F           *fHisChi2ZIO1;      //->chi2 z IO    - common
3338   fHisChi2ZIn1   = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
3339   fHisChi2ZOut1  = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
3340   fHisChi2ZIn2   = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
3341   fHisChi2ZOut2  = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
3342   fHisChi2ZIO1   = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
3343   //
3344   // Fit
3345   //
3346   //
3347   //   TH2F           *fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
3348   //   TH2F           *fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
3349   //   TH2F           *fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
3350   //   TH2F           *fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
3351   //   TH2F           *fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
3352   //   TH2F           *fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
3353   //   TH2F           *fHisPy2vP2In;   //-> Curv  P2inner - parabola
3354   //   TH2F           *fHisPy2vP2Out;  //-> Curv  P2outer - parabola
3355   //   TH2F           *fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
3356   fHisPy1vP0   = new TH2F("HisPy1vP0",   "HisPy1vP0",336,0,336,500,-0.3,0.3);
3357   fHisPy2vP0   = new TH2F("HisPy2vP0",   "HisPy2vP0",336,0,336,500,-0.3,0.3);
3358   fHisPy3vP0   = new TH2F("HisPy3vP0",   "HisPy3vP0",336,0,336,500,-0.3,0.3);
3359   fHisPy1vP1   = new TH2F("HisPy1vP1",   "HisPy1vP1",336,0,336,500,-0.01,0.01);
3360   fHisPy2vP1   = new TH2F("HisPy2vP1",   "HisPy2vP1",336,0,336,500,-0.01,0.01);
3361   fHisPy3vP1   = new TH2F("HisPy3vP1",   "HisPy3vP1",336,0,336,500,-0.01,0.01);
3362   fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
3363   fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
3364   fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
3365   //
3366   //
3367   //   TH2F           *fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
3368   //   TH2F           *fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
3369   //   TH2F           *fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
3370   //   TH2F           *fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
3371   //   TH2F           *fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
3372   //   TH2F           *fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
3373   //   TH2F           *fHisPz2vP2In;   //-> Curv  P2inner - parabola
3374   //   TH2F           *fHisPz2vP2Out;  //-> Curv  P2outer - parabola
3375   //   TH2F           *fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
3376   fHisPz1vP0   = new TH2F("HisPz1vP0",   "HisPz1vP0",336,0,336,500,-0.3,0.3);
3377   fHisPz2vP0   = new TH2F("HisPz2vP0",   "HisPz2vP0",336,0,336,500,-0.3,0.3);
3378   fHisPz3vP0   = new TH2F("HisPz3vP0",   "HisPz3vP0",336,0,336,500,-0.3,0.3);
3379   fHisPz1vP1   = new TH2F("HisPz1vP1",   "HisPz1vP1",336,0,336,500,-0.01,0.01);
3380   fHisPz2vP1   = new TH2F("HisPz2vP1",   "HisPz2vP1",336,0,336,500,-0.01,0.01);
3381   fHisPz3vP1   = new TH2F("HisPz3vP1",   "HisPz3vP1",336,0,336,500,-0.01,0.01);
3382   fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
3383   fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
3384   fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
3385
3386   fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
3387   fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
3388
3389   fHisNclIn->SetDirectory(0);      //->Number of clusters inner
3390   fHisNclOut->SetDirectory(0);     //->Number of clusters outer
3391   fHisNclIO->SetDirectory(0);      //->Number of cluster inner outer
3392   fHisLclIn->SetDirectory(0);      //->Level arm inner
3393   fHisLclOut->SetDirectory(0);     //->Level arm outer
3394   fHisLclIO->SetDirectory(0);      //->Number of cluster inner outer
3395   fHisdEdx->SetDirectory(0);      //->Number of cluster inner outer
3396   fHisdZfit->SetDirectory(0);      //->Number of cluster inner outer
3397   //
3398   //
3399   fHisChi2YIn1->SetDirectory(0);      //->chi2 y inner - line
3400   fHisChi2YOut1->SetDirectory(0);     //->chi2 y inner - line
3401   fHisChi2YIn2->SetDirectory(0);      //->chi2 y inner - parabola
3402   fHisChi2YOut2->SetDirectory(0);     //->chi2 y inner - parabola
3403   fHisChi2YIO1->SetDirectory(0);      //->chi2 y IO    - common
3404   fHisChi2ZIn1->SetDirectory(0);      //->chi2 z inner - line
3405   fHisChi2ZOut1->SetDirectory(0);     //->chi2 z inner - line
3406   fHisChi2ZIn2->SetDirectory(0);      //->chi2 z inner - parabola
3407   fHisChi2ZOut2->SetDirectory(0);     //->chi2 z inner - parabola
3408   fHisChi2ZIO1->SetDirectory(0);      //->chi2 z IO    - common
3409   //
3410   //
3411   fHisPy1vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - line
3412   fHisPy2vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - parabola
3413   fHisPy3vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - common parabola
3414   fHisPy1vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - line
3415   fHisPy2vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - parabola
3416   fHisPy3vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - common parabola
3417   fHisPy2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
3418   fHisPy2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
3419   fHisPy3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
3420   //
3421   //
3422   fHisPz1vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - line
3423   fHisPz2vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - parabola
3424   fHisPz3vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - common parabola
3425   fHisPz1vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - line
3426   fHisPz2vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - parabola
3427   fHisPz3vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - common parabola
3428   fHisPz2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
3429   fHisPz2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
3430   fHisPz3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
3431
3432   fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
3433   fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
3434
3435
3436
3437   //
3438   //
3439   //
3440   for (Int_t id=0; id<336;id++){
3441     TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
3442     TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
3443     TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
3444     TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
3445     TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3446     TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
3447     //    TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
3448     //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
3449     if (!profy){
3450       profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
3451       profy->SetDirectory(0);
3452       fDeltaYres.AddAt(profy,id);
3453       profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
3454       profy2->SetDirectory(0);
3455       fDeltaYres2.AddAt(profy2,id);
3456       if(!fUseFixedDriftV)
3457         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
3458       else
3459         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
3460       profyabs->SetDirectory(0);
3461       fDeltaYresAbs.AddAt(profyabs,id);
3462       //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3463       //profy3->SetDirectory(0);
3464       //fDeltaYres3.AddAt(profy3,id);
3465     }
3466     if (!profz){
3467       profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
3468       profz->SetDirectory(0);
3469       fDeltaZres.AddAt(profz,id);
3470       profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
3471       profz2->SetDirectory(0);
3472       fDeltaZres2.AddAt(profz2,id);
3473       if(!fUseFixedDriftV)
3474         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
3475       else
3476         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
3477       profzabs->SetDirectory(0);
3478       fDeltaZresAbs.AddAt(profzabs,id);
3479       //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
3480       //profz3->SetDirectory(0);
3481       //fDeltaZres3.AddAt(profz3,id);
3482     }
3483   }
3484   //
3485   //
3486   for (Int_t id=0; id<336;id++){
3487     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
3488     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
3489     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
3490
3491     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
3492     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
3493     TH1F * hisSignal = (TH1F*)fSignals.At(id);
3494
3495     if (!hisdz){
3496       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
3497       hisdz->SetDirectory(0);
3498       fDeltaZ.AddAt(hisdz,id);
3499
3500       hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
3501       hisP3->SetDirectory(0);
3502       fDeltaP3.AddAt(hisP3,id);
3503       //
3504       hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
3505       hisP4->SetDirectory(0);
3506       fDeltaP4.AddAt(hisP4,id);
3507
3508       //
3509       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
3510       hisdphi->SetDirectory(0);
3511       fDeltaPhi.AddAt(hisdphi,id);
3512       //
3513       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
3514       hisdphiP->SetDirectory(0);
3515       fDeltaPhiP.AddAt(hisdphiP,id);
3516       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
3517       hisSignal->SetDirectory(0);
3518       fSignals.AddAt(hisSignal,id);
3519     }
3520   }
3521
3522   SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
3523   SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
3524   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
3525   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
3526   //
3527   // Make THnSparse
3528   //
3529   //                   id   side   rod  bundle  beam  dP0  dP1  dP2  dP3  dP4 ncl dEdx
3530   Int_t binsLaser[12]= {336,  //id
3531                         2,    //side
3532                         6,    //rod
3533                         4,    //bundle
3534                         7,    //beam
3535                         300,  //dP0
3536                         300,  //dP1
3537                         300,  //dP2
3538                         300,  //dP3
3539                         300,  //dP4
3540                         80,   //ncl
3541                         50};   //dEdx
3542   Double_t xminLaser[12]= {0,  //id
3543                         0,     //side
3544                         0,     //rod
3545                         0,     //bundle
3546                         0,     //beam
3547                         -1,    //dP0
3548                         -1,    //dP1
3549                         -0.01,   //dP2
3550                         -0.01,  //dP3
3551                         -0.1,  //dP4
3552                         0,   //ncl
3553                         0};   //sqrt dEdx
3554   Double_t xmaxLaser[12]= {336,  //id
3555                         2,    //side
3556                         6,    //rod
3557                         4,    //bundle
3558                         7,    //beam
3559                         1,  //dP0
3560                         1,  //dP1
3561                         0.01,  //dP2
3562                         0.01,  //dP3
3563                         0.1,  //dP4
3564                         160,   //ncl
3565                         40};   //sqrt dEdx
3566
3567   TString nameLaser[12]= {"id",
3568                           "side",
3569                           "rod",
3570                           "bundle",
3571                           "beam",
3572                           "dP0",
3573                           "dP1",
3574                           "dP2",
3575                           "dP3",
3576                           "dP4",
3577                           "ncl",
3578                           "sqrt dEdx"};
3579   TString titleLaser[12]= {"id",
3580                           "side",
3581                           "rod",
3582                           "bundle",
3583                           "beam",
3584                           "#Delta_{P0}",
3585                           "#Delta_{P1}",
3586                           "#Delta_{P2}",
3587                           "#Delta_{P3}",
3588                           "#Delta_{P4}",
3589                           "N_{cl}",
3590                           "#sqrt{dEdx}"};
3591   fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
3592   for (Int_t iaxis=1; iaxis<12; iaxis++){
3593     fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
3594     fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
3595   }
3596 }
3597
3598 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
3599   //
3600   // Merge content of histograms
3601   //
3602   // Only first histogram is checked - all other should be the same
3603   if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
3604
3605   if (!laser->fHisNclIn) return;  // empty histograms
3606   if (!fHisNclIn) MakeFitHistos();
3607   //
3608   fHisNclIn->Add(laser->fHisNclIn  );      //->Number of clusters inner
3609   fHisNclOut->Add(laser->fHisNclOut  );     //->Number of clusters outer
3610   fHisNclIO->Add(laser->fHisNclIO  );      //->Number of cluster inner outer
3611   fHisLclIn->Add(laser->fHisLclIn  );      //->Level arm inner
3612   fHisLclOut->Add(laser->fHisLclOut  );     //->Level arm outer
3613   fHisLclIO->Add(laser->fHisLclIO  );      //->Number of cluster inner outer
3614   fHisdEdx->Add(laser->fHisdEdx  );      //->dedx
3615   fHisdZfit->Add(laser->fHisdZfit  );      //->dedx
3616   //
3617   //
3618   fHisChi2YIn1->Add(laser->fHisChi2YIn1  );      //->chi2 y inner - line
3619   fHisChi2YOut1->Add(laser->fHisChi2YOut1  );     //->chi2 y inner - line
3620   fHisChi2YIn2->Add(laser->fHisChi2YIn2  );      //->chi2 y inner - parabola
3621   fHisChi2YOut2->Add(laser->fHisChi2YOut2  );     //->chi2 y inner - parabola
3622   fHisChi2YIO1->Add(laser->fHisChi2YIO1  );      //->chi2 y IO    - common
3623   fHisChi2ZIn1->Add(laser->fHisChi2ZIn1  );      //->chi2 z inner - line
3624   fHisChi2ZOut1->Add(laser->fHisChi2ZOut1  );     //->chi2 z inner - line
3625   fHisChi2ZIn2->Add(laser->fHisChi2ZIn2  );      //->chi2 z inner - parabola
3626   fHisChi2ZOut2->Add(laser->fHisChi2ZOut2  );     //->chi2 z inner - parabola
3627   fHisChi2ZIO1->Add(laser->fHisChi2ZIO1  );      //->chi2 z IO    - common
3628   //
3629   //
3630   fHisPy1vP0->Add(laser->fHisPy1vP0  );     //-> delta y   P0outer-P0inner - line
3631   fHisPy2vP0->Add(laser->fHisPy2vP0  );     //-> delta y   P0outer-P0inner - parabola
3632   fHisPy3vP0->Add(laser->fHisPy3vP0  );     //-> delta y   P0outer-P0inner - common parabola
3633   fHisPy1vP1->Add(laser->fHisPy1vP1  );     //-> delta ky  P1outer-P1inner - line
3634   fHisPy2vP1->Add(laser->fHisPy2vP1  );     //-> delta ky  P1outer-P1inner - parabola
3635   fHisPy3vP1->Add(laser->fHisPy3vP1  );     //-> delta ky  P1outer-P1inner - common parabola
3636   fHisPy2vP2In->Add(laser-> fHisPy2vP2In  );   //-> Curv  P2inner - parabola
3637   fHisPy2vP2Out->Add(laser->fHisPy2vP2Out  );  //-> Curv  P2outer - parabola
3638   fHisPy3vP2IO->Add(laser->fHisPy3vP2IO  );   //-> Curv  P2outerinner - common parabola
3639   //
3640   //
3641   fHisPz1vP0->Add(laser->fHisPz1vP0  );     //-> delta z   P0outer-P0inner - line
3642   fHisPz2vP0->Add(laser->fHisPz2vP0  );     //-> delta z   P0outer-P0inner - parabola
3643   fHisPz3vP0->Add(laser->fHisPz3vP0  );     //-> delta z   P0outer-P0inner - common parabola
3644   fHisPz1vP1->Add(laser->fHisPz1vP1  );     //-> delta kz  P1outer-P1inner - line
3645   fHisPz2vP1->Add(laser->fHisPz2vP1  );     //-> delta kz  P1outer-P1inner - parabola
3646   fHisPz3vP1->Add(laser->fHisPz3vP1  );     //-> delta kz  P1outer-P1inner - common parabola
3647   fHisPz2vP2In->Add(laser->fHisPz2vP2In  );   //-> Curv  P2inner - parabola
3648   fHisPz2vP2Out->Add(laser->fHisPz2vP2Out  );  //-> Curv  P2outer - parabola
3649   fHisPz3vP2IO->Add(laser->fHisPz3vP2IO  );   //-> Curv  P2outerinner - common parabola
3650   fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
3651   fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
3652   //
3653   //
3654   //
3655
3656
3657
3658
3659 }
3660
3661
3662
3663
3664 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
3665   //
3666   // Dump fit information - collect information from the streamers
3667   //
3668   /*
3669     TChain * chainFit=0;
3670     TChain * chainTrack=0;
3671     TChain * chain=0;
3672     //
3673     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3674     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
3675     AliXRDPROOFtoolkit tool;
3676     chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3677     chainTrack->Lookup();
3678     chainTrack->SetProof(kTRUE);
3679     chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
3680     chainDrift->Lookup();
3681     chainDrift->SetProof(kTRUE);
3682
3683     chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3684     chain->Lookup();
3685     chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3686     chainFit->Lookup();
3687     chainFit->SetProof(kTRUE);
3688     chain->SetProof(kTRUE);
3689     AliTPCLaserTrack::LoadTracks();
3690     //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
3691
3692   */
3693   //
3694   // Fit cuts
3695   //
3696   TCut cutP3("abs(Tr.fP[3])<0.1");
3697   TCut cutP4("abs(Tr.fP[4])<0.5");
3698   TCut cutPx = cutP3+cutP4;
3699   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
3700   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
3701   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
3702   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
3703   //
3704   TCut cutdEdx("sqrt(dEdx)>3");
3705   TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
3706   TCut cutN("nclO>20&&nclI>20");
3707   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
3708   //
3709   // Cluster cuts
3710   //
3711   TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
3712   TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
3713   TCut cutClX("abs(Cl[].fX)>10");
3714   TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
3715   TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
3716   TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
3717   TCut cutQ("sqrt(Cl[].fMax)>4");
3718   TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
3719
3720
3721   TH1F * phisAl     = 0;
3722   TH1F * phisAccept = 0;
3723   TH1F * phisOut    = 0;
3724   TProfile * pdEdx  = 0;
3725
3726   TProfile * pP0    = 0;
3727   TProfile * pP1    = 0;
3728   TProfile * pP2    = 0;
3729   TProfile * pP3    = 0;
3730   TProfile * pP4    = 0;
3731   //
3732   TProfile * pNclI  = 0;
3733   TProfile * pNclO  = 0;
3734   //
3735   TProfile * pchi2YIn   =0;
3736   TProfile * pchi2ZIn   =0;
3737   TProfile * pchi2YOut  =0;
3738   TProfile * pchi2ZOut  =0;
3739   TProfile * pchi2YInOut =0;
3740   TProfile * pchi2ZInOut =0;;
3741   // laser counters
3742   chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
3743   phisAl = (TH1F*)gROOT->FindObject("hisAl");
3744   chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
3745   phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
3746   chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
3747   phisOut = (TH1F*)gROOT->FindObject("hisOut");
3748   //
3749   chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
3750   pdEdx   = (TProfile*)gROOT->FindObject("hdEdx");
3751   // track param
3752   //
3753   chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3754   pP0   = (TProfile*)gROOT->FindObject("hP0");
3755   chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3756   pP1   = (TProfile*)gROOT->FindObject("hP1");
3757   chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3758   pP2   = (TProfile*)gROOT->FindObject("hP2");
3759   chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3760   pP3   = (TProfile*)gROOT->FindObject("hP3");
3761   chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
3762   pP4   = (TProfile*)gROOT->FindObject("hP4");
3763   //
3764   chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
3765   pNclI   = (TProfile*)gROOT->FindObject("hNclI");
3766   chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
3767   pNclO   = (TProfile*)gROOT->FindObject("hNclO");
3768   //
3769   //
3770   chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
3771   pchi2YIn   = (TProfile*)gROOT->FindObject("hChi2YIn");
3772   chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
3773   pchi2YOut   = (TProfile*)gROOT->FindObject("hChi2YOut");
3774   chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
3775   pchi2YInOut   = (TProfile*)gROOT->FindObject("hChi2YInOut");
3776   chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
3777   pchi2ZIn   = (TProfile*)gROOT->FindObject("hChi2ZIn");
3778   chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
3779   pchi2ZOut   = (TProfile*)gROOT->FindObject("hChi2ZOut");
3780   chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
3781   pchi2ZInOut   = (TProfile*)gROOT->FindObject("hChi2ZInOut");
3782   //
3783   // second derivatives
3784   //
3785   TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3786   chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3787   TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3788   chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3789   TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3790   chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3791   //
3792   phisPy2In->FitSlicesY();
3793   TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3794   TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3795   TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3796   //
3797   phisPy2Out->FitSlicesY();
3798   TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3799   TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3800   TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3801   //
3802   phisPy2InOut->FitSlicesY();
3803   TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3804   TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3805   TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3806   //
3807   TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3808   chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3809   TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3810   chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3811   TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3812   chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3813   //
3814   phisPz2In->FitSlicesY();
3815   TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3816   TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3817   TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3818   //
3819   phisPz2Out->FitSlicesY();
3820   TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3821   TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3822   TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3823   //
3824   phisPz2InOut->FitSlicesY();
3825   TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3826   TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3827   TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3828   //
3829   //
3830   //
3831
3832
3833   {
3834     TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3835     for (Int_t ilaser=0; ilaser<336; ilaser++){
3836       Float_t all=phisAl->GetBinContent(ilaser+1);
3837       Float_t accept=phisAccept->GetBinContent(ilaser+1);
3838       Float_t out=phisOut->GetBinContent(ilaser+1);
3839       Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3840       Float_t mP0 = pP0->GetBinContent(ilaser+1);
3841       Float_t mP1 = pP1->GetBinContent(ilaser+1);
3842       Float_t mP2 = pP2->GetBinContent(ilaser+1);
3843       Float_t mP3 = pP3->GetBinContent(ilaser+1);
3844       Float_t mP4 = pP4->GetBinContent(ilaser+1);
3845
3846
3847       Float_t nclI  = pNclI->GetBinContent(ilaser+1);
3848       Float_t nclO  = pNclO->GetBinContent(ilaser+1);
3849       //
3850       Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3851       Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3852       Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3853       Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3854       Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3855       Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3856       //
3857       Float_t entriesPy2In  = phisPy2InEntries->GetBinContent(ilaser+1);
3858       Float_t meanPy2In     = phisPy2InMean->GetBinContent(ilaser+1);
3859       Float_t sigmaPy2In    = phisPy2InSigma->GetBinContent(ilaser+1);
3860       //
3861       Float_t entriesPy2Out  = phisPy2OutEntries->GetBinContent(ilaser+1);
3862       Float_t meanPy2Out     = phisPy2OutMean->GetBinContent(ilaser+1);
3863       Float_t sigmaPy2Out    = phisPy2OutSigma->GetBinContent(ilaser+1);
3864       //
3865       Float_t entriesPy2InOut  = phisPy2InOutEntries->GetBinContent(ilaser+1);
3866       Float_t meanPy2InOut     = phisPy2InOutMean->GetBinContent(ilaser+1);
3867       Float_t sigmaPy2InOut    = phisPy2InOutSigma->GetBinContent(ilaser+1);
3868       //
3869       Float_t entriesPz2In  = phisPz2InEntries->GetBinContent(ilaser+1);
3870       Float_t meanPz2In     = phisPz2InMean->GetBinContent(ilaser+1);
3871       Float_t sigmaPz2In    = phisPz2InSigma->GetBinContent(ilaser+1);
3872       //
3873       Float_t entriesPz2Out  = phisPz2OutEntries->GetBinContent(ilaser+1);
3874       Float_t meanPz2Out     = phisPz2OutMean->GetBinContent(ilaser+1);
3875       Float_t sigmaPz2Out    = phisPz2OutSigma->GetBinContent(ilaser+1);
3876       //
3877       Float_t entriesPz2InOut  = phisPz2InOutEntries->GetBinContent(ilaser+1);
3878       Float_t meanPz2InOut     = phisPz2InOutMean->GetBinContent(ilaser+1);
3879       Float_t sigmaPz2InOut    = phisPz2InOutSigma->GetBinContent(ilaser+1);
3880
3881       AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3882       (*pcstream)<<"Scan"<<
3883         "Run="<<id<<
3884         "LTr.="<<ltrp<<
3885         "all="<<all<<
3886         "accept="<<accept<<
3887         "out="<<out<<
3888         "sdedx="<<sdedx<<
3889         "mP0="<<mP0<<
3890         "mP1="<<mP1<<
3891         "mP2="<<mP2<<
3892         "mP3="<<mP3<<
3893         "mP4="<<mP4<<
3894         "nclI="<<nclI<<
3895         "nclO="<<nclO<<
3896         "chi2YIn="<<chi2YIn<<
3897         "chi2YOut="<<chi2YOut<<
3898         "chi2YInOut="<<chi2YInOut<<
3899         "chi2ZIn="<<chi2ZIn<<
3900         "chi2ZOut="<<chi2ZOut<<
3901         "chi2ZInOut="<<chi2ZInOut<<
3902         //
3903         "nPy2In="<<entriesPy2In<<
3904         "mPy2In="<<meanPy2In<<
3905         "sPy2In="<<sigmaPy2In<<
3906         //
3907         "nPy2Out="<<entriesPy2Out<<
3908         "mPy2Out="<<meanPy2Out<<
3909         "sPy2Out="<<sigmaPy2Out<<
3910         //
3911         "nPy2InOut="<<entriesPy2InOut<<
3912         "mPy2InOut="<<meanPy2InOut<<
3913         "sPy2InOut="<<sigmaPy2InOut<<
3914         //
3915         "nPz2In="<<entriesPz2In<<
3916         "mPz2In="<<meanPz2In<<
3917         "sPz2In="<<sigmaPz2In<<
3918         //
3919         "nPz2Out="<<entriesPz2Out<<
3920         "mPz2Out="<<meanPz2Out<<
3921         "sPz2Out="<<sigmaPz2Out<<
3922         //
3923         "nPz2InOut="<<entriesPz2InOut<<
3924         "mPz2InOut="<<meanPz2InOut<<
3925         "sPz2InOut="<<sigmaPz2InOut<<
3926         "\n";
3927     }
3928
3929     delete pcstream;
3930   }
3931 }
3932
3933 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3934                                          TVectorD& meanSlope,
3935                                          TVectorD& sectorArray,
3936                                          Int_t option)
3937 {
3938   // This method should ideally go in AliTPCLaser
3939   // option == 0 (pads outer - closest to beam)
3940   // option == 1 (pads inner)
3941   // option == 2 (time outer)
3942   // option == 3 (time inner)
3943   Int_t nFailures = 0;
3944
3945   for(Int_t id = 0; id < 336; id++) {
3946
3947     if (!AliTPCLaserTrack::GetTracks())
3948       AliTPCLaserTrack::LoadTracks();
3949     AliTPCLaserTrack *ltrp =
3950       (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3951
3952     AliExternalTrackParam trackParam(*ltrp);
3953
3954     Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3955     if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3956       deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3957
3958     Double_t angle = trackParam.GetAlpha();
3959     if(angle<0)
3960       angle += 2*TMath::Pi();
3961     if(trackParam.GetSnp()>0) // track points to sector "before"
3962       angle -= deltaangle*TMath::DegToRad();
3963     else // track points to sector "after"
3964       angle += deltaangle*TMath::DegToRad();
3965
3966     Bool_t success = trackParam.Rotate(angle);
3967
3968     if(!success) {
3969       //      cout << "WARNING: Rotate failed for ID: " << id << endl;
3970       nFailures++;
3971     }
3972
3973     angle *= TMath::RadToDeg();
3974
3975     Int_t sector = TMath::Nint((angle-10.0)/20.0);
3976     if(sector<0)
3977       sector += 18;
3978     else if(sector>=18)
3979       sector -= 18;
3980     if(ltrp->GetSide()==1) // C side
3981       sector += 18;
3982     if(option==0 || option==2)
3983       sector += 36;
3984     if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
3985       sector += 36;
3986
3987     sectorArray[id] = sector;
3988
3989     const Double_t x0 = 0;
3990
3991     Double_t slopey  = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3992     Double_t slopez  = trackParam.GetTgl();
3993     //  One needs a factor sqrt(1+slopey**2) to take into account the
3994     //  longer path length
3995     slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3996     if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3997       slopez *= -1;
3998     //    Double_t offsetz = trackParam.GetZ();
3999     Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
4000     Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
4001     if(option==2 || option==3) {
4002       meanOffset[id] = offsetz; meanSlope[id] = slopez;
4003     } else {
4004       meanOffset[id] = offsety; meanSlope[id] = slopey;
4005     }
4006   }
4007
4008   if(nFailures>0)
4009     AliWarning(Form("Rotate method failed %d times", nFailures));
4010 }
4011
4012
4013   /*
4014     TFile f("vscan.root");
4015    */
4016
4017   /*
4018     pad binning effect
4019    chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
4020    //
4021    chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
4022    //
4023
4024 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);
4025
4026
4027 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)
4028
4029   */
4030
4031
4032
4033
4034
4035   /*
4036   // check edge effects
4037   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4038   //
4039   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4040
4041   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)
4042
4043
4044
4045   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4046   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
4047
4048 */
4049
4050
4051
4052
4053
4054   /*
4055     Edge y effect
4056
4057     dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
4058
4059
4060     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)
4061
4062     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)
4063
4064
4065
4066
4067
4068     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)
4069
4070     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)
4071
4072
4073
4074     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)
4075
4076     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)
4077
4078
4079
4080 */
4081
4082
4083 /*
4084
4085 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4086
4087 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
4088
4089
4090
4091 chainFit->Draw("LTr.fId","nclI>10",100000)
4092
4093 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
4094
4095 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
4096
4097 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
4098
4099 */
4100
4101
4102
4103
4104
4105
4106 /*
4107  gSystem->Load("libSTAT.so")
4108  TStatToolkit toolkit;
4109  Double_t chi2;
4110  TVectorD fitParam;
4111  TMatrixD covMatrix;
4112  Int_t npoints;
4113
4114  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
4115
4116
4117 TString fstring="";
4118 //
4119 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
4120 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
4121 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
4122 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4
4123 //
4124 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
4125 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
4126 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
4127 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
4128 //
4129 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
4130 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
4131 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
4132 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
4133
4134
4135
4136
4137  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4138
4139  treeT->SetAlias("fit",strq0->Data());
4140
4141
4142  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
4143
4144  treeT->SetAlias("fitP",strqP->Data());
4145
4146
4147  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
4148  treeT->SetAlias("fitD",strqDrift->Data());
4149
4150
4151 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
4152 {
4153 for (Int_t i=0; i<6;i++){
4154 treeT->SetLineColor(i+2);
4155 treeT->SetMarkerSize(1);
4156 treeT->SetMarkerStyle(22+i);
4157 treeT->SetMarkerColor(i+2);
4158
4159 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
4160 }
4161 }
4162  */
4163
4164
4165
4166 /*
4167   TTree * tree = (TTree*)f.Get("FitModels");
4168
4169   TEventList listLFit0("listLFit0","listLFit0");
4170   TEventList listLFit1("listLFit1","listLFit1");
4171   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
4172   tree->SetEventList(&listLFit0);
4173
4174
4175
4176
4177   gSystem->Load("libSTAT.so")
4178   TStatToolkit toolkit;
4179   Double_t chi2;
4180   TVectorD fitParam;
4181   TMatrixD covMatrix;
4182   Int_t npoints;
4183
4184   chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
4185   chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
4186
4187
4188   TString fstring="";
4189   fstring+="cos(dp)++";
4190   fstring+="sin(dp)++";
4191   fstring+="cos(dt)++";
4192   fstring+="sin(dt)++";
4193
4194   TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
4195
4196
4197
4198 */
4199
4200
4201
4202 /*
4203   Edge effects
4204   //
4205   //
4206   //
4207   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
4208   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
4209   AliXRDPROOFtoolkit tool;
4210   TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
4211   chainTrack->Lookup();
4212   chainTrack->SetProof(kTRUE);
4213
4214   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
4215   chain->Lookup();
4216   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
4217   chainFit->Lookup();
4218   chainFit->SetProof(kTRUE);
4219   chain->SetProof(kTRUE);
4220   //
4221   // Fit cuts
4222   //
4223   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
4224   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
4225   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
4226   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
4227   //
4228   TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
4229   TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
4230   TCut cutN("nclO>20&&nclI>20");
4231   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
4232   //
4233   // Cluster cuts
4234   //
4235   TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
4236   TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
4237   TCut cutClX("abs(Cl.fX)>10");
4238   TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
4239   TCut cutCl=cutClY+cutClZ+cutClX;
4240
4241
4242   // check edge effects
4243   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
4244   //
4245   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
4246
4247   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)
4248
4249
4250
4251   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
4252   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
4253
4254 */
4255