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