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