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