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