TPC/AliTPCcalibLaser.cxx

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