TPC/AliTPCcalibLaser.cxx

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