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