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("CalibObjects.root");
  41   TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  42   AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC");
  43   laser->DumpMeanInfo(-0,0)
  44   TFile fmean("laserMean.root")
  45   //
  46   //  laser track clasification;
  47   //
  48   TCut cutT("cutT","abs(Tr.fP[3])<0.06");
  49   TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
  50   TCut cutN("cutN","fTPCncls>70");
  51   TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
  52   TCut cutA = cutT+cutPt+cutP;
  53   TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
  54
  55   //
  56   //
  57   // Analyze  LASER scan
  58   //
  59
  60   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  61   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  62   AliXRDPROOFtoolkit tool;
  63   TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
  64   chain->Lookup();
  65   AliTPCcalibLaser::DumpScanInfo(chain)
  66   TFile fscan("laserScan.root")
  67   TTree * treeT = (TTree*)fscan.Get("Mean")
  68   //
  69   // Analyze laser
  70   //
  71   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  72   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  73   AliXRDPROOFtoolkit tool;
  74   TChain * chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
  75   chainDrift->Lookup();
  76
  77   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
  78   chain->Lookup();
  79   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
  80   chainFit->Lookup();
  81
  82 */
  83
  84
  85
  86 #include "TLinearFitter.h"
  87 #include "AliTPCcalibLaser.h"
  88 #include "AliExternalTrackParam.h"
  89 #include "AliESDEvent.h"
  90 #include "AliESDfriend.h"
  91 #include "AliESDtrack.h"
  92 #include "AliTPCTracklet.h"
  93 #include "TH1D.h"
  94 #include "TH1F.h"
  95 #include "TProfile.h"
  96 #include "TVectorD.h"
  97 #include "TTreeStream.h"
  98 #include "TFile.h"
  99 #include "TF1.h"
 100 #include "TGraphErrors.h"
 101 #include "AliTPCclusterMI.h"
 102 #include "AliTPCseed.h"
 103 #include "AliTracker.h"
 104 #include "AliLog.h"
 105 #include "TClonesArray.h"
 106 #include "TPad.h"
 107 #include "TSystem.h"
 108 #include "TCut.h"
 109 #include "TChain.h"
 110 #include "TH2F.h"
 111 #include "TStatToolkit.h"
 112 #include "TROOT.h"
 113
 114
 115 #include "TTreeStream.h"
 116 #include <iostream>
 117 #include <sstream>
 118 #include "AliTPCLaserTrack.h"
 119 #include "AliTPCcalibDB.h"
 120 #include "AliTPCParam.h"
 121 #include "TTimeStamp.h"
 122 #include "AliDCSSensorArray.h"
 123 #include "AliDCSSensor.h"
 124
 125 using namespace std;
 126
 127 ClassImp(AliTPCcalibLaser)
 128
 129 AliTPCcalibLaser::AliTPCcalibLaser():
 130   AliTPCcalibBase(),
 131   fESD(0),
 132   fESDfriend(0),
 133   fNtracks(0),
 134   fTracksMirror(336),
 135   fTracksEsd(336),
 136   fTracksEsdParam(336),
 137   fTracksTPC(336),
 138   fFullCalib(kTRUE),
 139   fDeltaZ(336),
 140   fDeltaP3(336),
 141   fDeltaP4(336),
 142   fDeltaPhi(336),
 143   fDeltaPhiP(336),
 144   fSignals(336),
 145   //
 146   fHisLaser(0),      //  N dim histogram of laser
 147   fHisNclIn(0),      //->Number of clusters inner
 148   fHisNclOut(0),     //->Number of clusters outer
 149   fHisNclIO(0),      //->Number of cluster inner outer
 150   fHisLclIn(0),      //->Level arm inner
 151   fHisLclOut(0),     //->Level arm outer
 152   fHisLclIO(0),      //->Number of cluster inner outer
 153   fHisdEdx(0),     //->dEdx histo
 154   fHisdZfit(0),     //->distance to the mirror after linear fit
 155   //
 156   //
 157   fHisChi2YIn1(0),      //->chi2 y inner - line
 158   fHisChi2YOut1(0),     //->chi2 y inner - line
 159   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 160   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 161   fHisChi2YIO1(0),      //->chi2 y IO    - common
 162   fHisChi2ZIn1(0),      //->chi2 z inner - line
 163   fHisChi2ZOut1(0),     //->chi2 z inner - line
 164   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 165   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 166   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 167   //
 168   //
 169   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 170   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 171   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 172   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 173   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 174   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 175   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 176   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 177   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 178   //
 179   //
 180   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 181   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 182   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 183   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 184   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 185   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 186   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 187   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 188   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 189   //
 190   fDeltaYres(336),   //->2D histo of residuals
 191   fDeltaZres(336),   //->2D histo fo residuals
 192   fDeltaYres2(336),   //->2D histo of residuals
 193   fDeltaZres2(336),   //->2D histo fo residuals
 194   fDeltaYresAbs(336),   //->2D histo of residuals
 195   fDeltaZresAbs(336),   //->2D histo of residuals
 196   //fDeltaYres3(336),   //->2D histo of residuals
 197   //fDeltaZres3(336),   //->2D histo fo residuals
 198   fFitAside(new TVectorD(5)),
 199   fFitCside(new TVectorD(5)),
 200   fFitACside(new TVectorD(6)),
 201   fEdgeXcuts(3),
 202   fEdgeYcuts(3),
 203   fNClCuts(5),
 204   fNcuts(0),
 205   fBeamSectorOuter(336),
 206   fBeamSectorInner(336),
 207   fBeamOffsetYOuter(336),
 208   fBeamSlopeYOuter(336),
 209   fBeamOffsetYInner(336),
 210   fBeamSlopeYInner(336),
 211   fBeamOffsetZOuter(336),
 212   fBeamSlopeZOuter(336),
 213   fBeamOffsetZInner(336),
 214   fBeamSlopeZInner(336),
 215   fInverseSlopeZ(kTRUE)
 216 {
 217   //
 218   // Constructor
 219   //
 220   fTracksEsdParam.SetOwner(kTRUE);
 221 }
 222
 223 AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
 224   AliTPCcalibBase(),
 225   fESD(0),
 226   fESDfriend(0),
 227   fNtracks(0),
 228   fTracksMirror(336),
 229   fTracksEsd(336),
 230   fTracksEsdParam(336),
 231   fTracksTPC(336),
 232   fFullCalib(full),
 233   //
 234   fDeltaZ(336),          // array of histograms of delta z for each track
 235   fDeltaP3(336),          // array of histograms of delta z for each track
 236   fDeltaP4(336),          // array of histograms of P3 for each track
 237   fDeltaPhi(336),          // array of histograms of P4 for each track
 238   fDeltaPhiP(336),          // array of histograms of delta z for each track
 239   fSignals(336),           // array of dedx signals
 240   //
 241   //
 242   fHisLaser(0),      //  N dim histogram of laser
 243   fHisNclIn(0),      //->Number of clusters inner
 244   fHisNclOut(0),     //->Number of clusters outer
 245   fHisNclIO(0),      //->Number of cluster inner outer
 246   fHisLclIn(0),      //->Level arm inner
 247   fHisLclOut(0),     //->Level arm outer
 248   fHisLclIO(0),      //->Number of cluster inner outer
 249   fHisdEdx(0), //->dEdx histo
 250   fHisdZfit(0),     //->distance to the mirror after linear fit
 251   //
 252   //
 253   fHisChi2YIn1(0),      //->chi2 y inner - line
 254   fHisChi2YOut1(0),     //->chi2 y inner - line
 255   fHisChi2YIn2(0),      //->chi2 y inner - parabola
 256   fHisChi2YOut2(0),     //->chi2 y inner - parabola
 257   fHisChi2YIO1(0),      //->chi2 y IO    - common
 258   fHisChi2ZIn1(0),      //->chi2 z inner - line
 259   fHisChi2ZOut1(0),     //->chi2 z inner - line
 260   fHisChi2ZIn2(0),      //->chi2 z inner - parabola
 261   fHisChi2ZOut2(0),     //->chi2 z inner - parabola
 262   fHisChi2ZIO1(0),      //->chi2 z IO    - common
 263   //
 264   //
 265   fHisPy1vP0(0),     //-> delta y   P0outer-P0inner - line
 266   fHisPy2vP0(0),     //-> delta y   P0outer-P0inner - parabola
 267   fHisPy3vP0(0),     //-> delta y   P0outer-P0inner - common parabola
 268   fHisPy1vP1(0),     //-> delta ky  P1outer-P1inner - line
 269   fHisPy2vP1(0),     //-> delta ky  P1outer-P1inner - parabola
 270   fHisPy3vP1(0),     //-> delta ky  P1outer-P1inner - common parabola
 271   fHisPy2vP2In(0),   //-> Curv  P2inner - parabola
 272   fHisPy2vP2Out(0),  //-> Curv  P2outer - parabola
 273   fHisPy3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 274   //
 275   //
 276   fHisPz1vP0(0),     //-> delta z   P0outer-P0inner - line
 277   fHisPz2vP0(0),     //-> delta z   P0outer-P0inner - parabola
 278   fHisPz3vP0(0),     //-> delta z   P0outer-P0inner - common parabola
 279   fHisPz1vP1(0),     //-> delta kz  P1outer-P1inner - line
 280   fHisPz2vP1(0),     //-> delta kz  P1outer-P1inner - parabola
 281   fHisPz3vP1(0),     //-> delta kz  P1outer-P1inner - common parabola
 282   fHisPz2vP2In(0),   //-> Curv  P2inner - parabola
 283   fHisPz2vP2Out(0),  //-> Curv  P2outer - parabola
 284   fHisPz3vP2IO(0),   //-> Curv  P2outerinner - common parabola
 285   //
 286   //
 287   //
 288   fDeltaYres(336),
 289   fDeltaZres(336),
 290   fDeltaYres2(336),
 291   fDeltaZres2(336),
 292   fDeltaYresAbs(336),
 293   fDeltaZresAbs(336),
 294   //  fDeltaYres3(336),
 295   //fDeltaZres3(336),
 296   fFitAside(new TVectorD(5)),        // drift fit - A side
 297   fFitCside(new TVectorD(5)),        // drift fit - C- side
 298   fFitACside(new TVectorD(6)),        // drift fit - AC- side
 299   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 300   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 301   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 302   fNcuts(0),           // number of cuts
 303   fBeamSectorOuter(336),
 304   fBeamSectorInner(336),
 305   fBeamOffsetYOuter(336),
 306   fBeamSlopeYOuter(336),
 307   fBeamOffsetYInner(336),
 308   fBeamSlopeYInner(336),
 309   fBeamOffsetZOuter(336),
 310   fBeamSlopeZOuter(336),
 311   fBeamOffsetZInner(336),
 312   fBeamSlopeZInner(336),
 313   fInverseSlopeZ(kTRUE)
 314 {
 315   SetName(name);
 316   SetTitle(title);
 317   //
 318   // Constructor
 319   //
 320   fTracksEsdParam.SetOwner(kTRUE);
 321 }
 322
 323 AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
 324   AliTPCcalibBase(calibLaser),
 325   fESD(0),
 326   fESDfriend(0),
 327   fNtracks(0),
 328   fTracksMirror(336),
 329   fTracksEsd(336),
 330   fTracksEsdParam(336),
 331   fTracksTPC(336),
 332   fFullCalib(calibLaser.fFullCalib),
 333   //
 334   fDeltaZ(calibLaser.fDeltaZ),          // array of histograms of delta z for each track
 335   fDeltaP3(((calibLaser.fDeltaP3))),          // array of histograms of delta z for each track
 336   fDeltaP4(((calibLaser.fDeltaP4))),          // array of histograms of P3 for each track
 337   fDeltaPhi(((calibLaser.fDeltaPhi))),          // array of histograms of P4 for each track
 338   fDeltaPhiP(((calibLaser.fDeltaPhiP))),          // array of histograms of delta z for each track
 339   fSignals(((calibLaser.fSignals))),           // array of dedx signals
 340   //
 341   //
 342   fHisLaser(0),      //  N dim histogram of laser
 343   fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))),      //->Number of clusters inner
 344   fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))),     //->Number of clusters outer
 345   fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))),      //->Number of cluster inner outer
 346   fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))),      //->Level arm inner
 347   fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))),     //->Level arm outer
 348   fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))),      //->Number of cluster inner outer
 349   fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))),      //
 350   fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))),     //->distance to the mirror after linear fit
 351   //
 352   //
 353   fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))),      //->chi2 y inner - line
 354   fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))),     //->chi2 y inner - line
 355   fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))),      //->chi2 y inner - parabola
 356   fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))),     //->chi2 y inner - parabola
 357   fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))),      //->chi2 y IO    - common
 358   fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))),      //->chi2 z inner - line
 359   fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))),     //->chi2 z inner - line
 360   fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))),      //->chi2 z inner - parabola
 361   fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))),     //->chi2 z inner - parabola
 362   fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))),      //->chi2 z IO    - common
 363   //
 364   //
 365   fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))),     //-> delta y   P0outer-P0inner - line
 366   fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))),     //-> delta y   P0outer-P0inner - parabola
 367   fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))),     //-> delta y   P0outer-P0inner - common parabola
 368   fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))),     //-> delta ky  P1outer-P1inner - line
 369   fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))),     //-> delta ky  P1outer-P1inner - parabola
 370   fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))),     //-> delta ky  P1outer-P1inner - common parabola
 371   fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))),   //-> Curv  P2inner - parabola
 372   fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))),  //-> Curv  P2outer - parabola
 373   fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))),   //-> Curv  P2outerinner - common parabola
 374   //
 375   //
 376   fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))),     //-> delta z   P0outer-P0inner - line
 377   fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))),     //-> delta z   P0outer-P0inner - parabola
 378   fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))),     //-> delta z   P0outer-P0inner - common parabola
 379   fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))),     //-> delta kz  P1outer-P1inner - line
 380   fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))),     //-> delta kz  P1outer-P1inner - parabola
 381   fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))),     //-> delta kz  P1outer-P1inner - common parabola
 382   fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))),   //-> Curv  P2inner - parabola
 383   fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))),  //-> Curv  P2outer - parabola
 384   fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))),   //-> Curv  P2outerinner - common parabola
 385   //
 386   //
 387   fDeltaYres(((calibLaser.fDeltaYres))),
 388   fDeltaZres(((calibLaser.fDeltaZres))),
 389   fDeltaYres2(((calibLaser.fDeltaYres))),
 390   fDeltaZres2(((calibLaser.fDeltaZres))),
 391   fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
 392   fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
 393   //  fDeltaYres3(((calibLaser.fDeltaYres))),
 394   //fDeltaZres3(((calibLaser.fDeltaZres))),
 395   fFitAside(new TVectorD(5)),        // drift fit - A side
 396   fFitCside(new TVectorD(5)),        // drift fit - C- side
 397   fFitACside(new TVectorD(6)),        // drift fit - C- side
 398   fEdgeXcuts(3),       // cuts in local x direction; used in the refit of the laser tracks
 399   fEdgeYcuts(3),       // cuts in local y direction; used in the refit of the laser tracks
 400   fNClCuts(5),         // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
 401   fNcuts(0),           // number of cuts
 402   fBeamSectorOuter(336),
 403   fBeamSectorInner(336),
 404   fBeamOffsetYOuter(336),
 405   fBeamSlopeYOuter(336),
 406   fBeamOffsetYInner(336),
 407   fBeamSlopeYInner(336),
 408   fBeamOffsetZOuter(336),
 409   fBeamSlopeZOuter(336),
 410   fBeamOffsetZInner(336),
 411   fBeamSlopeZInner(336),
 412   fInverseSlopeZ(calibLaser.fInverseSlopeZ)
 413
 414 {
 415   //
 416   // copy constructor
 417   //
 418 }
 419
 420
 421
 422 AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
 423   //
 424   // assgnment operator
 425   //
 426   if (this != &calibLaser) {
 427     new (this) AliTPCcalibLaser(calibLaser);
 428   }
 429   return *this;
 430
 431 }
 432
 433
 434
 435
 436 AliTPCcalibLaser::~AliTPCcalibLaser() {
 437   //
 438   // destructor
 439   //
 440   if ( fHisNclIn){
 441     delete fHisLaser;      //->
 442     delete fHisNclIn;      //->Number of clusters inner
 443     delete fHisNclOut;     //->Number of clusters outer
 444     delete fHisNclIO;      //->Number of cluster inner outer
 445     delete fHisLclIn;      //->Level arm inner
 446     delete fHisLclOut;     //->Level arm outer
 447     delete fHisLclIO;      //->Number of cluster inner outer
 448     delete fHisdEdx;
 449     delete fHisdZfit;
 450     //
 451     //
 452     delete fHisChi2YIn1;      //->chi2 y inner - line
 453     delete fHisChi2YOut1;     //->chi2 y inner - line
 454     delete fHisChi2YIn2;      //->chi2 y inner - parabola
 455     delete fHisChi2YOut2;     //->chi2 y inner - parabola
 456     delete fHisChi2YIO1;      //->chi2 y IO    - common
 457     delete fHisChi2ZIn1;      //->chi2 z inner - line
 458     delete fHisChi2ZOut1;     //->chi2 z inner - line
 459     delete fHisChi2ZIn2;      //->chi2 z inner - parabola
 460     delete fHisChi2ZOut2;     //->chi2 z inner - parabola
 461     delete fHisChi2ZIO1;      //->chi2 z IO    - common
 462     //
 463     //
 464     delete fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
 465     delete fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
 466     delete fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
 467     delete fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
 468     delete fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
 469     delete fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
 470     delete fHisPy2vP2In;   //-> Curv  P2inner - parabola
 471     delete fHisPy2vP2Out;  //-> Curv  P2outer - parabola
 472     delete fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
 473     //
 474     //
 475     delete fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
 476     delete fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
 477     delete fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
 478     delete fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
 479     delete fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
 480     delete fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
 481     delete fHisPz2vP2In;   //-> Curv  P2inner - parabola
 482     delete fHisPz2vP2Out;  //-> Curv  P2outer - parabola
 483     delete fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
 484   }
 485   //
 486   //
 487   //
 488   fDeltaZ.SetOwner();          //-> array of histograms of delta z for each track
 489   fDeltaP3.SetOwner();         //-> array of histograms of P3      for each track
 490   fDeltaP4.SetOwner();         //-> array of histograms of P4      for each track
 491   fDeltaPhi.SetOwner();        //-> array of histograms of delta z for each track
 492   fDeltaPhiP.SetOwner();       //-> array of histograms of delta z for each track
 493   fSignals.SetOwner();         //->Array of dedx signals
 494
 495   fDeltaZ.Delete();          //-> array of histograms of delta z for each track
 496   fDeltaP3.Delete();         //-> array of histograms of P3      for each track
 497   fDeltaP4.Delete();         //-> array of histograms of P4      for each track
 498   fDeltaPhi.Delete();        //-> array of histograms of delta z for each track
 499   fDeltaPhiP.Delete();       //-> array of histograms of delta z for each track
 500   fSignals.Delete();         //->Array of dedx signals
 501
 502   fDeltaYres.SetOwner();
 503   fDeltaYres.Delete();
 504   fDeltaZres.SetOwner();
 505   fDeltaZres.Delete();
 506   fDeltaYres2.SetOwner();
 507   fDeltaYres2.Delete();
 508   fDeltaZres2.SetOwner();
 509   fDeltaZres2.Delete();
 510
 511   fDeltaYresAbs.SetOwner();
 512   fDeltaYresAbs.Delete();
 513   fDeltaZresAbs.SetOwner();
 514   fDeltaZresAbs.Delete();
 515 }
 516
 517
 518
 519 void AliTPCcalibLaser::Process(AliESDEvent * event) {
 520   //
 521   //
 522   // Loop over tracks and call  Process function
 523   //
 524   Int_t kMinTracks=20;
 525   fESD = event;
 526   if (!fESD) {
 527     return;
 528   }
 529   fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
 530   if (!fESDfriend) {
 531     return;
 532   }
 533   if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
 534   AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
 535   fTracksTPC.Clear();
 536   fTracksEsd.Clear();
 537   fTracksEsdParam.Delete();
 538   for (Int_t id=0; id<336;id++) {
 539     fCounter[id]=0;
 540     fClusterCounter[id]=0;
 541     fClusterSatur[id]=0;
 542   }
 543   //
 544   Int_t n=fESD->GetNumberOfTracks();
 545   Int_t counter=0;
 546   for (Int_t i=0;i<n;++i) {
 547     AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
 548     if (!friendTrack) continue;
 549     AliESDtrack *track=fESD->GetTrack(i);
 550     TObject *calibObject=0;
 551     AliTPCseed *seed=0;
 552     for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
 553       if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
 554         break;
 555     if (track&&seed &&TMath::Abs(track->Pt()) >1 ) {
 556       //filter CE tracks
 557       Int_t id = FindMirror(track,seed);
 558       if (id>0) counter++;
 559     }
 560     //
 561   }
 562   fNtracks=counter;
 563   if (counter<kMinTracks) return;
 564
 565   FitDriftV();
 566   if (!fFullCalib) return;
 567   static Bool_t init=kFALSE;
 568   if (!init){
 569     init = kTRUE;  // way around for PROOF - to be investigated
 570     MakeFitHistos();
 571   }
 572   //
 573   for (Int_t id=0; id<336; id++){
 574     //
 575     if (!fTracksEsdParam.At(id)) continue;
 576     if (fClusterSatur[id]>0.3)   continue;  // tracks in saturation
 577     DumpLaser(id);
 578     if ( AcceptLaser(id) && fFitZ[id]<0.5){
 579       RefitLaserJW(id);
 580       MakeDistHisto(id);
 581     }
 582   }
 583
 584 }
 585
 586 void AliTPCcalibLaser::MakeDistHisto(Int_t id){
 587   //
 588   //
 589   //
 590   //  for (Int_t id=0; id<336; id++){
 591     //
 592     //
 593     if (!fTracksEsdParam.At(id)) return;
 594     if (!AcceptLaser(id)) return;
 595     Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
 596     //
 597     //
 598     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
 599     if (!hisdz) MakeFitHistos();
 600     hisdz = (TH1F*)fDeltaZ.At(id);
 601     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
 602     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
 603     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
 604     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
 605     TH1F * hisSignal = (TH1F*)fSignals.At(id);
 606     //
 607
 608     AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 609     AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 610     AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 611     if (!param) return;
 612     if (!ltrp) return;
 613     if (!track) return;
 614     Double_t xyz[3];
 615     Double_t pxyz[3];
 616     Double_t lxyz[3];
 617     Double_t lpxyz[3];
 618     param->GetXYZ(xyz);
 619     param->GetPxPyPz(pxyz);
 620     ltrp->GetXYZ(lxyz);
 621     ltrp->GetPxPyPz(lpxyz);
 622     //
 623     Float_t dz   = param->GetZ()-ltrp->GetZ();
 624     Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
 625     Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
 626     if (hisdz) hisdz->Fill(dz);
 627     if (hisP3) hisP3->Fill(param->GetParameter()[3]);
 628     if (hisP4) hisP4->Fill(param->GetParameter()[4]);
 629     if (hisdphi) hisdphi->Fill(dphi);
 630     if (hisdphiP) hisdphiP->Fill(dphiP);
 631     if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
 632     // fill HisLaser
 633     xhis[0] = ltrp->GetId();
 634     xhis[1] = ltrp->GetSide();
 635     xhis[2] = ltrp->GetRod();
 636     xhis[3] = ltrp->GetBundle();
 637     xhis[4] = ltrp->GetBeam();
 638     xhis[5] = dphi;
 639     xhis[6] = fFitZ[id];
 640     xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
 641     xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
 642     xhis[9] = param->GetParameter()[4];
 643     xhis[10]= track->GetTPCNcls();
 644     xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
 645     // }
 646     fHisLaser->Fill(xhis);
 647 }
 648
 649 void AliTPCcalibLaser::FitDriftV(){
 650   //
 651   // Fit corrections to the drift velocity - linear approximation in the z and global y
 652   //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
 653   //
 654   /*
 655     Formulas:
 656
 657     z  = s* (z0 - vd*(t-t0))
 658
 659     s  - side -1 and +1
 660     t0 - time 0
 661     vd - nominal drift velocity
 662     zs - miscalibrated position
 663
 664     zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
 665     vr  - relative change of the drift velocity
 666     dzt - vd*dt
 667     dr  = zz0-s*z
 668     ..
 669     ==>
 670     zs ~ z - s*vr*(z0-s*z)+s*dzt
 671     --------------------------------
 672     1. Correction function vr constant:
 673
 674
 675     dz = zs-z = -s*vr *(z0-s*z)+s*dzt
 676     dzs/dl = dz/dl +s*s*vr*dz/dl
 677     d(dz/dl) = vr*dz/dl
 678   */
 679
 680   const Float_t kZCut        = 240;    // remove the closest laser beam
 681   const Float_t kSaturCut    = 0.05;   // remove saturated lasers - cut on fraction of saturated
 682   const Float_t kDistCut     = 3;    // distance sigma cut
 683   const Float_t kDistCutAbs  = 0.25;
 684   const Float_t kMinClusters = 60;  // minimal amount of the clusters
 685   const Float_t kMinSignal   = 16;  // minimal mean height of the signal
 686   const Float_t kChi2Cut     = 0.1; // chi2 cut to accept drift fit
 687   static TLinearFitter fdriftA(3,"hyp2");
 688   static TLinearFitter fdriftC(3,"hyp2");
 689   static TLinearFitter fdriftAC(4,"hyp3");
 690   TVectorD fitA(3),fitC(3),fitAC(4);
 691
 692   AliTPCcalibDB*  calib=AliTPCcalibDB::Instance();
 693   AliTPCParam  * tpcparam    = calib->GetParameters();
 694
 695   //
 696   for (Int_t id=0; id<336; id++) fFitZ[id]=0;
 697
 698   Float_t chi2A = 10;
 699   Float_t chi2C = 10;
 700   Float_t chi2AC = 10;
 701   Int_t npointsA=0;
 702   Int_t npointsC=0;
 703   Int_t npointsAC=0;
 704
 705
 706   for (Int_t iter=0; iter<3; iter++){
 707     fdriftA.ClearPoints();
 708     fdriftC.ClearPoints();
 709     fdriftAC.ClearPoints();
 710     //
 711     for (Int_t id=0; id<336; id++){
 712       if (!fTracksEsdParam.At(id)) continue;
 713       if (!AcceptLaser(id)) continue;
 714       if ( fClusterSatur[id]>kSaturCut)  continue;
 715       if ( fClusterCounter[id]<kMinClusters)  continue;
 716       AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 717       if (track->GetTPCsignal()<kMinSignal) continue;
 718       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 719       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 720
 721       Double_t xyz[3];
 722       Double_t pxyz[3];
 723       Double_t lxyz[3];
 724       Double_t lpxyz[3];
 725       param->GetXYZ(xyz);
 726       param->GetPxPyPz(pxyz);
 727       ltrp->GetXYZ(lxyz);
 728       ltrp->GetPxPyPz(lpxyz);
 729       if (TMath::Abs(lxyz[2])>kZCut) continue;
 730       Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
 731       if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
 732       if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
 733       if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
 734
 735       // drift distance
 736       Double_t zlength =  tpcparam->GetZLength(0);
 737       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
 738       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
 739       Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
 740       if (ltrp->GetSide()==0){
 741         fdriftA.AddPoint(xxx,mdrift,1);
 742       }else{
 743         fdriftC.AddPoint(xxx,mdrift,1);
 744       }
 745       Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
 746       fdriftAC.AddPoint(xxx2,mdrift,1);
 747     }
 748     //
 749     if (fdriftA.GetNpoints()>10){
 750       //
 751       fdriftA.Eval();
 752       if (iter==0) fdriftA.EvalRobust(0.7);
 753       else fdriftA.EvalRobust(0.8);
 754       fdriftA.GetParameters(fitA);
 755       npointsA= fdriftA.GetNpoints();
 756       chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
 757       if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
 758         if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
 759         (*fFitAside)[0] = fitA[0];
 760         (*fFitAside)[1] = fitA[1];
 761         (*fFitAside)[2] = fitA[2];
 762         (*fFitAside)[3] = fdriftA.GetNpoints();
 763         (*fFitAside)[4] = chi2A;
 764       }
 765     }
 766     if (fdriftC.GetNpoints()>10){
 767       fdriftC.Eval();
 768       if (iter==0) fdriftC.EvalRobust(0.7);
 769       else fdriftC.EvalRobust(0.8);
 770       //
 771       fdriftC.GetParameters(fitC);
 772       npointsC= fdriftC.GetNpoints();
 773       chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
 774       if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
 775         if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
 776         (*fFitCside)[0] = fitC[0];
 777         (*fFitCside)[1] = fitC[1];
 778         (*fFitCside)[2] = fitC[2];
 779         (*fFitCside)[3] = fdriftC.GetNpoints();
 780         (*fFitCside)[4] = chi2C;
 781       }
 782     }
 783
 784     if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
 785       fdriftAC.Eval();
 786       if (iter==0) fdriftAC.EvalRobust(0.7);
 787       else fdriftAC.EvalRobust(0.8);
 788       //
 789       fdriftAC.GetParameters(fitAC);
 790       npointsAC= fdriftAC.GetNpoints();
 791       chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
 792       if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
 793     }
 794
 795     for (Int_t id=0; id<336; id++){
 796       if (!fTracksEsdParam.At(id)) continue;
 797       //
 798       AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
 799       AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 800       Double_t xyz[3];
 801       Double_t pxyz[3];
 802       Double_t lxyz[3];
 803       Double_t lpxyz[3];
 804       param->GetXYZ(xyz);
 805       param->GetPxPyPz(pxyz);
 806       ltrp->GetXYZ(lxyz);
 807       ltrp->GetPxPyPz(lpxyz);
 808       Double_t zlength =  tpcparam->GetZLength(0);
 809       Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
 810       Double_t mdrift = zlength-TMath::Abs(xyz[2]);
 811
 812       Float_t fz =0;
 813       if (ltrp->GetSide()==0){
 814         fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
 815       }else{
 816         fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
 817       }
 818       if (npointsAC>10){
 819         fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
 820       }
 821       fFitZ[id]=mdrift-fz;
 822     }
 823     if (fStreamLevel>0){
 824       TTreeSRedirector *cstream = GetDebugStreamer();
 825       TTimeStamp tstamp(fTime);
 826       Float_t valuePressure0  = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
 827       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
 828       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
 829       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
 830       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
 831       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
 832       TVectorD vecGoofie(20);
 833       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
 834       if (goofieArray)
 835         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
 836           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
 837           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
 838         }
 839
 840       if (cstream){
 841         (*cstream)<<"driftv"<<
 842           "run="<<fRun<<              //  run number
 843           "event="<<fEvent<<          //  event number
 844           "time="<<fTime<<            //  time stamp of event
 845           "trigger="<<fTrigger<<      //  trigger
 846           "mag="<<fMagF<<             //  magnetic field
 847           // Environment values
 848           "press0="<<valuePressure0<<
 849           "press1="<<valuePressure1<<
 850           "pt0="<<ptrelative0<<
 851           "pt1="<<ptrelative1<<
 852           "temp0="<<temp0<<
 853           "temp1="<<temp1<<
 854           "vecGoofie.="<<&vecGoofie<<
 855           //
 856           //
 857           "iter="<<iter<<
 858           "driftA.="<<fFitAside<<
 859           "driftC.="<<fFitCside<<
 860           "driftAC.="<<fFitACside<<
 861           "chi2A="<<chi2A<<
 862           "chi2C="<<chi2C<<
 863           "chi2AC="<<chi2AC<<
 864           "nA="<<npointsA<<
 865           "nC="<<npointsC<<
 866           "nAC="<<npointsAC<<
 867           "\n";
 868       }
 869     }
 870   }
 871 }
 872
 873 Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
 874   //
 875   // get distance between mirror and laser track
 876   //
 877   //
 878   Double_t xyz[3];
 879   Double_t lxyz[3];
 880   param->GetXYZ(xyz);
 881   ltrp->GetXYZ(lxyz);
 882   //
 883   //
 884   Double_t dist = 0;
 885   //radial distance
 886   dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
 887   //
 888   // z distance
 889   // apply drift correction if already exist
 890   //
 891   Float_t dz = 0;
 892   if (ltrp->GetSide()==0){
 893     if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
 894   }else{
 895     if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
 896   }
 897   if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
 898   dist+=TMath::Abs(dz);
 899   //
 900   // phi dist - divergence on 50 cm
 901   //
 902   dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
 903   return dist;
 904 }
 905
 906
 907 Bool_t  AliTPCcalibLaser::AcceptLaser(Int_t id){
 908   //
 909   //
 910   //
 911   /*
 912   TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
 913   TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
 914   TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
 915   TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
 916   TCut cutP4("cutPt","abs(Tr.fP[4])<0.1");
 917
 918   TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
 919   */
 920   AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
 921   AliTPCLaserTrack *ltrp  = ( AliTPCLaserTrack*)fTracksMirror.At(id);
 922   //AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
 923   Double_t xyz[3];
 924   Double_t lxyz[3];
 925   param->GetXYZ(xyz);
 926   ltrp->GetXYZ(lxyz);
 927   if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
 928   if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE;    // cutZ -P1
 929   if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE;  // cut -P2
 930   if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE;   // cut Tl -P3
 931   if (TMath::Abs(param->GetParameter()[4])>0.1) return kFALSE;   // cut Pt  -P4
 932   //
 933   //
 934
 935   return kTRUE;
 936 }
 937
 938 Int_t  AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
 939   //
 940   // Find corresponding mirror
 941   // add the corresponding tracks
 942
 943
 944   Float_t kRadius0  = 252;
 945   Float_t kRadius   = 253.4;
 946
 947   if (!track->GetOuterParam()) return -1;
 948   AliExternalTrackParam param(*(track->GetOuterParam()));
 949   AliTracker::PropagateTrackTo(&param,kRadius0,0.10566,3,kTRUE);
 950   AliTracker::PropagateTrackTo(&param,kRadius,0.10566,0.1,kTRUE);
 951   AliTPCLaserTrack ltr;
 952   AliTPCLaserTrack *ltrp=0x0;
 953   //  AliTPCLaserTrack *ltrpjw=0x0;
 954   //
 955   Int_t id   = AliTPCLaserTrack::IdentifyTrack(&param);
 956  // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(&param);
 957   //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
 958
 959   if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
 960     ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
 961   else
 962     ltrp=&ltr;
 963
 964   if (id<0) return id;
 965   fCounter[id]++;
 966   //
 967   //
 968   //
 969   Int_t countercl=0;
 970   Float_t counterSatur=0;
 971   for (Int_t irow=158;irow>-1;--irow) {
 972     AliTPCclusterMI *c=seed->GetClusterPointer(irow);
 973     if (!c) continue;
 974     Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
 975     Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
 976     if (pedgeY<3) continue;
 977     if (pedgeX<3) continue;
 978     countercl++;
 979     if (c->GetMax()>900) counterSatur++;
 980   }
 981   counterSatur/=(countercl+1);
 982   //
 983   if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
 984   //
 985   //
 986   Float_t radius=TMath::Abs(ltrp->GetX());
 987   AliTracker::PropagateTrackTo(&param,radius,0.10566,0.01,kTRUE);
 988   param.Rotate(ltrp->GetAlpha());
 989   //
 990   if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
 991   Bool_t accept=kTRUE;
 992   //
 993   // choose closer track
 994   //
 995   AliExternalTrackParam * param0  = (AliExternalTrackParam *)fTracksEsdParam.At(id);
 996   if (param0){
 997     Float_t dist0=GetDistance(param0,ltrp);
 998     Float_t dist1=GetDistance(&param,ltrp);
 999     if (dist0<dist1)    accept=kFALSE;
1000   }
1001
1002   if (accept){
1003     fClusterCounter[id]=countercl;
1004     fTracksEsdParam.AddAt(param.Clone(),id);
1005     fTracksEsd.AddAt(track,id);
1006     fTracksTPC.AddAt(seed,id);
1007   }
1008   return id;
1009 }
1010
1011
1012
1013 void AliTPCcalibLaser::DumpLaser(Int_t id) {
1014   //
1015   //  Dump Laser info to the tree
1016   //
1017   AliESDtrack   *track    = (AliESDtrack*)fTracksEsd.At(id);
1018   AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1019   AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
1020   //
1021   // Fast laser ID
1022   //
1023   Double_t xyz[3];
1024   Double_t pxyz[3];
1025   Double_t lxyz[3];
1026   Double_t lpxyz[3];
1027   param->GetXYZ(xyz);
1028   param->GetPxPyPz(pxyz);
1029   ltrp->GetXYZ(lxyz);
1030   ltrp->GetPxPyPz(lpxyz);
1031   Float_t dist3D   = GetDistance(param,ltrp);
1032   Float_t dist0    = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
1033   Float_t distphi  = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
1034
1035
1036   if (fStreamLevel>0){
1037     TTreeSRedirector *cstream = GetDebugStreamer();
1038     Int_t time = fESD->GetTimeStamp();
1039     Bool_t accept = AcceptLaser(id);
1040     if (cstream){
1041       (*cstream)<<"Track"<<
1042         //
1043         "run="<<fRun<<              //  run number
1044         "event="<<fEvent<<          //  event number
1045         "time="<<fTime<<            //  time stamp of event
1046         "trigger="<<fTrigger<<      //  trigger
1047         "mag="<<fMagF<<             //  magnetic field
1048
1049         "id="<<id<<
1050         "accept="<<accept<<
1051         "driftA.="<<fFitAside<<
1052         "driftC.="<<fFitCside<<
1053         "time="<<time<<
1054         "dist3D="<<dist3D<<
1055         "dist0="<<dist0<<
1056         "distphi="<<distphi<<
1057         //
1058         //
1059         "counter="<<fCounter[id]<<
1060         "clcounter="<<fClusterCounter[id]<<
1061         "satur="<<fClusterSatur[id]<<
1062         "fz="<<fFitZ[id]<<
1063         //
1064         "LTr.="<<ltrp<<
1065         "Esd.="<<track<<
1066         "Tr.="<<param<<
1067         "x0="<<xyz[0]<<
1068         "x1="<<xyz[1]<<
1069         "x2="<<xyz[2]<<
1070         "px0="<<pxyz[0]<<
1071         "px1="<<pxyz[1]<<
1072         "px2="<<pxyz[2]<<
1073         //
1074         "lx0="<<lxyz[0]<<
1075         "lx1="<<lxyz[1]<<
1076         "lx2="<<lxyz[2]<<
1077         "lpx0="<<lpxyz[0]<<
1078         "lpx1="<<lpxyz[1]<<
1079         "lpx2="<<lpxyz[2]<<
1080         "\n";
1081     }
1082   }
1083 }
1084
1085 void AliTPCcalibLaser::RefitLaserJW(Int_t id){
1086   //
1087   // Refit the track with different tracklet models:
1088   // 1. Per ROC using the kalman filter, different edge cuts
1089   // 2. Per ROC linear in y and z
1090   // 3. Per ROC quadratic in y and z
1091   // 4. Per track offset for each sector, linear for each sector, common quadratic
1092   // store x, y, z information for all models and the cluster to calculate the residuals
1093   //
1094
1095   // The clusters which do not fulfill given  criteria are skipped
1096   //
1097   // Cluters removed from fit
1098   //  1. Extended shape     > kShapeCut
1099   //  2. In saturation  Max > kMax
1100   //  3. Distance to edge   < cutEdge
1101   //
1102   // Clusters not used for the calibration:
1103   //
1104   //  1. Extended shape     > kShapeCut
1105   //  2. In saturation  Max > kMax
1106
1107
1108   AliTPCseed *track      = (AliTPCseed*)fTracksTPC.At(id);
1109   AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
1110   AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1111
1112   AliTPCclusterMI dummyCl;
1113
1114   //two tracklets
1115   Int_t kMaxTracklets=2;
1116   Float_t kShapeCut = 1.3;
1117   Float_t kRatioCut = 2.;
1118
1119   Float_t kMax      = 900.;
1120
1121
1122   //=============================================//
1123   // Linear Fitters for the Different Approaches //
1124   //=============================================//
1125   //linear fit model in y and z; inner - outer sector, combined with offset
1126   static TLinearFitter fy1I(2,"hyp1");
1127   static TLinearFitter fy1O(2,"hyp1");
1128   static TLinearFitter fz1I(2,"hyp1");
1129   static TLinearFitter fz1O(2,"hyp1");
1130   static TLinearFitter fy1IO(3,"hyp2");
1131   static TLinearFitter fz1IO(3,"hyp2");
1132   //quadratic fit model in y and z; inner - sector
1133   static TLinearFitter fy2I(3,"hyp2");
1134   static TLinearFitter fy2O(3,"hyp2");
1135   static TLinearFitter fz2I(3,"hyp2");
1136   static TLinearFitter fz2O(3,"hyp2");
1137   //common quadratic fit for IROC and OROC in y and z
1138   static TLinearFitter fy4(5,"hyp4");
1139   static TLinearFitter fz4(5,"hyp4");
1140
1141
1142   //set standard cuts
1143   if ( fNcuts==0 ){
1144       fNcuts=1;
1145       fEdgeXcuts[0]=4;
1146       fEdgeYcuts[0]=3;
1147       fNClCuts[0]=20;
1148   }
1149   //=============================//
1150   // Loop over all Tracklet Cuts //
1151   //=============================//
1152   for (Int_t icut=0; icut<fNcuts; icut++){
1153     Float_t xinMin = 2500, xinMax=-90;
1154     Float_t xoutMin=2500, xoutMax=-90;
1155     Float_t msigmaYIn=0;
1156     Float_t msigmaYOut=0;
1157     Float_t mqratioIn=0;
1158     Float_t mqratioOut=0;
1159
1160
1161       AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
1162       //cut parameters
1163       Double_t edgeCutX = fEdgeXcuts[icut];
1164       Double_t edgeCutY = fEdgeYcuts[icut];
1165       Int_t    nclCut   = (Int_t)fNClCuts[icut];
1166       //=========================//
1167       // Parameters to calculate //
1168       //=========================//
1169       //fit parameter inner, outer tracklet and combined fit
1170       TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
1171       TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
1172       //
1173       TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
1174       TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
1175       TVectorD vecy4res(5),vecz4res(5);
1176       TVectorD vecy1resIO(3),vecz1resIO(3);
1177       // cluster and track positions for each row - used for residuals
1178       TVectorD vecgX(159);        // global X
1179       TVectorD vecgY(159);        // global Y
1180       TVectorD vecgZ(159);        // global Z
1181
1182       TVectorD vecX(159);        // x is the same for all (row center)
1183       TVectorD vecYkalman(159);  // y from kalman fit
1184       TVectorD vecZkalman(159);  // z from kalman fit
1185       TVectorD vecY1(159);       // y from pol1 fit per ROC
1186       TVectorD vecZ1(159);       // z from pol1 fit per ROC
1187       TVectorD vecY1IO(159);     // y from pol1 fit per ROC
1188       TVectorD vecZ1IO(159);     // z from pol1 fit per ROC
1189       TVectorD vecY2(159);       // y from pol2 fit per ROC
1190       TVectorD vecZ2(159);       // z from pol2 fit per ROC
1191       TVectorD vecY4(159);       // y from sector fit
1192       TVectorD vecZ4(159);       // z from sector fit
1193       TVectorD vecClY(159);      // y cluster position
1194       TVectorD vecClZ(159);      // z cluster position
1195       TVectorD vecSec(159);      // sector for each row
1196       TVectorD  isReject(159);    // flag - cluster to be rejected
1197       //chi2 of fits
1198       Double_t chi2I1z=0;       // chi2 of pol1 fit in z (inner)
1199       Double_t chi2I1y=0;       // chi2 of pol1 fit in y (inner)
1200       Double_t chi2O1z=0;       // chi2 of pol1 fit in z (outer)
1201       Double_t chi2O1y=0;       // chi2 of pol1 fit in y (outer)
1202       Double_t chi2IO1z=0;       // chi2 of pol1 fit in z (outer)
1203       Double_t chi2IO1y=0;       // chi2 of pol1 fit in y (outer)
1204       Double_t chi2I2z=0;       // chi2 of pol2 fit in z (inner)
1205       Double_t chi2I2y=0;       // chi2 of pol2 fit in y (inner)
1206       Double_t chi2O2z=0;       // chi2 of pol2 fit in z (outer)
1207       Double_t chi2O2y=0;       // chi2 of pol2 fit in y (outer)
1208       Double_t chi2IOz=0;       // chi2 of hyp4 fit in z (inner+outer)
1209       Double_t chi2IOy=0;       // chi2 of hyp4 fit in y (inner+outer)
1210       //more
1211       Int_t innerSector = -1;    // number of inner sector
1212       Int_t outerSector = -1;    // number of outer sector
1213       Int_t nclI=0;              // number of clusters (inner)
1214       Int_t nclO=0;              // number of clusters (outer)
1215       //=================================================//
1216       // Perform the Kalman Fit using the Tracklet Class //
1217       //=================================================//
1218       AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
1219       TObjArray tracklets=
1220           AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
1221                                           kFALSE,nclCut,kMaxTracklets);
1222       // tracklet pointers
1223       AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
1224       AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
1225       AliTPCTracklet *tr=0x0;
1226       AliTPCTracklet dummy;
1227       //continue if we didn't find a tracklet
1228       if ( !trInner && !trOuter ) continue;
1229       //================================================================================//
1230       // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
1231       //================================================================================//
1232       if ( trInner && trOuter ){
1233           if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
1234           if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
1235               tr = trInner;
1236               trInner=trOuter;
1237               trOuter=tr;
1238               AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
1239           }
1240       } else {
1241           if ( trInner ){
1242               if ( !trInner->GetInner() ) continue;
1243               trOuter=&dummy;
1244               if ( trInner->GetSector()>35 ){
1245                   trOuter=trInner;
1246                   trInner=&dummy;
1247               }
1248           } else { //trOuter
1249               if ( !trOuter->GetInner() ) continue;
1250               trInner=&dummy;
1251               if ( trOuter->GetSector()<36 ){
1252                   trInner=trOuter;
1253                   trOuter=&dummy;
1254               }
1255           }
1256       }
1257       innerSector = trInner->GetSector();
1258       if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
1259       outerSector = trOuter->GetSector();
1260       if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
1261
1262       // array of clusters
1263       TClonesArray arrCl("AliTPCclusterMI",159);
1264       arrCl.ExpandCreateFast(159);
1265       //=======================================//
1266       // fill fitters with cluster information //
1267       //=======================================//
1268       AliDebug(3,"Filing Cluster Information");
1269
1270       //
1271       //
1272
1273       for (Int_t irow=158;irow>-1;--irow) {
1274           AliTPCclusterMI *c=track->GetClusterPointer(irow);
1275           AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
1276           cl=dummyCl;
1277           vecX[irow]   = 0;
1278           vecClY[irow] = 0;
1279           vecClZ[irow] = 0;
1280           Float_t meanY=0, sumY=0;
1281           for (Int_t drow=-1;drow<=1;drow++) {
1282             if (irow+drow<0) continue;
1283             if (irow+drow>158) continue;
1284             AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
1285             if (!ccurrent) continue;
1286             Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
1287             if ( roc!=innerSector && roc!=outerSector ) continue;
1288             if (ccurrent->GetX()<10) continue;
1289             if (ccurrent->GetY()==0) continue;
1290             meanY+=ccurrent->GetY();
1291             sumY++;
1292           }
1293           if (sumY>0)  meanY/=sumY;
1294
1295           //
1296           vecSec[irow]=-1;
1297           if (!c) continue;
1298           Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
1299           Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
1300
1301           //
1302           Int_t roc = static_cast<Int_t>(c->GetDetector());
1303           if ( roc!=innerSector && roc!=outerSector ) continue;
1304           vecSec[irow]=roc;
1305           //store clusters in clones array
1306           cl=*c;
1307           //
1308           if (c->GetMax()<4)   continue;  // noise cluster?
1309           if (TMath::Abs(c->GetY())<0.0001)   continue;  // noise cluster?
1310           //cluster position
1311           vecX[irow]   = c->GetX();
1312           vecClY[irow] = c->GetY();
1313           vecClZ[irow] = c->GetZ();
1314           //
1315 //        Float_t gxyz[3];
1316 //        c->GetGlobalXYZ(gxyz);
1317 //        vecgX[irow]   = gxyz[0];
1318 //        vecgY[irow]   = gxyz[1];
1319 //        vecgZ[irow]   = gxyz[2];
1320           //
1321           Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
1322           Double_t y = vecClY[irow];
1323           Double_t z = vecClZ[irow];
1324           //
1325           Double_t x2[2]={x,x*x};   //linear and parabolic parameters
1326           Double_t x4[4]={0,0,0,0}; //hyp4 parameters
1327           Double_t xIO[2]={0,x};    //common linear + offset IROC-OROC
1328           if ( roc == innerSector ) {
1329               x4[0]=1; //offset inner - outer sector
1330               x4[1]=x; //slope parameter inner sector
1331               xIO[0]=1;
1332           } else {
1333               x4[2]=x; //slope parameter outer sector
1334           }
1335           x4[3]=x*x;   //common parabolic shape
1336           if (pedgeX < fEdgeXcuts[icut]) continue;
1337           if (pedgeY < fEdgeYcuts[icut]) continue;
1338           if (c->GetMax()>900) continue;  // cluster in saturation
1339           //
1340           if ( roc==innerSector ){
1341               fy1I.AddPoint(x2,y);
1342               fz1I.AddPoint(x2,z);
1343               fy2I.AddPoint(x2,y);
1344               fz2I.AddPoint(x2,z);
1345               ++nclI;
1346               if (c->GetX()<xinMin) xinMin=c->GetX();
1347               if (c->GetX()>xinMax) xinMax=c->GetX();
1348
1349               msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
1350               mqratioIn +=c->GetMax()/c->GetQ();
1351
1352           }
1353           if ( roc==outerSector ){
1354               fy1O.AddPoint(x2,y);
1355               fz1O.AddPoint(x2,z);
1356               fy2O.AddPoint(x2,y);
1357               fz2O.AddPoint(x2,z);
1358               ++nclO;
1359               if (c->GetX()<xoutMin) xoutMin=c->GetX();
1360               if (c->GetX()>xoutMax) xoutMax=c->GetX();
1361               msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
1362               mqratioOut +=c->GetMax()/c->GetQ();
1363
1364           }
1365           fy4.AddPoint(x4,y);
1366           fz4.AddPoint(x4,z);
1367           fy1IO.AddPoint(xIO,y);
1368           fz1IO.AddPoint(xIO,z);
1369       }
1370       if (nclI>0)  {
1371         msigmaYIn/=nclI;
1372         mqratioIn/=nclI;
1373       }
1374       if (nclO>0)  {
1375         msigmaYOut/=nclO;
1376         mqratioOut/=nclO;
1377       }
1378       //======================================//
1379       // Evaluate and retrieve fit parameters //
1380       //======================================//
1381       AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
1382       //inner sector
1383       if (  (innerSector>-1) && (fy1I.GetNpoints()>0) ){
1384           fy1I.Eval();
1385           fz1I.Eval();
1386           fy2I.Eval();
1387           fz2I.Eval();
1388           fy1I.GetParameters(vecy1resInner);
1389           fz1I.GetParameters(vecz1resInner);
1390           fy2I.GetParameters(vecy2resInner);
1391           fz2I.GetParameters(vecz2resInner);
1392           chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
1393           chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
1394           chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
1395           chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
1396       }
1397       //outer sector
1398       if (  (outerSector>-1) && (fy1O.GetNpoints()>0) ){
1399           fy1O.Eval();
1400           fz1O.Eval();
1401           fy2O.Eval();
1402           fz2O.Eval();
1403           fy1O.GetParameters(vecy1resOuter);
1404           fz1O.GetParameters(vecz1resOuter);
1405           fy2O.GetParameters(vecy2resOuter);
1406           fz2O.GetParameters(vecz2resOuter);
1407           chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
1408           chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
1409           chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
1410           chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
1411       }
1412       //combined hyp4 fit
1413       if ( innerSector>0 && outerSector>0 ){
1414           if (fy4.GetNpoints()>0) {
1415               fy4.Eval();
1416               fy4.GetParameters(vecy4res);
1417               chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
1418           }
1419           if (fz4.GetNpoints()>0) {
1420               fz4.Eval();
1421               fz4.GetParameters(vecz4res);
1422               chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
1423           }
1424           if (fy1IO.GetNpoints()>0) {
1425             fy1IO.Eval();
1426             fy1IO.GetParameters(vecy1resIO);
1427             chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
1428           }
1429           if (fz1IO.GetNpoints()>0) {
1430             fz1IO.Eval();
1431             fz1IO.GetParameters(vecz1resIO);
1432             chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
1433           }
1434       }
1435       //clear points
1436       fy4.ClearPoints();  fz4.ClearPoints();
1437       fy1I.ClearPoints(); fy1O.ClearPoints();
1438       fz1I.ClearPoints(); fz1O.ClearPoints();
1439       fy2I.ClearPoints(); fy2O.ClearPoints();
1440       fz2I.ClearPoints(); fz2O.ClearPoints();
1441       fy1IO.ClearPoints(); fz1IO.ClearPoints();
1442       //==============================//
1443       // calculate tracklet positions //
1444       //==============================//
1445       AliDebug(4,"Calculate tracklet positions");
1446       for (Int_t irow=158;irow>-1;--irow) {
1447         isReject[irow]=0;
1448          AliTPCclusterMI *c=track->GetClusterPointer(irow);
1449          if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
1450            isReject[irow]+=1;
1451          }
1452
1453          if (!c) { //no cluster
1454            isReject[irow]+=2;
1455          }else{
1456            if (c->GetMax()>kMax){   //saturation
1457              isReject[irow]+=4;
1458            }
1459            if ( vecSec[irow] == outerSector ) {   //extended shape
1460              if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut)                isReject[irow]+=8;
1461              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut)    isReject[irow]+=16;
1462            }else{
1463              if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut)                isReject[irow]+=8;
1464              if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut)    isReject[irow]+=16;
1465            }
1466          }
1467
1468
1469
1470           if ( vecSec[irow]==-1 ) continue;  //no cluster info
1471           if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
1472           tr=&dummy;
1473           Double_t x=vecX[irow];
1474           Double_t xref=x-133.4;
1475           //
1476           Double_t yoffInner=0;
1477           Double_t zoffInner=0;
1478           Double_t yoffInner1=0;
1479           Double_t zoffInner1=0;
1480           Double_t yslopeInner=0;
1481           Double_t yslopeOuter=0;
1482           Double_t zslopeInner=0;
1483           Double_t zslopeOuter=0;
1484           //positions of hyperplane fits
1485           if ( vecSec[irow] == outerSector ) {
1486               tr=trOuter;
1487               vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
1488               vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
1489               vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
1490               vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
1491               yslopeOuter=vecy4res[3];
1492               zslopeOuter=vecz4res[3];
1493           } else {
1494               tr=trInner;
1495               vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
1496               vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
1497               vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
1498               vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
1499               yoffInner=vecy4res[1];
1500               zoffInner=vecz4res[1];
1501               yoffInner1=vecy1resIO[1];
1502               zoffInner1=vecz1resIO[1];
1503               yslopeInner=vecy4res[2];
1504               zslopeInner=vecz4res[2];
1505           }
1506           vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
1507           vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
1508           vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
1509           vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
1510           //positions of kalman fits
1511           Double_t gxyz[3],xyz[3];
1512           AliExternalTrackParam *param = 0x0;
1513           //
1514           param=tr->GetInner();
1515           if (param){
1516               param->GetXYZ(gxyz);
1517               Float_t bz = AliTracker::GetBz(gxyz);
1518               param->GetYAt(x, bz, xyz[1]);
1519               param->GetZAt(x, bz, xyz[2]);
1520               vecYkalman[irow]=xyz[1];
1521               vecZkalman[irow]=xyz[2];
1522           }
1523           //
1524           //
1525           //
1526
1527       }
1528       //=====================================================================//
1529       // write results from the different tracklet fits with debug streamers //
1530       //=====================================================================//
1531       if (fStreamLevel>4){
1532           TTreeSRedirector *cstream = GetDebugStreamer();
1533           if (cstream){
1534             Float_t dedx = track->GetdEdx();
1535             (*cstream)<<"FitModels"<<
1536               "run="<<fRun<<              //  run number
1537               "event="<<fEvent<<          //  event number
1538               "time="<<fTime<<            //  time stamp of event
1539               "trigger="<<fTrigger<<      //  trigger
1540               "mag="<<fMagF<<             //  magnetic field
1541               //
1542               "cutNr="      << icut <<
1543               "edgeCutX="   << edgeCutX <<
1544               "edgeCutY="   << edgeCutY <<
1545               "nclCut="     << nclCut <<
1546               "innerSector="<< innerSector <<
1547               "outerSector="<< outerSector <<
1548               "dEdx="       << dedx <<
1549               "LTr.="       << ltrp <<
1550               "Tr.="        << extparam <<
1551               "yPol1In.="   << &vecy1resInner <<
1552               "zPol1In.="   << &vecz1resInner <<
1553               "yPol1InOut.="<< &vecy1resIO <<
1554               "zPol1InOut.="<< &vecz1resIO <<
1555               "yPol2In.="   << &vecy2resInner <<
1556               "zPol2In.="   << &vecz2resInner <<
1557               "yPol1Out.="  << &vecy1resOuter <<
1558               "zPol1Out.="  << &vecz1resOuter <<
1559               "yPol2Out.="  << &vecy2resOuter <<
1560               "zPol2Out.="  << &vecz2resOuter <<
1561               "yInOut.="    << &vecy4res <<
1562               "zInOut.="    << &vecz4res <<
1563               "chi2y1In="   << chi2I1y <<
1564               "chi2z1In="   << chi2I1z <<
1565               "chi2y1InOut="<< chi2IO1y <<
1566               "chi2z1InOut="<< chi2IO1z <<
1567               "chi2y1Out="  << chi2O1y <<
1568               "chi2z1Out="  << chi2O1z <<
1569               "chi2y2In="   << chi2I2y <<
1570               "chi2z2In="   << chi2I2z <<
1571               "chi2y2Out="  << chi2O2y <<
1572               "chi2z2Out="  << chi2O2z <<
1573               "chi2yInOut=" << chi2IOy <<
1574               "chi2zInOut=" << chi2IOz <<
1575               "trletIn.="   << trInner <<
1576               "trletOut.="  << trOuter <<
1577               "nclI="       << nclI <<
1578               "nclO="       << nclO <<
1579               "xinMin="     << xinMin<<
1580               "xinMax="     << xinMax<<
1581               "xoutMin="    << xoutMin<<
1582               "xoutMax="     << xoutMax<<
1583               "msigmaYIn=" <<msigmaYIn<<
1584               "msigmaYOut=" <<msigmaYOut<<
1585               "mqratioIn=" <<mqratioIn<<
1586               "mqratioOut=" <<  mqratioOut      <<
1587               "\n";
1588           }
1589       }
1590
1591       // wirte residuals information
1592       if (fStreamLevel>5){
1593           TTreeSRedirector *cstream = GetDebugStreamer();
1594           if (cstream){
1595             Float_t dedx = track->GetdEdx();
1596             (*cstream)<<"Residuals"<<
1597               "run="<<fRun<<              //  run number
1598               "event="<<fEvent<<          //  event number
1599               "time="<<fTime<<            //  time stamp of event
1600               "trigger="<<fTrigger<<      //  trigger
1601               "mag="<<fMagF<<             //  magnetic field
1602               //
1603               "cutNr="      << icut <<
1604               "edgeCutX="   << edgeCutX <<
1605               "edgeCutY="   << edgeCutY <<
1606               "nclCut="     << nclCut   <<
1607               "LTr.="       << ltrp <<
1608               "Tr.="        << extparam<<
1609               "dEdx="       << dedx <<
1610               "Cl.="        << &arrCl <<
1611               "vX.="        << &vecgX<<   // global x
1612               "vY.="        << &vecgY<<   // global y
1613               "vZ.="        << &vecgZ<<   // global z
1614               "TrX.="       << &vecX <<
1615               "TrYpol1.="   << &vecY1 <<
1616               "TrZpol1.="   << &vecZ1 <<
1617               "TrYpol2.="   << &vecY2 <<
1618               "TrZpol2.="   << &vecZ2 <<
1619               "TrYpol1InOut.="<< &vecY1IO <<
1620               "TrZpol1InOut.="<< &vecZ1IO <<
1621               "TrYInOut.="  << &vecY4 <<
1622               "TrZInOut.="  << &vecZ4 <<
1623               "ClY.="       << &vecClY <<
1624               "ClZ.="       << &vecClZ <<
1625               "isReject.="  << &isReject<<
1626               "sec.="       << &vecSec <<
1627               "nclI="       << nclI <<
1628               "nclO="       << nclO <<
1629               "xinMin="     << xinMin<<
1630               "xinMax="     << xinMax<<
1631               "xoutMin="    << xoutMin<<
1632               "xoutMax="     << xoutMax<<
1633               "msigmaYIn=" <<msigmaYIn<<
1634               "msigmaYOut=" <<msigmaYOut<<
1635               "mqratioIn=" <<mqratioIn<<
1636               "mqratioOut=" <<  mqratioOut      <<
1637               "yInOut.="    << &vecy4res <<
1638               "zInOut.="    << &vecz4res <<
1639               //chi2s
1640               "chi2y1In="   << chi2I1y <<    //
1641               "chi2z1In="   << chi2I1z <<
1642               "chi2y1Out="  << chi2O1y <<
1643               "chi2z1Out="  << chi2O1z <<
1644               "chi2y1InOut="<< chi2IO1y <<
1645               "chi2z1InOut="<< chi2IO1z <<
1646               "chi2y2In="   << chi2I2y <<
1647               "chi2z2In="   << chi2I2z <<
1648               "chi2y2Out="  << chi2O2y <<
1649               "chi2z2Out="  << chi2O2z <<
1650               "chi2yInOut=" << chi2IOy <<
1651               "chi2zInOut=" << chi2IOz <<
1652               // fit parameters
1653               "yPol1In.="   << &vecy1resInner <<
1654               "zPol1In.="   << &vecz1resInner <<
1655               "yPol2In.="   << &vecy2resInner <<
1656               "zPol2In.="   << &vecz2resInner <<
1657               "yPol1Out.="  << &vecy1resOuter <<
1658               "zPol1Out.="  << &vecz1resOuter <<
1659               "yPol1InOut.="<< &vecy1resIO <<
1660               "zPol1InOut.="<< &vecz1resIO <<
1661               "yPol2Out.="  << &vecy2resOuter <<
1662               "zPol2Out.="  << &vecz2resOuter <<
1663
1664               "\n";
1665
1666           }
1667       }
1668       //==========================//
1669       // Fill Residual Histograms //
1670       //==========================//
1671       if (!fHisNclIn) MakeFitHistos();
1672
1673       TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
1674       TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
1675       TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
1676       TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
1677       TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
1678       TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
1679       //      TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
1680       //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
1681       //
1682       for (Int_t irow=158;irow>-1;--irow) {
1683         if (vecSec[irow]==-1)continue; // no cluster info
1684         if (isReject[irow]>0.5) continue;  //
1685         Double_t ycl  = vecClY[irow];
1686         Double_t yfit = vecY1[irow];
1687         Double_t yfit2 = vecY2[irow];
1688         Double_t x    = vecX[irow];
1689         Double_t yabsbeam = -1000;
1690         if(vecSec[irow]==outerSector && (outerSector%36)==fBeamSectorOuter[id])
1691           yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
1692         else if((innerSector%36)==fBeamSectorInner[id])
1693           yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
1694
1695         //      Double_t yfit3 = vecY2[irow];
1696         Double_t zcl  = vecClZ[irow];
1697         Double_t zfit = vecZ1[irow];
1698         Double_t zfit2 = vecZ2[irow];
1699         //Double_t zfit3 = vecZ2[irow];
1700
1701         // dz abs
1702         // The expressions for zcorrected has been obtained by
1703         // inverting the fits in the FitDriftV() method (ignoring the
1704         // global y dependence for now):
1705         // A side:
1706         // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
1707         // =>
1708         // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
1709         //
1710         // C side:
1711         // 250 + zmeasured  = [0] + [1]*(250+zreal) + .... (yglobal)
1712         // =>
1713         // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
1714
1715         Double_t dzabs = -1000;
1716         if (ltrp->GetSide()==0){
1717           if ((*fFitAside)[1]>0.) {
1718             // ignore global y dependence for now
1719             Double_t zcorrected = (zcl + (*fFitAside)[0] -
1720                                    (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
1721             //      zcorrected = zcl;
1722             if(vecSec[irow]==outerSector && (outerSector%36)==fBeamSectorOuter[id])
1723               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
1724             else if((innerSector%36)==fBeamSectorInner[id])
1725               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
1726           }
1727         } else {
1728           if ((*fFitCside)[1]>0.) {
1729             Double_t zcorrected = (zcl - (*fFitCside)[0] +
1730                                    (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
1731             //      zcorrected = zcl;
1732             if(vecSec[irow]==outerSector && (outerSector%36)==fBeamSectorOuter[id])
1733               dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
1734             else if((innerSector%36)==fBeamSectorInner[id])
1735               dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
1736           }
1737         }
1738
1739         if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
1740           if (profy){
1741               profy->Fill(irow,ycl-yfit);
1742               profy2->Fill(irow,ycl-yfit2);
1743               if(yabsbeam<-100)
1744                 profyabs->Fill(irow,-0.99);
1745               else
1746                 profyabs->Fill(irow,ycl-yabsbeam);
1747
1748               //              profy3->Fill(irow,ycl-yfit3);
1749           }
1750           if (profz) {
1751               profz->Fill(irow,zcl-zfit);
1752               profz2->Fill(irow,zcl-zfit2);
1753               //profz3->Fill(irow,zcl-zfit3);
1754               if(dzabs<-100)
1755                 profzabs->Fill(irow,-0.99);
1756               else
1757                 profzabs->Fill(irow,dzabs);
1758           }
1759         }
1760       }
1761       //
1762       //
1763       // Fill laser fit histograms
1764       //
1765       Float_t dedx = track->GetdEdx();
1766       if (nclI>20&&nclO>20){
1767         fHisNclIn->Fill(id,nclI);      //->Number of clusters inner
1768         fHisNclOut->Fill(id,nclO);     //->Number of clusters outer
1769         fHisNclIO->Fill(id,nclI+nclO);      //->Number of cluster inner outer
1770         //
1771         fHisLclIn->Fill(id,xinMax-xinMin);        //->Level arm inner
1772         fHisLclOut->Fill(id,xoutMax-xoutMin);     //->Level arm outer
1773         fHisLclIO->Fill(id,xoutMax-xinMin);       //->Number of cluster inner outer
1774         //
1775         fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
1776         fHisdZfit->Fill(id,fFitZ[id]);
1777         //
1778         //
1779         fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y));      //->chi2 y inner - line
1780         fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y));     //->chi2 y inner - line
1781         fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y));      //->chi2 y inner - parabola
1782         fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y));     //->chi2 y inner - parabola
1783         fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy));      //->chi2 y IO    - common
1784
1785
1786         fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z));      //->chi2 z inner - line
1787         fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z));     //->chi2 z inner - line
1788         fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y));      //->chi2 z inner - parabola
1789         fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z));     //->chi2 z inner - parabola
1790         fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz));      //->chi2 z IO    - common
1791         //
1792         //
1793         fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]);     //-> delta y   P0outer-P0inner - line
1794         fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]);     //-> delta y   P0outer-P0inner - parabola
1795         fHisPy3vP0->Fill(id,vecy4res[1]);                           //-> delta y   P0outer-P0inner - common parabola
1796         //
1797         fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]);     //-> delta ky  P1outer-P1inner - line
1798         fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]);     //-> delta ky  P1outer-P1inner - parabola
1799         fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]);               //-> delta ky  P1outer-P1inner - common parabola
1800         //
1801         fHisPy3vP2IO->Fill(id,vecy4res[4]);        //-> Curv  P2outerinner - common parabola
1802         fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]);     //-> delta z   P0outer-P0inner - line
1803         fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]);     //-> delta z   P0outer-P0inner - parabola
1804         fHisPz3vP0->Fill(id,vecz4res[1]);                           //-> delta z   P0outer-P0inner - common parabola
1805         //
1806         fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]);     //-> delta kz  P1outer-P1inner - line
1807         fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]);     //-> delta kz  P1outer-P1inner - parabola
1808         fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]);               //-> delta kz  P1outer-P1inner - common parabola
1809         fHisPz3vP2IO->Fill(id,vecz4res[4]);        //-> Curv  P2outerinner - common parabola
1810       }
1811       if(nclI>20){
1812         fHisPy2vP2In->Fill(id,vecy2resInner[2]);   //-> Curv  P2inner - parabola
1813         fHisPz2vP2In->Fill(id,vecz2resInner[2]);   //-> Curv  P2inner - parabola
1814       }
1815       //
1816       if (nclO>20){
1817         fHisPz2vP2Out->Fill(id,vecz2resOuter[2]);  //-> Curv  P2outer - parabola
1818         fHisPy2vP2Out->Fill(id,vecy2resOuter[2]);  //-> Curv  P2outer - parabola
1819       }
1820
1821   }
1822 }
1823
1824
1825
1826 void AliTPCcalibLaser::DumpMeanInfo(Float_t bfield, Int_t run){
1827   //
1828   //  Dump information about laser beams
1829   //  isOK variable indicates usability of the beam
1830   //  Beam is not usable if:
1831   //     a.  No entries in range (krmsCut0)
1832   //     b.  Big sperad          (krmscut1)
1833   //     c.  RMSto fit sigma bigger then (kmultiCut)
1834   //     d.  Too big angular spread
1835   //
1836
1837   const Float_t krmsCut0=0.001;
1838   const Float_t krmsCut1=0.16;
1839   const Float_t kmultiCut=2;
1840   const Float_t kcutP0=0.002;
1841   //
1842   AliTPCcalibLaser *laser = this;
1843   TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
1844   TF1 fg("fg","gaus");
1845   //
1846   //
1847   for (Int_t id=0; id<336; id++){
1848     Bool_t isOK=kTRUE;
1849     TH1F * hisphi  = (TH1F*)laser->fDeltaPhi.At(id);
1850     TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
1851     TH1F * hisZ    = (TH1F*)laser->fDeltaZ.At(id);
1852     TH1F * hisP3    = (TH1F*)laser->fDeltaP3.At(id);
1853     TH1F * hisP4    = (TH1F*)laser->fDeltaP4.At(id);
1854     TH1F * hisS    = (TH1F*)laser->fSignals.At(id);
1855     //if (!hisphi) continue;
1856     Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
1857     //if (entries<minEntries) continue;
1858     //
1859     AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
1860     if (!ltrp) {
1861      AliTPCLaserTrack::LoadTracks();
1862       ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
1863     }
1864     pcstream->GetFile()->cd();
1865     if (hisphi)  hisphi->Write();
1866     if (hisphiP) hisphiP->Write();
1867     if (hisZ)    hisZ->Write();
1868     if (hisP3)    hisP3->Write();
1869     if (hisP4)    hisP4->Write();
1870
1871     Float_t meanphi = hisphi->GetMean();
1872     Float_t rmsphi = hisphi->GetRMS();
1873     //
1874     Float_t meanphiP = hisphiP->GetMean();
1875     Float_t rmsphiP = hisphiP->GetRMS();
1876     Float_t meanZ = hisZ->GetMean();
1877     Float_t rmsZ = hisZ->GetRMS();
1878     if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
1879       hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
1880     Double_t gphi1 = fg.GetParameter(1);
1881     Double_t gphi2 = fg.GetParameter(2);
1882     if (hisphiP->GetRMS()>0)
1883       hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
1884     Double_t gphiP1 = fg.GetParameter(1);
1885     Double_t gphiP2 = fg.GetParameter(2);
1886     //
1887     if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
1888       hisZ->Fit(&fg,"","");
1889     Double_t gz1 = fg.GetParameter(1);
1890     Double_t gz2 = fg.GetParameter(2);
1891     //
1892     if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
1893       hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
1894     Double_t gp31 = fg.GetParameter(1);
1895     Double_t gp32 = fg.GetParameter(2);
1896     //
1897     if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
1898       hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
1899     Double_t gp41 = fg.GetParameter(1);
1900     Double_t gp42 = fg.GetParameter(2);
1901     //
1902     Float_t meanS=hisS->GetMean();
1903     //
1904     Double_t lxyz[3];
1905     Double_t lpxyz[3];
1906     ltrp->GetXYZ(lxyz);
1907     ltrp->GetPxPyPz(lpxyz);
1908
1909     if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
1910     if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
1911     if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE;   // multi peak structure
1912     if (gphiP2>kcutP0) isOK=kFALSE;
1913     //
1914     //
1915     //
1916     TH1 *his =0;
1917     //
1918     his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
1919     Float_t mnclIn =  his->GetMean();
1920     delete his;
1921     his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
1922     Float_t mnclOut =  his->GetMean();
1923     delete his;
1924     his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
1925     Float_t mnclIO =  his->GetMean();
1926     delete his;
1927     his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
1928     Float_t mLclIn =  his->GetMean();
1929     delete his;
1930     his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
1931     Float_t mLclOut =  his->GetMean();
1932     delete his;
1933     his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
1934     Float_t mLclIO =  his->GetMean();
1935     delete his;
1936     //
1937     his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
1938     Float_t mdEdx =  his->GetMean();
1939     delete his;
1940     //
1941     //
1942     //
1943     //
1944     his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - line
1945     Float_t  mChi2YIn1=  his->GetMean();
1946     delete his;
1947     his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - line
1948     Float_t mChi2YOut1 =  his->GetMean();
1949     delete his;
1950     his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 y inner - parabola
1951     Float_t mChi2YIn2 =  his->GetMean();
1952     delete his;
1953     his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 y inner - parabola
1954     Float_t mChi2YOut2 =  his->GetMean();
1955     delete his;
1956     his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 y IO    - common
1957     Float_t mChi2YIO1 =  his->GetMean();
1958     delete his;
1959     his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - line
1960     Float_t mChi2ZIn1 =  his->GetMean();
1961     delete his;
1962     his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - line
1963     Float_t mChi2ZOut1 =  his->GetMean();
1964     delete his;
1965     his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1);      //->chi2 z inner - parabola
1966     Float_t mChi2ZIn2 =  his->GetMean();
1967     delete his;
1968     his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1);     //->chi2 z inner - parabola
1969     Float_t mChi2ZOut2 =  his->GetMean();
1970     delete his;
1971     his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1);      //->chi2 z IO    - common
1972     Float_t mChi2ZIO1 =  his->GetMean();
1973     delete his;
1974     //
1975     // fit res. histos
1976     //
1977     his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
1978     Float_t  edZfit  = his->GetEntries();
1979     Float_t  mdZfit =  his->GetMean();
1980     Float_t  rdZfit  = his->GetRMS();
1981     delete his;
1982
1983     his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - line
1984     Float_t  ePy1vP0  = his->GetEntries();
1985     Float_t  mPy1vP0  = his->GetMean();
1986     Float_t  rPy1vP0  = his->GetRMS();
1987     delete his;
1988
1989     his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - parabola
1990     Float_t  ePy2vP0  = his->GetEntries();
1991     Float_t  mPy2vP0  = his->GetMean();
1992     Float_t  rPy2vP0  = his->GetRMS();
1993     delete his;
1994
1995     his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta y   P0outer-P0inner - common parabola
1996     Float_t  ePy3vP0  = his->GetEntries();
1997     Float_t  mPy3vP0  = his->GetMean();
1998     Float_t  rPy3vP0  = his->GetRMS();
1999     delete his;
2000
2001     his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - line
2002     Float_t  ePy1vP1  = his->GetEntries();
2003     Float_t  mPy1vP1  = his->GetMean();
2004     Float_t  rPy1vP1  = his->GetRMS();
2005     delete his;
2006
2007     his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - parabola
2008     Float_t  ePy2vP1  = his->GetEntries();
2009     Float_t  mPy2vP1  = his->GetMean();
2010     Float_t  rPy2vP1  = his->GetRMS();
2011     delete his;
2012
2013     his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta ky  P1outer-P1inner - common parabola
2014     Float_t  ePy3vP1  = his->GetEntries();
2015     Float_t  mPy3vP1  = his->GetMean();
2016     Float_t  rPy3vP1  = his->GetRMS();
2017     delete his;
2018
2019     his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2020     Float_t  ePy2vP2In  = his->GetEntries();
2021     Float_t  mPy2vP2In  = his->GetMean();
2022     Float_t  rPy2vP2In  = his->GetRMS();
2023     delete his;
2024
2025     his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2026     Float_t  ePy2vP2Out  = his->GetEntries();
2027     Float_t  mPy2vP2Out  = his->GetMean();
2028     Float_t  rPy2vP2Out  = his->GetRMS();
2029     delete his;
2030
2031     his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2032     Float_t  ePy3vP2IO  = his->GetEntries();
2033     Float_t  mPy3vP2IO  = his->GetMean();
2034     Float_t  rPy3vP2IO  = his->GetRMS();
2035     delete his;
2036
2037     //
2038     //
2039     his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - line
2040     Float_t  ePz1vP0  = his->GetEntries();
2041     Float_t  mPz1vP0  = his->GetMean();
2042     Float_t  rPz1vP0  = his->GetRMS();
2043     delete his;
2044
2045     his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - parabola
2046     Float_t  ePz2vP0  = his->GetEntries();
2047     Float_t  mPz2vP0  = his->GetMean();
2048     Float_t  rPz2vP0  = his->GetRMS();
2049     delete his;
2050
2051     his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1);     //-> delta z   P0outer-P0inner - common parabola
2052     Float_t  ePz3vP0  = his->GetEntries();
2053     Float_t  mPz3vP0  = his->GetMean();
2054     Float_t  rPz3vP0  = his->GetRMS();
2055     delete his;
2056
2057     his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - line
2058     Float_t  ePz1vP1  = his->GetEntries();
2059     Float_t  mPz1vP1  = his->GetMean();
2060     Float_t  rPz1vP1  = his->GetRMS();
2061     delete his;
2062
2063     his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - parabola
2064     Float_t  ePz2vP1  = his->GetEntries();
2065     Float_t  mPz2vP1  = his->GetMean();
2066     Float_t  rPz2vP1  = his->GetRMS();
2067     delete his;
2068
2069     his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1);     //-> delta kz  P1outer-P1inner - common parabola
2070     Float_t  ePz3vP1  = his->GetEntries();
2071     Float_t  mPz3vP1  = his->GetMean();
2072     Float_t  rPz3vP1  = his->GetRMS();
2073     delete his;
2074
2075     his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2inner - parabola
2076     Float_t  ePz2vP2In  = his->GetEntries();
2077     Float_t  mPz2vP2In  = his->GetMean();
2078     Float_t  rPz2vP2In  = his->GetRMS();
2079     delete his;
2080
2081     his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1);  //-> Curv  P2outer - parabola
2082     Float_t  ePz2vP2Out  = his->GetEntries();
2083     Float_t  mPz2vP2Out  = his->GetMean();
2084     Float_t  rPz2vP2Out  = his->GetRMS();
2085     delete his;
2086
2087     his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1);   //-> Curv  P2outerinner - common parabola
2088     Float_t  ePz3vP2IO  = his->GetEntries();
2089     Float_t  mPz3vP2IO  = his->GetMean();
2090     Float_t  rPz3vP2IO  = his->GetRMS();
2091     delete his;
2092
2093     //
2094     (*pcstream)<<"Mean"<<
2095       "run="<<run<<               //
2096       "isOK="<<isOK<<             //
2097       "id="<<id<<                 // track id
2098       "entries="<<entries<<       // number of entries
2099       "bz="<<bfield<<             // bfield
2100       "LTr.="<<ltrp<<             // refernece track
2101       //
2102       "mnclIn="<<mnclIn<<         // mean number of clusters in inner
2103       "mnclOut="<<mnclOut<<       // mean number of clusters in outer
2104       "mnclIO="<<mnclIO<<         // mean number of clusters in inner+outer
2105       "mLclIn="<<mLclIn<<         // mean number of clusters in inner
2106       "mLclOut="<<mLclOut<<       // mean number of clusters in outer
2107       "mLclIO="<<mLclIO<<         // mean number of clusters in inner+outer
2108       "mdEdx="<<mdEdx<<           // mean dEdx
2109       "edZfit="<<edZfit<<           // entries z fit
2110       "mdZfit="<<mdZfit<<           // mean z fit
2111       "rdZfit="<<rdZfit<<           // RMS z fit
2112       //
2113       //
2114       "mChi2YIn1="<<mChi2YIn1<<       //->chi2 y inner - line
2115       "mChi2YOut1="<<mChi2YOut1<<     //->chi2 y inner - line
2116       "mChi2YIn2="<<mChi2YIn2<<       //->chi2 y inner - parabola
2117       "mChi2YOut2="<<mChi2YOut2<<     //->chi2 y inner - parabola
2118       "mChi2YIO1="<<mChi2YIO1<<       //->chi2 y IO    - common
2119       "mChi2ZIn1="<<mChi2ZIn1<<       //->chi2 z inner - line
2120       "mChi2ZOut1="<<mChi2ZOut1<<     //->chi2 z inner - line
2121       "mChi2ZIn2="<<mChi2ZIn2<<       //->chi2 z inner - parabola
2122       "mChi2ZOut2="<<mChi2ZOut2<<     //->chi2 z inner - parabola
2123       "mChi2ZIO1="<<mChi2ZIO1<<       //->chi2 z IO    - common
2124       //
2125       //
2126       "ePy1vP0="<<ePy1vP0<<
2127       "mPy1vP0="<<mPy1vP0<<
2128       "rPy1vP0="<<rPy1vP0<<
2129       "ePy2vP0="<<ePy2vP0<<
2130       "mPy2vP0="<<mPy2vP0<<
2131       "rPy2vP0="<<rPy2vP0<<
2132       "ePy3vP0="<<ePy3vP0<<
2133       "mPy3vP0="<<mPy3vP0<<
2134       "rPy3vP0="<<rPy3vP0<<
2135       "ePy1vP1="<<ePy1vP1<<
2136       "mPy1vP1="<<mPy1vP1<<
2137       "rPy1vP1="<<rPy1vP1<<
2138       "ePy2vP1="<<ePy2vP1<<
2139       "mPy2vP1="<<mPy2vP1<<
2140       "rPy2vP1="<<rPy2vP1<<
2141       "ePy3vP1="<<ePy3vP1<<
2142       "mPy3vP1="<<mPy3vP1<<
2143       "rPy3vP1="<<rPy3vP1<<
2144       "ePy2vP2In="<<ePy2vP2In<<
2145       "mPy2vP2In="<<mPy2vP2In<<
2146       "rPy2vP2In="<<rPy2vP2In<<
2147       "ePy2vP2Out="<<ePy2vP2Out<<
2148       "mPy2vP2Out="<<mPy2vP2Out<<
2149       "rPy2vP2Out="<<rPy2vP2Out<<
2150       "ePy3vP2IO="<<ePy3vP2IO<<
2151       "mPy3vP2IO="<<mPy3vP2IO<<
2152       "rPy3vP2IO="<<rPy3vP2IO<<
2153       //
2154       //
2155       "ePz1vP0="<<ePz1vP0<<
2156       "mPz1vP0="<<mPz1vP0<<
2157       "rPz1vP0="<<rPz1vP0<<
2158       "ePz2vP0="<<ePz2vP0<<
2159       "mPz2vP0="<<mPz2vP0<<
2160       "rPz2vP0="<<rPz2vP0<<
2161       "ePz3vP0="<<ePz3vP0<<
2162       "mPz3vP0="<<mPz3vP0<<
2163       "rPz3vP0="<<rPz3vP0<<
2164       "ePz1vP1="<<ePz1vP1<<
2165       "mPz1vP1="<<mPz1vP1<<
2166       "rPz1vP1="<<rPz1vP1<<
2167       "ePz2vP1="<<ePz2vP1<<
2168       "mPz2vP1="<<mPz2vP1<<
2169       "rPz2vP1="<<rPz2vP1<<
2170       "ePz3vP1="<<ePz3vP1<<
2171       "mPz3vP1="<<mPz3vP1<<
2172       "rPz3vP1="<<rPz3vP1<<
2173       "ePz2vP2In="<<ePz2vP2In<<
2174       "mPz2vP2In="<<mPz2vP2In<<
2175       "rPz2vP2In="<<rPz2vP2In<<
2176       "ePz2vP2Out="<<ePz2vP2Out<<
2177       "mPz2vP2Out="<<mPz2vP2Out<<
2178       "rPz2vP2Out="<<rPz2vP2Out<<
2179       "ePz3vP2IO="<<ePz3vP2IO<<
2180       "mPz3vP2IO="<<mPz3vP2IO<<
2181       "rPz3vP2IO="<<rPz3vP2IO<<
2182       //
2183       //
2184       //
2185       "lx0="<<lxyz[0]<<            // reference x
2186       "lx1="<<lxyz[1]<<            // reference y
2187       "lx2="<<lxyz[2]<<            // refernece z
2188       "lpx0="<<lpxyz[0]<<           // reference x
2189       "lpx1="<<lpxyz[1]<<          // reference y
2190       "lpx2="<<lpxyz[2]<<          // refernece z
2191       //
2192       "msig="<<meanS<<
2193       //
2194       "mphi="<<meanphi<<         //
2195       "rmsphi="<<rmsphi<<        //
2196       "gphi1="<<gphi1<<
2197       "gphi2="<<gphi2<<
2198       //
2199       "mphiP="<<meanphiP<<       //
2200       "rmsphiP="<<rmsphiP<<      //
2201       "gphiP1="<<gphiP1<<
2202       "gphiP2="<<gphiP2<<
2203       //
2204       "meanZ="<<meanZ<<
2205       "rmsZ="<<rmsZ<<
2206       "gz1="<<gz1<<
2207       "gz2="<<gz2<<
2208       //
2209       "gp31="<<gp31<<            //gaus mean - tgl
2210       "gp32="<<gp32<<            //gaus rms  -tgl
2211       "gp41="<<gp41<<            //gaus mean - P4
2212       "gp42="<<gp42<<            //gaus rms  - P4
2213
2214       "\n";
2215   }
2216   delete pcstream;
2217   /*
2218     Browse the content
2219     TFile fmean("laserMean.root");
2220
2221
2222    */
2223
2224
2225 }
2226
2227
2228
2229 void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
2230   //
2231   //
2232   //
2233   TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
2234   TFile * f = pcstream->GetFile();
2235   f->mkdir("dirphi");
2236   f->mkdir("dirphiP");
2237   f->mkdir("dirZ");
2238   TF1 fp("p1","pol1");
2239   //
2240   //
2241   char    cut[1000];
2242   char grname[1000];
2243   char grnamefull[1000];
2244
2245   Double_t mphi[100];
2246   Double_t mphiP[100];
2247   Double_t smphi[100];
2248   Double_t smphiP[100];
2249   Double_t mZ[100];
2250   Double_t smZ[100];
2251   Double_t bz[100];
2252   Double_t sbz[100];
2253   // fit parameters
2254   Double_t pphi[3];
2255   Double_t pphiP[3];
2256   Double_t pmZ[3];
2257
2258   //
2259   for (Int_t id=0; id<336; id++){
2260     // id =205;
2261     sprintf(cut,"fId==%d&&%s",id,cutUser);
2262     Int_t entries = chain->Draw("bz",cut,"goff");
2263     if (entries<3) continue;
2264     AliTPCLaserTrack *ltrp = 0;
2265     if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
2266     ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
2267     Double_t lxyz[3];
2268     Double_t lpxyz[3];
2269     ltrp->GetXYZ(lxyz);
2270     ltrp->GetPxPyPz(lpxyz);
2271
2272     chain->Draw("bz",cut,"goff");
2273     memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
2274     chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
2275     memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
2276     //
2277     chain->Draw("gphi1",cut,"goff");
2278     memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
2279     chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
2280     memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
2281     //
2282     chain->Draw("gphiP1",cut,"goff");
2283     memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
2284     chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
2285     memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
2286     //
2287     chain->Draw("gz1",cut,"goff");
2288     memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
2289     chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
2290     memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
2291     //
2292     //
2293     sprintf(grnamefull,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
2294             ltrp->GetSide(),  ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
2295     // store data
2296     // phi
2297     f->cd("dirphi");
2298     Float_t phiB0 =0;
2299     Float_t phiPB0=0;
2300     Float_t zB0=0;
2301     for (Int_t ientry=0; ientry<entries; ientry++){
2302       if (TMath::Abs(bz[ientry])<0.05){
2303         phiB0  = mphi[ientry];
2304         phiPB0 = mphiP[ientry];
2305         zB0    = mZ[ientry];
2306       }
2307     }
2308     TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
2309     grphi->Draw("a*");
2310     grphi->Fit(&fp);
2311     pphi[0] = fp.GetParameter(0);                          // offset
2312     pphi[1] = fp.GetParameter(1);                          // slope
2313     pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
2314     sprintf(grname,"phi_id%d",id);
2315     grphi->SetName(grname);  grphi->SetTitle(grnamefull);
2316     grphi->GetXaxis()->SetTitle("b_{z} (T)");
2317     grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
2318     grphi->SetMaximum(1.2);
2319     grphi->SetMinimum(-1.2);
2320     grphi->Draw("a*");
2321
2322     grphi->Write();
2323     gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
2324     // phiP
2325     f->cd("dirphiP)");
2326     TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
2327     grphiP->Draw("a*");
2328     grphiP->Fit(&fp);
2329     pphiP[0] = fp.GetParameter(0);                          // offset
2330     pphiP[1] = fp.GetParameter(1);                          // slope
2331     pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
2332     sprintf(grname,"phiP_id%d",id);
2333     grphiP->SetName(grname);  grphiP->SetTitle(grnamefull);
2334     grphiP->GetXaxis()->SetTitle("b_{z} (T)");
2335     grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
2336     grphiP->SetMaximum(pphiP[0]+0.005);
2337     grphiP->SetMinimum(pphiP[0]-0.005);
2338
2339     gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
2340     grphiP->Write();
2341     //
2342     //Z
2343     f->cd("dirZ");
2344     TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
2345     grmZ->Draw("a*");
2346     grmZ->Fit(&fp);
2347     pmZ[0] = fp.GetParameter(0);                          // offset
2348     pmZ[1] = fp.GetParameter(1);                          // slope
2349     pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.));  // normalized chi2
2350     sprintf(grname,"mZ_id%d",id);
2351     grmZ->SetName(grname);  grmZ->SetTitle(grnamefull);
2352     grmZ->GetXaxis()->SetTitle("b_{z} (T)");
2353     grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
2354
2355     gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
2356     grmZ->Write();
2357     //
2358     // P4
2359     //
2360
2361     for (Int_t ientry=0; ientry<entries; ientry++){
2362       (*pcstream)<<"Mean"<<
2363         "id="<<id<<
2364         "LTr.="<<ltrp<<
2365         "entries="<<entries<<
2366         "bz="<<bz[ientry]<<
2367         "lx0="<<lxyz[0]<<          // reference x
2368         "lx1="<<lxyz[1]<<          // reference y
2369         "lx2="<<lxyz[2]<<          // refernece z
2370         "lpx0="<<lpxyz[0]<<          // reference x
2371         "lpx1="<<lpxyz[1]<<          // reference y
2372         "lpx2="<<lpxyz[2]<<          // refernece z
2373         //values
2374         "phiB0="<<phiB0<<          // position shift at 0 field
2375         "phiPB0="<<phiPB0<<        // angular  shift at 0 field
2376         "zB0="<<zB0<<              // z shift for 0 field
2377         //
2378         "gphi1="<<mphi[ientry]<< // mean - from gaus fit
2379         "pphi0="<<pphi[0]<<   // offset
2380         "pphi1="<<pphi[1]<<   // slope
2381         "pphi2="<<pphi[2]<<   // norm chi2
2382         //
2383         "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
2384         "pphiP0="<<pphiP[0]<< // offset
2385         "pphiP1="<<pphiP[1]<< // slope
2386         "pphiP2="<<pphiP[2]<< // norm chi2
2387         //
2388         "gz1="<<mZ[ientry]<<
2389         "pmZ0="<<pmZ[0]<<     // offset
2390         "pmZ1="<<pmZ[1]<<     // slope
2391         "pmZ2="<<pmZ[2]<<     // norm chi2
2392         "\n";
2393     }
2394   }
2395
2396   delete pcstream;
2397
2398 }
2399
2400
2401 void AliTPCcalibLaser::Analyze(){
2402   //
2403   //
2404   //
2405 }
2406
2407
2408 Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
2409
2410   TIterator* iter = li->MakeIterator();
2411   AliTPCcalibLaser* cal = 0;
2412   static Int_t counter0=0;
2413   while ((cal = (AliTPCcalibLaser*)iter->Next())) {
2414     if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
2415       //      Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
2416       return -1;
2417     }
2418     printf("Marging number %d\n", counter0);
2419     counter0++;
2420     //
2421     MergeFitHistos(cal);
2422     TH1F *h=0x0;
2423     TH1F *hm=0x0;
2424     TH2F *h2=0x0;
2425     TH2F *h2m=0x0;
2426     //    TProfile *hp=0x0;
2427     //TProfile *hpm=0x0;
2428
2429     for (Int_t id=0; id<336; id++){
2430       // merge fDeltaZ histograms
2431       hm = (TH1F*)cal->fDeltaZ.At(id);
2432       h  = (TH1F*)fDeltaZ.At(id);
2433       if (!h) {
2434         h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
2435         h->SetDirectory(0);
2436         fDeltaZ.AddAt(h,id);
2437       }
2438       if (hm) h->Add(hm);
2439       // merge fP3 histograms
2440       hm = (TH1F*)cal->fDeltaP3.At(id);
2441       h  = (TH1F*)fDeltaP3.At(id);
2442       if (!h) {
2443         h=new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
2444         h->SetDirectory(0);
2445         fDeltaP3.AddAt(h,id);
2446       }
2447       if (hm) h->Add(hm);
2448       // merge fP4 histograms
2449       hm = (TH1F*)cal->fDeltaP4.At(id);
2450       h  = (TH1F*)fDeltaP4.At(id);
2451       if (!h) {
2452         h=new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
2453         h->SetDirectory(0);
2454         fDeltaP4.AddAt(h,id);
2455       }
2456       if (hm) h->Add(hm);
2457
2458       //
2459       // merge fDeltaPhi histograms
2460       hm = (TH1F*)cal->fDeltaPhi.At(id);
2461       h  = (TH1F*)fDeltaPhi.At(id);
2462       if (!h) {
2463         h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
2464         h->SetDirectory(0);
2465         fDeltaPhi.AddAt(h,id);
2466       }
2467       if (hm) h->Add(hm);
2468       // merge fDeltaPhiP histograms
2469       hm = (TH1F*)cal->fDeltaPhiP.At(id);
2470       h  = (TH1F*)fDeltaPhiP.At(id);
2471       if (!h) {
2472         h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
2473         h->SetDirectory(0);
2474         fDeltaPhiP.AddAt(h,id);
2475       }
2476       if (hm) h->Add(hm);
2477       // merge fSignals histograms
2478       hm = (TH1F*)cal->fSignals.At(id);
2479       h  = (TH1F*)fSignals.At(id);
2480       if (!h) {
2481         h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
2482         h->SetDirectory(0);
2483         fSignals.AddAt(h,id);
2484       }
2485       if (hm) h->Add(hm);
2486       //
2487       //
2488       // merge ProfileY histograms -0
2489       h2m = (TH2F*)cal->fDeltaYres.At(id);
2490       h2  = (TH2F*)fDeltaYres.At(id);
2491       if (h2m&&h2) h2->Add(h2m);
2492       //
2493       h2m = (TH2F*)cal->fDeltaZres.At(id);
2494       h2  = (TH2F*)fDeltaZres.At(id);
2495       if (h2m&&h2) h2->Add(h2m);
2496       // merge ProfileY histograms - 2
2497       h2m = (TH2F*)cal->fDeltaYres2.At(id);
2498       h2  = (TH2F*)fDeltaYres2.At(id);
2499       if (h2m&&h2) h2->Add(h2m);
2500       //
2501       h2m = (TH2F*)cal->fDeltaZres2.At(id);
2502       h2  = (TH2F*)fDeltaZres2.At(id);
2503       if (h2m&&h2) h2->Add(h2m);
2504
2505       // merge ProfileY histograms - abs
2506       h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
2507       h2  = (TH2F*)fDeltaYresAbs.At(id);
2508       if (h2m&&h2) h2->Add(h2m);
2509
2510       h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
2511       h2  = (TH2F*)fDeltaZresAbs.At(id);
2512       if (h2m&&h2) h2->Add(h2m);
2513       // merge ProfileY histograms - 3
2514       //h2m = (TH2F*)cal->fDeltaYres3.At(id);
2515       //h2  = (TH2F*)fDeltaYres3.At(id);
2516       //if (h2m) h2->Add(h2m);
2517       //
2518       //h2m = (TH2F*)cal->fDeltaZres3.At(id);
2519       //h2  = (TH2F*)fDeltaZres3.At(id);
2520       //if (h2m) h->Add(h2m);
2521       //
2522       //
2523     }
2524   }
2525   return 0;
2526 }
2527
2528 void   AliTPCcalibLaser::MakeFitHistos(){
2529   //
2530   // Make a fit histograms
2531   //
2532   // Number of clusters
2533   //
2534   //TH2F           *fHisNclIn;      //->Number of clusters inner
2535   //TH2F           *fHisNclOut;     //->Number of clusters outer
2536   //TH2F           *fHisNclIO;      //->Number of cluster inner outer
2537   // TH2F           *fHisdEdx;       //->dEdx histo
2538   fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
2539   fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
2540   fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
2541   //
2542   fHisLclIn  = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
2543   fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
2544   fHisLclIO  = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
2545   //
2546   fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
2547   fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
2548
2549   //
2550   //  Chi2
2551   //
2552   //   TH2F           *fHisChi2YIn1;      //->chi2 y inner - line
2553   //   TH2F           *fHisChi2YOut1;     //->chi2 y inner - line
2554   //   TH2F           *fHisChi2YIn2;      //->chi2 y inner - parabola
2555   //   TH2F           *fHisChi2YOut2;     //->chi2 y inner - parabola
2556   //   TH2F           *fHisChi2YIO1;      //->chi2 y IO    - common
2557   fHisChi2YIn1   = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
2558   fHisChi2YOut1  = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
2559   fHisChi2YIn2   = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
2560   fHisChi2YOut2  = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
2561   fHisChi2YIO1   = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
2562   //   TH2F           *fHisChi2ZIn1;      //->chi2 z inner - line
2563   //   TH2F           *fHisChi2ZOut1;     //->chi2 z inner - line
2564   //   TH2F           *fHisChi2ZIn2;      //->chi2 z inner - parabola
2565   //   TH2F           *fHisChi2ZOut2;     //->chi2 z inner - parabola
2566   //   TH2F           *fHisChi2ZIO1;      //->chi2 z IO    - common
2567   fHisChi2ZIn1   = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
2568   fHisChi2ZOut1  = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
2569   fHisChi2ZIn2   = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
2570   fHisChi2ZOut2  = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
2571   fHisChi2ZIO1   = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
2572   //
2573   // Fit
2574   //
2575   //
2576   //   TH2F           *fHisPy1vP0;     //-> delta y   P0outer-P0inner - line
2577   //   TH2F           *fHisPy2vP0;     //-> delta y   P0outer-P0inner - parabola
2578   //   TH2F           *fHisPy3vP0;     //-> delta y   P0outer-P0inner - common parabola
2579   //   TH2F           *fHisPy1vP1;     //-> delta ky  P1outer-P1inner - line
2580   //   TH2F           *fHisPy2vP1;     //-> delta ky  P1outer-P1inner - parabola
2581   //   TH2F           *fHisPy3vP1;     //-> delta ky  P1outer-P1inner - common parabola
2582   //   TH2F           *fHisPy2vP2In;   //-> Curv  P2inner - parabola
2583   //   TH2F           *fHisPy2vP2Out;  //-> Curv  P2outer - parabola
2584   //   TH2F           *fHisPy3vP2IO;   //-> Curv  P2outerinner - common parabola
2585   fHisPy1vP0   = new TH2F("HisPy1vP0",   "HisPy1vP0",336,0,336,500,-0.3,0.3);
2586   fHisPy2vP0   = new TH2F("HisPy2vP0",   "HisPy2vP0",336,0,336,500,-0.3,0.3);
2587   fHisPy3vP0   = new TH2F("HisPy3vP0",   "HisPy3vP0",336,0,336,500,-0.3,0.3);
2588   fHisPy1vP1   = new TH2F("HisPy1vP1",   "HisPy1vP1",336,0,336,500,-0.01,0.01);
2589   fHisPy2vP1   = new TH2F("HisPy2vP1",   "HisPy2vP1",336,0,336,500,-0.01,0.01);
2590   fHisPy3vP1   = new TH2F("HisPy3vP1",   "HisPy3vP1",336,0,336,500,-0.01,0.01);
2591   fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
2592   fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
2593   fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
2594   //
2595   //
2596   //   TH2F           *fHisPz1vP0;     //-> delta z   P0outer-P0inner - line
2597   //   TH2F           *fHisPz2vP0;     //-> delta z   P0outer-P0inner - parabola
2598   //   TH2F           *fHisPz3vP0;     //-> delta z   P0outer-P0inner - common parabola
2599   //   TH2F           *fHisPz1vP1;     //-> delta kz  P1outer-P1inner - line
2600   //   TH2F           *fHisPz2vP1;     //-> delta kz  P1outer-P1inner - parabola
2601   //   TH2F           *fHisPz3vP1;     //-> delta kz  P1outer-P1inner - common parabola
2602   //   TH2F           *fHisPz2vP2In;   //-> Curv  P2inner - parabola
2603   //   TH2F           *fHisPz2vP2Out;  //-> Curv  P2outer - parabola
2604   //   TH2F           *fHisPz3vP2IO;   //-> Curv  P2outerinner - common parabola
2605   fHisPz1vP0   = new TH2F("HisPz1vP0",   "HisPz1vP0",336,0,336,500,-0.3,0.3);
2606   fHisPz2vP0   = new TH2F("HisPz2vP0",   "HisPz2vP0",336,0,336,500,-0.3,0.3);
2607   fHisPz3vP0   = new TH2F("HisPz3vP0",   "HisPz3vP0",336,0,336,500,-0.3,0.3);
2608   fHisPz1vP1   = new TH2F("HisPz1vP1",   "HisPz1vP1",336,0,336,500,-0.01,0.01);
2609   fHisPz2vP1   = new TH2F("HisPz2vP1",   "HisPz2vP1",336,0,336,500,-0.01,0.01);
2610   fHisPz3vP1   = new TH2F("HisPz3vP1",   "HisPz3vP1",336,0,336,500,-0.01,0.01);
2611   fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
2612   fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
2613   fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
2614
2615   fHisNclIn->SetDirectory(0);      //->Number of clusters inner
2616   fHisNclOut->SetDirectory(0);     //->Number of clusters outer
2617   fHisNclIO->SetDirectory(0);      //->Number of cluster inner outer
2618   fHisLclIn->SetDirectory(0);      //->Level arm inner
2619   fHisLclOut->SetDirectory(0);     //->Level arm outer
2620   fHisLclIO->SetDirectory(0);      //->Number of cluster inner outer
2621   fHisdEdx->SetDirectory(0);      //->Number of cluster inner outer
2622   fHisdZfit->SetDirectory(0);      //->Number of cluster inner outer
2623   //
2624   //
2625   fHisChi2YIn1->SetDirectory(0);      //->chi2 y inner - line
2626   fHisChi2YOut1->SetDirectory(0);     //->chi2 y inner - line
2627   fHisChi2YIn2->SetDirectory(0);      //->chi2 y inner - parabola
2628   fHisChi2YOut2->SetDirectory(0);     //->chi2 y inner - parabola
2629   fHisChi2YIO1->SetDirectory(0);      //->chi2 y IO    - common
2630   fHisChi2ZIn1->SetDirectory(0);      //->chi2 z inner - line
2631   fHisChi2ZOut1->SetDirectory(0);     //->chi2 z inner - line
2632   fHisChi2ZIn2->SetDirectory(0);      //->chi2 z inner - parabola
2633   fHisChi2ZOut2->SetDirectory(0);     //->chi2 z inner - parabola
2634   fHisChi2ZIO1->SetDirectory(0);      //->chi2 z IO    - common
2635   //
2636   //
2637   fHisPy1vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - line
2638   fHisPy2vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - parabola
2639   fHisPy3vP0->SetDirectory(0);     //-> delta y   P0outer-P0inner - common parabola
2640   fHisPy1vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - line
2641   fHisPy2vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - parabola
2642   fHisPy3vP1->SetDirectory(0);     //-> delta ky  P1outer-P1inner - common parabola
2643   fHisPy2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
2644   fHisPy2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
2645   fHisPy3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
2646   //
2647   //
2648   fHisPz1vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - line
2649   fHisPz2vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - parabola
2650   fHisPz3vP0->SetDirectory(0);     //-> delta z   P0outer-P0inner - common parabola
2651   fHisPz1vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - line
2652   fHisPz2vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - parabola
2653   fHisPz3vP1->SetDirectory(0);     //-> delta kz  P1outer-P1inner - common parabola
2654   fHisPz2vP2In->SetDirectory(0);   //-> Curv  P2inner - parabola
2655   fHisPz2vP2Out->SetDirectory(0);  //-> Curv  P2outer - parabola
2656   fHisPz3vP2IO->SetDirectory(0);   //-> Curv  P2outerinner - common parabola
2657
2658
2659   //
2660   //
2661   //
2662   for (Int_t id=0; id<336;id++){
2663     TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
2664     TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
2665     TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
2666     TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
2667     TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2668     TH2F *profzabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
2669     //    TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
2670     //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
2671     if (!profy){
2672       profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
2673       profy->SetDirectory(0);
2674       fDeltaYres.AddAt(profy,id);
2675       profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
2676       profy2->SetDirectory(0);
2677       fDeltaYres2.AddAt(profy2,id);
2678       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
2679       profyabs->SetDirectory(0);
2680       fDeltaYresAbs.AddAt(profyabs,id);
2681       //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
2682       //profy3->SetDirectory(0);
2683       //fDeltaYres3.AddAt(profy3,id);
2684     }
2685     if (!profz){
2686       profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
2687       profz->SetDirectory(0);
2688       fDeltaZres.AddAt(profz,id);
2689       profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
2690       profz2->SetDirectory(0);
2691       fDeltaZres2.AddAt(profz2,id);
2692       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
2693       profzabs->SetDirectory(0);
2694       fDeltaZresAbs.AddAt(profzabs,id);
2695       //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
2696       //profz3->SetDirectory(0);
2697       //fDeltaZres3.AddAt(profz3,id);
2698     }
2699   }
2700   //
2701   //
2702   for (Int_t id=0; id<336;id++){
2703     TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
2704     TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
2705     TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
2706
2707     TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
2708     TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
2709     TH1F * hisSignal = (TH1F*)fSignals.At(id);
2710
2711     if (!hisdz){
2712       hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
2713       hisdz->SetDirectory(0);
2714       fDeltaZ.AddAt(hisdz,id);
2715
2716       hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
2717       hisP3->SetDirectory(0);
2718       fDeltaP3.AddAt(hisP3,id);
2719       //
2720       hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
2721       hisP4->SetDirectory(0);
2722       fDeltaP4.AddAt(hisP4,id);
2723
2724       //
2725       hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
2726       hisdphi->SetDirectory(0);
2727       fDeltaPhi.AddAt(hisdphi,id);
2728       //
2729       hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
2730       hisdphiP->SetDirectory(0);
2731       fDeltaPhiP.AddAt(hisdphiP,id);
2732       hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
2733       hisSignal->SetDirectory(0);
2734       fSignals.AddAt(hisSignal,id);
2735     }
2736   }
2737
2738   SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
2739   SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
2740   SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
2741   SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
2742   //
2743   // Make THnSparse
2744   //
2745   //                   id   side   rod  bundle  beam  dP0  dP1  dP2  dP3  dP4 ncl dEdx
2746   Int_t binsLaser[12]= {336,  //id
2747                         2,    //side
2748                         6,    //rod
2749                         4,    //bundle
2750                         7,    //beam
2751                         300,  //dP0
2752                         300,  //dP1
2753                         300,  //dP2
2754                         300,  //dP3
2755                         300,  //dP4
2756                         80,   //ncl
2757                         50};   //dEdx
2758   Double_t xminLaser[12]= {0,  //id
2759                         0,     //side
2760                         0,     //rod
2761                         0,     //bundle
2762                         0,     //beam
2763                         -1,    //dP0
2764                         -1,    //dP1
2765                         -0.01,   //dP2
2766                         -0.01,  //dP3
2767                         -0.1,  //dP4
2768                         0,   //ncl
2769                         0};   //sqrt dEdx
2770   Double_t xmaxLaser[12]= {336,  //id
2771                         2,    //side
2772                         6,    //rod
2773                         4,    //bundle
2774                         7,    //beam
2775                         1,  //dP0
2776                         1,  //dP1
2777                         0.01,  //dP2
2778                         0.01,  //dP3
2779                         0.1,  //dP4
2780                         160,   //ncl
2781                         40};   //sqrt dEdx
2782
2783   TString nameLaser[12]= {"id",
2784                           "side",
2785                           "rod",
2786                           "bundle",
2787                           "beam",
2788                           "dP0",
2789                           "dP1",
2790                           "dP2",
2791                           "dP3",
2792                           "dP4",
2793                           "ncl",
2794                           "sqrt dEdx"};
2795   TString titleLaser[12]= {"id",
2796                           "side",
2797                           "rod",
2798                           "bundle",
2799                           "beam",
2800                           "#Delta_{P0}",
2801                           "#Delta_{P1}",
2802                           "#Delta_{P2}",
2803                           "#Delta_{P3}",
2804                           "#Delta_{P4}",
2805                           "N_{cl}",
2806                           "#sqrt{dEdx}"};
2807   fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
2808   for (Int_t iaxis=1; iaxis<12; iaxis++){
2809     fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
2810     fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
2811   }
2812 }
2813
2814 void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
2815   //
2816   // Merge content of histograms
2817   //
2818   // Only first histogram is checked - all other should be the same
2819   if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
2820
2821   if (!laser->fHisNclIn) return;  // empty histograms
2822   if (!fHisNclIn) MakeFitHistos();
2823   //
2824   fHisNclIn->Add(laser->fHisNclIn  );      //->Number of clusters inner
2825   fHisNclOut->Add(laser->fHisNclOut  );     //->Number of clusters outer
2826   fHisNclIO->Add(laser->fHisNclIO  );      //->Number of cluster inner outer
2827   fHisLclIn->Add(laser->fHisLclIn  );      //->Level arm inner
2828   fHisLclOut->Add(laser->fHisLclOut  );     //->Level arm outer
2829   fHisLclIO->Add(laser->fHisLclIO  );      //->Number of cluster inner outer
2830   fHisdEdx->Add(laser->fHisdEdx  );      //->dedx
2831   fHisdZfit->Add(laser->fHisdZfit  );      //->dedx
2832   //
2833   //
2834   fHisChi2YIn1->Add(laser->fHisChi2YIn1  );      //->chi2 y inner - line
2835   fHisChi2YOut1->Add(laser->fHisChi2YOut1  );     //->chi2 y inner - line
2836   fHisChi2YIn2->Add(laser->fHisChi2YIn2  );      //->chi2 y inner - parabola
2837   fHisChi2YOut2->Add(laser->fHisChi2YOut2  );     //->chi2 y inner - parabola
2838   fHisChi2YIO1->Add(laser->fHisChi2YIO1  );      //->chi2 y IO    - common
2839   fHisChi2ZIn1->Add(laser->fHisChi2ZIn1  );      //->chi2 z inner - line
2840   fHisChi2ZOut1->Add(laser->fHisChi2ZOut1  );     //->chi2 z inner - line
2841   fHisChi2ZIn2->Add(laser->fHisChi2ZIn2  );      //->chi2 z inner - parabola
2842   fHisChi2ZOut2->Add(laser->fHisChi2ZOut2  );     //->chi2 z inner - parabola
2843   fHisChi2ZIO1->Add(laser->fHisChi2ZIO1  );      //->chi2 z IO    - common
2844   //
2845   //
2846   fHisPy1vP0->Add(laser->fHisPy1vP0  );     //-> delta y   P0outer-P0inner - line
2847   fHisPy2vP0->Add(laser->fHisPy2vP0  );     //-> delta y   P0outer-P0inner - parabola
2848   fHisPy3vP0->Add(laser->fHisPy3vP0  );     //-> delta y   P0outer-P0inner - common parabola
2849   fHisPy1vP1->Add(laser->fHisPy1vP1  );     //-> delta ky  P1outer-P1inner - line
2850   fHisPy2vP1->Add(laser->fHisPy2vP1  );     //-> delta ky  P1outer-P1inner - parabola
2851   fHisPy3vP1->Add(laser->fHisPy3vP1  );     //-> delta ky  P1outer-P1inner - common parabola
2852   fHisPy2vP2In->Add(laser-> fHisPy2vP2In  );   //-> Curv  P2inner - parabola
2853   fHisPy2vP2Out->Add(laser->fHisPy2vP2Out  );  //-> Curv  P2outer - parabola
2854   fHisPy3vP2IO->Add(laser->fHisPy3vP2IO  );   //-> Curv  P2outerinner - common parabola
2855   //
2856   //
2857   fHisPz1vP0->Add(laser->fHisPz1vP0  );     //-> delta z   P0outer-P0inner - line
2858   fHisPz2vP0->Add(laser->fHisPz2vP0  );     //-> delta z   P0outer-P0inner - parabola
2859   fHisPz3vP0->Add(laser->fHisPz3vP0  );     //-> delta z   P0outer-P0inner - common parabola
2860   fHisPz1vP1->Add(laser->fHisPz1vP1  );     //-> delta kz  P1outer-P1inner - line
2861   fHisPz2vP1->Add(laser->fHisPz2vP1  );     //-> delta kz  P1outer-P1inner - parabola
2862   fHisPz3vP1->Add(laser->fHisPz3vP1  );     //-> delta kz  P1outer-P1inner - common parabola
2863   fHisPz2vP2In->Add(laser->fHisPz2vP2In  );   //-> Curv  P2inner - parabola
2864   fHisPz2vP2Out->Add(laser->fHisPz2vP2Out  );  //-> Curv  P2outer - parabola
2865   fHisPz3vP2IO->Add(laser->fHisPz3vP2IO  );   //-> Curv  P2outerinner - common parabola
2866   //
2867   //
2868   //
2869
2870
2871
2872
2873 }
2874
2875
2876
2877
2878 void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
2879   //
2880   // Dump fit information - collect information from the streamers
2881   //
2882   /*
2883     TChain * chainFit=0;
2884     TChain * chainTrack=0;
2885     TChain * chain=0;
2886     //
2887     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
2888     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
2889     AliXRDPROOFtoolkit tool;
2890     chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
2891     chainTrack->Lookup();
2892     chainTrack->SetProof(kTRUE);
2893     chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
2894     chainDrift->Lookup();
2895     chainDrift->SetProof(kTRUE);
2896
2897     chain = tool.MakeChain("laser.txt","Residuals",0,10200);
2898     chain->Lookup();
2899     chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
2900     chainFit->Lookup();
2901     chainFit->SetProof(kTRUE);
2902     chain->SetProof(kTRUE);
2903     AliTPCLaserTrack::LoadTracks();
2904     //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
2905
2906   */
2907   //
2908   // Fit cuts
2909   //
2910   TCut cutP3("abs(Tr.fP[3])<0.1");
2911   TCut cutP4("abs(Tr.fP[4])<0.5");
2912   TCut cutPx = cutP3+cutP4;
2913   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
2914   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
2915   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
2916   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
2917   //
2918   TCut cutdEdx("sqrt(dEdx)>3");
2919   TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
2920   TCut cutN("nclO>20&&nclI>20");
2921   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
2922   //
2923   // Cluster cuts
2924   //
2925   TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
2926   TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
2927   TCut cutClX("abs(Cl[].fX)>10");
2928   TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
2929   TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
2930   TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
2931   TCut cutQ("sqrt(Cl[].fMax)>4");
2932   TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
2933
2934
2935   TH1F * phisAl     = 0;
2936   TH1F * phisAccept = 0;
2937   TH1F * phisOut    = 0;
2938   TProfile * pdEdx  = 0;
2939
2940   TProfile * pP0    = 0;
2941   TProfile * pP1    = 0;
2942   TProfile * pP2    = 0;
2943   TProfile * pP3    = 0;
2944   TProfile * pP4    = 0;
2945   //
2946   TProfile * pNclI  = 0;
2947   TProfile * pNclO  = 0;
2948   //
2949   TProfile * pchi2YIn   =0;
2950   TProfile * pchi2ZIn   =0;
2951   TProfile * pchi2YOut  =0;
2952   TProfile * pchi2ZOut  =0;
2953   TProfile * pchi2YInOut =0;
2954   TProfile * pchi2ZInOut =0;;
2955   // laser counters
2956   chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
2957   phisAl = (TH1F*)gROOT->FindObject("hisAl");
2958   chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
2959   phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
2960   chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
2961   phisOut = (TH1F*)gROOT->FindObject("hisOut");
2962   //
2963   chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
2964   pdEdx   = (TProfile*)gROOT->FindObject("hdEdx");
2965   // track param
2966   //
2967   chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2968   pP0   = (TProfile*)gROOT->FindObject("hP0");
2969   chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2970   pP1   = (TProfile*)gROOT->FindObject("hP1");
2971   chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2972   pP2   = (TProfile*)gROOT->FindObject("hP2");
2973   chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2974   pP3   = (TProfile*)gROOT->FindObject("hP3");
2975   chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
2976   pP4   = (TProfile*)gROOT->FindObject("hP4");
2977   //
2978   chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
2979   pNclI   = (TProfile*)gROOT->FindObject("hNclI");
2980   chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
2981   pNclO   = (TProfile*)gROOT->FindObject("hNclO");
2982   //
2983   //
2984   chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
2985   pchi2YIn   = (TProfile*)gROOT->FindObject("hChi2YIn");
2986   chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
2987   pchi2YOut   = (TProfile*)gROOT->FindObject("hChi2YOut");
2988   chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
2989   pchi2YInOut   = (TProfile*)gROOT->FindObject("hChi2YInOut");
2990   chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
2991   pchi2ZIn   = (TProfile*)gROOT->FindObject("hChi2ZIn");
2992   chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
2993   pchi2ZOut   = (TProfile*)gROOT->FindObject("hChi2ZOut");
2994   chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
2995   pchi2ZInOut   = (TProfile*)gROOT->FindObject("hChi2ZInOut");
2996   //
2997   // second derivatives
2998   //
2999   TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
3000   chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
3001   TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
3002   chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
3003   TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
3004   chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
3005   //
3006   phisPy2In->FitSlicesY();
3007   TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
3008   TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
3009   TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
3010   //
3011   phisPy2Out->FitSlicesY();
3012   TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
3013   TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
3014   TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
3015   //
3016   phisPy2InOut->FitSlicesY();
3017   TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
3018   TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
3019   TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
3020   //
3021   TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
3022   chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
3023   TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
3024   chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
3025   TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
3026   chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
3027   //
3028   phisPz2In->FitSlicesY();
3029   TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
3030   TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
3031   TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
3032   //
3033   phisPz2Out->FitSlicesY();
3034   TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
3035   TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
3036   TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
3037   //
3038   phisPz2InOut->FitSlicesY();
3039   TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
3040   TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
3041   TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
3042   //
3043   //
3044   //
3045
3046
3047   {
3048     TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
3049     for (Int_t ilaser=0; ilaser<336; ilaser++){
3050       Float_t all=phisAl->GetBinContent(ilaser+1);
3051       Float_t accept=phisAccept->GetBinContent(ilaser+1);
3052       Float_t out=phisOut->GetBinContent(ilaser+1);
3053       Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
3054       Float_t mP0 = pP0->GetBinContent(ilaser+1);
3055       Float_t mP1 = pP1->GetBinContent(ilaser+1);
3056       Float_t mP2 = pP2->GetBinContent(ilaser+1);
3057       Float_t mP3 = pP3->GetBinContent(ilaser+1);
3058       Float_t mP4 = pP4->GetBinContent(ilaser+1);
3059
3060
3061       Float_t nclI  = pNclI->GetBinContent(ilaser+1);
3062       Float_t nclO  = pNclO->GetBinContent(ilaser+1);
3063       //
3064       Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
3065       Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
3066       Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
3067       Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
3068       Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
3069       Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
3070       //
3071       Float_t entriesPy2In  = phisPy2InEntries->GetBinContent(ilaser+1);
3072       Float_t meanPy2In     = phisPy2InMean->GetBinContent(ilaser+1);
3073       Float_t sigmaPy2In    = phisPy2InSigma->GetBinContent(ilaser+1);
3074       //
3075       Float_t entriesPy2Out  = phisPy2OutEntries->GetBinContent(ilaser+1);
3076       Float_t meanPy2Out     = phisPy2OutMean->GetBinContent(ilaser+1);
3077       Float_t sigmaPy2Out    = phisPy2OutSigma->GetBinContent(ilaser+1);
3078       //
3079       Float_t entriesPy2InOut  = phisPy2InOutEntries->GetBinContent(ilaser+1);
3080       Float_t meanPy2InOut     = phisPy2InOutMean->GetBinContent(ilaser+1);
3081       Float_t sigmaPy2InOut    = phisPy2InOutSigma->GetBinContent(ilaser+1);
3082       //
3083       Float_t entriesPz2In  = phisPz2InEntries->GetBinContent(ilaser+1);
3084       Float_t meanPz2In     = phisPz2InMean->GetBinContent(ilaser+1);
3085       Float_t sigmaPz2In    = phisPz2InSigma->GetBinContent(ilaser+1);
3086       //
3087       Float_t entriesPz2Out  = phisPz2OutEntries->GetBinContent(ilaser+1);
3088       Float_t meanPz2Out     = phisPz2OutMean->GetBinContent(ilaser+1);
3089       Float_t sigmaPz2Out    = phisPz2OutSigma->GetBinContent(ilaser+1);
3090       //
3091       Float_t entriesPz2InOut  = phisPz2InOutEntries->GetBinContent(ilaser+1);
3092       Float_t meanPz2InOut     = phisPz2InOutMean->GetBinContent(ilaser+1);
3093       Float_t sigmaPz2InOut    = phisPz2InOutSigma->GetBinContent(ilaser+1);
3094
3095       AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
3096       (*pcstream)<<"Scan"<<
3097         "Run="<<id<<
3098         "LTr.="<<ltrp<<
3099         "all="<<all<<
3100         "accept="<<accept<<
3101         "out="<<out<<
3102         "sdedx="<<sdedx<<
3103         "mP0="<<mP0<<
3104         "mP1="<<mP1<<
3105         "mP2="<<mP2<<
3106         "mP3="<<mP3<<
3107         "mP4="<<mP4<<
3108         "nclI="<<nclI<<
3109         "nclO="<<nclO<<
3110         "chi2YIn="<<chi2YIn<<
3111         "chi2YOut="<<chi2YOut<<
3112         "chi2YInOut="<<chi2YInOut<<
3113         "chi2ZIn="<<chi2ZIn<<
3114         "chi2ZOut="<<chi2ZOut<<
3115         "chi2ZInOut="<<chi2ZInOut<<
3116         //
3117         "nPy2In="<<entriesPy2In<<
3118         "mPy2In="<<meanPy2In<<
3119         "sPy2In="<<sigmaPy2In<<
3120         //
3121         "nPy2Out="<<entriesPy2Out<<
3122         "mPy2Out="<<meanPy2Out<<
3123         "sPy2Out="<<sigmaPy2Out<<
3124         //
3125         "nPy2InOut="<<entriesPy2InOut<<
3126         "mPy2InOut="<<meanPy2InOut<<
3127         "sPy2InOut="<<sigmaPy2InOut<<
3128         //
3129         "nPz2In="<<entriesPz2In<<
3130         "mPz2In="<<meanPz2In<<
3131         "sPz2In="<<sigmaPz2In<<
3132         //
3133         "nPz2Out="<<entriesPz2Out<<
3134         "mPz2Out="<<meanPz2Out<<
3135         "sPz2Out="<<sigmaPz2Out<<
3136         //
3137         "nPz2InOut="<<entriesPz2InOut<<
3138         "mPz2InOut="<<meanPz2InOut<<
3139         "sPz2InOut="<<sigmaPz2InOut<<
3140         "\n";
3141     }
3142
3143     delete pcstream;
3144   }
3145 }
3146
3147 void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
3148                                          TVectorD& meanSlope,
3149                                          TVectorD& sectorArray,
3150                                          Int_t option)
3151 {
3152   // This method should ideally go in AliTPCLaser
3153   // option == 0 (pads outer - closest to beam)
3154   // option == 1 (pads inner)
3155   // option == 2 (time outer)
3156   // option == 3 (time inner)
3157   Int_t nFailures = 0;
3158
3159   for(Int_t id = 0; id < 336; id++) {
3160
3161     if (!AliTPCLaserTrack::GetTracks())
3162       AliTPCLaserTrack::LoadTracks();
3163     AliTPCLaserTrack *ltrp =
3164       (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
3165
3166     AliExternalTrackParam trackParam(*ltrp);
3167
3168     Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
3169     if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
3170       deltaangle = 30.0; // inner sector is 1 sector further away from mirror
3171
3172     Double_t angle = trackParam.GetAlpha();
3173     if(angle<0)
3174       angle += 2*TMath::Pi();
3175     if(trackParam.GetSnp()>0) // track points to sector "before"
3176       angle -= deltaangle*TMath::DegToRad();
3177     else // track points to sector "after"
3178       angle += deltaangle*TMath::DegToRad();
3179
3180     Bool_t success = trackParam.Rotate(angle);
3181
3182     if(!success) {
3183       //      cout << "WARNING: Rotate failed for ID: " << id << endl;
3184       nFailures++;
3185     }
3186
3187     angle *= TMath::RadToDeg();
3188
3189     Int_t sector = TMath::Nint((angle-10.0)/20.0);
3190     if(sector<0)
3191       sector += 18;
3192     else if(sector>=18)
3193       sector -= 18;
3194     if(ltrp->GetSide()==1) // C side
3195       sector += 18;
3196     sectorArray[id] = sector;
3197
3198     const Double_t x0 = 0;
3199
3200     Double_t slopey  = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
3201     Double_t slopez  = trackParam.GetTgl();
3202     //  One needs a factor sqrt(1+slopey**2) to take into account the
3203     //  longer path length
3204     slopez *= TMath::Sqrt(1.0 + slopey*slopey);
3205     if(fInverseSlopeZ) // wrong sign in database, should be fixed there
3206       slopez *= -1;
3207     //    Double_t offsetz = trackParam.GetZ();
3208     Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
3209     Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
3210     if(option==2 || option==3) {
3211       meanOffset[id] = offsetz; meanSlope[id] = slopez;
3212     } else {
3213       meanOffset[id] = offsety; meanSlope[id] = slopey;
3214     }
3215   }
3216
3217   if(nFailures>0)
3218     AliWarning(Form("Rotate method failed %d times", nFailures));
3219 }
3220
3221
3222   /*
3223     TFile f("vscan.root");
3224    */
3225
3226   /*
3227     pad binning effect
3228    chain->Draw("Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","",10000);
3229    //
3230    chain->Draw("Cl[].fY-TrYpol1.fElements:Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
3231    //
3232
3233 chain->Draw("Cl.fY-TrYpol1.fElements-AliTPCClusterParam::SPosCorrection(0,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fPad-int(Cl[].fPad)",cutA+cutCl+"Cl[].fZ>0&&Cl[].fPad>1","prof",10000);
3234
3235
3236 chain->Draw("Cl[].fZ-TrZpol1.fElements-0*AliTPCClusterParam::SPosCorrection(1,1,Cl[].fPad,Cl[].fTimeBin,Cl[].fZ,Cl[].fSigmaY2,Cl[].fSigmaZ2,Cl[].fMax):Cl[].fTimeBin-int(Cl[].fTimeBin)",cutA+cutCl+"Cl[].fZ>0","prof",10000)
3237
3238   */
3239
3240
3241
3242
3243
3244   /*
3245   // check edge effects
3246   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
3247   //
3248   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
3249
3250   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
3251
3252
3253
3254   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
3255   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
3256
3257 */
3258
3259
3260
3261
3262
3263   /*
3264     Edge y effect
3265
3266     dedge = sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))
3267
3268
3269     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):pi/18-abs(Cl.fY/Cl.fX)>>hisYdphi(100,0,0.03)",""+cutA+cutCl,"prof",10000)
3270
3271     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdy(100,0,5)",""+cutA+cutCl,"prof",10000)
3272
3273
3274
3275
3276
3277     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyIROC(100,0,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
3278
3279     chain->Draw("sign(Cl.fY)*(Cl.fY-TrYpol1.fElements):Cl.fX*(pi/18-abs(Cl.fY/Cl.fX))>>hisYdyOROC(100,0,5)","Cl.fDetector>36"+cutA+cutCl,"prof",100000)
3280
3281
3282
3283     chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>his(100,-5,5)",""+cutA+cutCl,"prof",100000)
3284
3285     chain->Draw("Cl.fY-TrYpol1.fElements:sign(Cl.fY)*(Cl.fX*tan(pi/18)-abs(Cl.fY))>>hisdyInner(100,-5,5)","Cl.fDetector<36"+cutA+cutCl,"prof",100000)
3286
3287
3288
3289 */
3290
3291
3292 /*
3293
3294 chainFit->Draw("yPol2Out.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutDY,"prof")
3295
3296 chainFit->Draw("yPol2In.fElements[2]*64*64/4.:LTr.fP[2]","nclI>20&&LTr.fP[1]<0"+cutA+cutDY,"prof")
3297
3298
3299
3300 chainFit->Draw("LTr.fId","nclI>10",100000)
3301
3302 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>his(350,0,350,100,-0.002,0.002)","nclI>20","")
3303
3304 chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>hisPy2In0(350,0,350,100,-0.002,0.002)","nclI>20","");
3305
3306 TH2 * phisPy2In = (TH2*) gROOT->FindObject("hisPy2In0")
3307
3308 */
3309
3310
3311
3312
3313
3314
3315 /*
3316  gSystem->Load("libSTAT.so")
3317  TStatToolkit toolkit;
3318  Double_t chi2;
3319  TVectorD fitParam;
3320  TMatrixD covMatrix;
3321  Int_t npoints;
3322
3323  TCut cutA("entries>2&&pphi2<3&&abs(gphiP1-pphiP0)<0.003&&abs(gz1)<6");
3324
3325
3326 TString fstring="";
3327 //
3328 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz++";                               //1
3329 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*LTr.fP[2]++";                     //2
3330 fstring+="(abs(LTr.fP[1]/250)^1-1)*bz++";                               //3
3331 fstring+="(abs(LTr.fP[1]/250)-1)*bz*LTr.fP[2]++";                       //4
3332 //
3333 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))++"            //5
3334 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"  //6
3335 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))++"              //7
3336 fstring+="(abs(LTr.fP[1]/250)-1)*bz*sin(atan2(lx1,lx0))*LTr.fP[2]++"    //8
3337 //
3338 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))++"            //9
3339 fstring+="(abs(LTr.fP[1]/250)^3-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"  //10
3340 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))++"              //11
3341 fstring+="(abs(LTr.fP[1]/250)-1)*bz*cos(atan2(lx1,lx0))*LTr.fP[2]++"    //12
3342
3343
3344
3345
3346  TString *strq0 = toolkit.FitPlane(treeT,"gphi1-pphi0",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
3347
3348  treeT->SetAlias("fit",strq0->Data());
3349
3350
3351  TString *strqP = toolkit.FitPlane(treeT,"1000*(gphiP1-pphiP0)",fstring->Data(), "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix);
3352
3353  treeT->SetAlias("fitP",strqP->Data());
3354
3355
3356  TString *strqDrift = toolkit.FitPlane(treeT,"gz1","LTr.fP[1]++(1-2*(fSide==1))++lx1", cutA, chi2,npoints,fitParam,covMatrix);
3357  treeT->SetAlias("fitD",strqDrift->Data());
3358
3359
3360 treeT->Draw("fit:LTr.fP[1]","abs(bz+0.4)<0.05"+cutA,"");
3361 {
3362 for (Int_t i=0; i<6;i++){
3363 treeT->SetLineColor(i+2);
3364 treeT->SetMarkerSize(1);
3365 treeT->SetMarkerStyle(22+i);
3366 treeT->SetMarkerColor(i+2);
3367
3368 treeT->Draw("fit:LTr.fP[1]",Form("abs(bz+0.4)<0.05&fRod==%d",i)+cutA,"same");
3369 }
3370 }
3371  */
3372
3373
3374
3375 /*
3376   TTree * tree = (TTree*)f.Get("FitModels");
3377
3378   TEventList listLFit0("listLFit0","listLFit0");
3379   TEventList listLFit1("listLFit1","listLFit1");
3380   tree->Draw(">>listLFit0","seed.fdEdx<200&&seed.fdEdx>40");
3381   tree->SetEventList(&listLFit0);
3382
3383
3384
3385
3386   gSystem->Load("libSTAT.so")
3387   TStatToolkit toolkit;
3388   Double_t chi2;
3389   TVectorD fitParam;
3390   TMatrixD covMatrix;
3391   Int_t npoints;
3392
3393   chain->SetAlias("dp","((Cl.fPad-int(Cl.fPad))*pi)");
3394   chain->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin))*pi)");
3395
3396
3397   TString fstring="";
3398   fstring+="cos(dp)++";
3399   fstring+="sin(dp)++";
3400   fstring+="cos(dt)++";
3401   fstring+="sin(dt)++";
3402
3403   TString *str = toolkit.FitPlane(chain,"Cl.fZ-TrZInOut.fElements",fstring->Data(), "Cl.fDetector>35", chi2,npoints,fitParam,covMatrix,-1,0,200);
3404
3405
3406
3407 */
3408
3409
3410
3411 /*
3412   Edge effects
3413   //
3414   //
3415   //
3416   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
3417   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
3418   AliXRDPROOFtoolkit tool;
3419   TChain * chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
3420   chainTrack->Lookup();
3421   chainTrack->SetProof(kTRUE);
3422
3423   TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
3424   chain->Lookup();
3425   TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
3426   chainFit->Lookup();
3427   chainFit->SetProof(kTRUE);
3428   chain->SetProof(kTRUE);
3429   //
3430   // Fit cuts
3431   //
3432   TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<10");
3433   TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<10");
3434   TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<10");
3435   TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<10");
3436   //
3437   TCut cutdEdx("sqrt(dEdx)<30&&sqrt(dEdx)>3");
3438   TCut cutDY("abs(yPol2In.fElements[2]*nclO*nclO/4.)<3")
3439   TCut cutN("nclO>20&&nclI>20");
3440   TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx;
3441   //
3442   // Cluster cuts
3443   //
3444   TCut cutClY("abs(Cl.fY-TrYpol2.fElements)<0.2");
3445   TCut cutClZ("abs(Cl.fZ-TrZpol2.fElements)<0.4");
3446   TCut cutClX("abs(Cl.fX)>10");
3447   TCut cutE("abs(Cl.fY/Cl.fX)<0.14");
3448   TCut cutCl=cutClY+cutClZ+cutClX;
3449
3450
3451   // check edge effects
3452   chain->Draw("Cl.fY-TrYpol1.fElements:Cl.fY/Cl.fX",""+cutA+cutCl,"prof",10000)
3453   //
3454   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fZ>0"+cutA+cutCl+cutE,"prof",100000)
3455
3456   chain->Draw("Cl.fY-TrYpol2.fElements:Cl.fPad-int(Cl.fPad)","Cl.fX>80&&Cl.fZ>0&&Cl.fDetector>35"+cutA+cutCl+cutE,"prof",100000)
3457
3458
3459
3460   chainFit->Draw("yInOut.fElements[4]:LTr.fP[2]","LTr.fP[1]<0"+cutA,"prof",1000);
3461   chainFit->Draw("yPol2In.fElements[2]*90*90/4.:LTr.fP[2]","nclO>40&&LTr.fP[1]<0"+cutA+cutD,"prof")
3462
3463 */
3464