**************************************************************************/
/*
- laser track clasification;
+ //
+ // FUNCTIONALITY:
+ //
+ // 1. The laser track is associated with the mirror
+ // see function FindMirror
+ //
+ // 2. The laser track is accepted for the analysis under certain condition
+ // (see function Accpet laser)
+ //
+ // 3. The drift velocity and jitter is calculated event by event
+ // (see function drift velocity)
+ //
+ // 4. The tracks are refitted at different sectors
+ // Fit model
+ // 4.a) line
+ // 4.b) parabola
+ // 4.c) parabola with common P2 for inner and outer
+ //
+ // To make laser scan the user interaction neccessary
+ //
+ .x ~/NimStyle.C
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libTPCcalib");
+ TFile fcalib("CalibObjectsTrain2.root");
+ AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)fcalib->Get("laserTPC");
+ laser->DumpMeanInfo(run)
+ TFile fmean("laserMean.root")
+ //
+ // laser track clasification;
+ //
TCut cutT("cutT","abs(Tr.fP[3])<0.06");
- TCut cutPt("cutPt","abs(Tr.fP[4])<0.1");
- TCut cutN("cutN","fTPCncls>100");
+ TCut cutN("cutN","fTPCncls>70");
+ TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03")
+ TCut cutA = cutT+cutPt+cutP;
+ TChain * chainTrL = tool.MakeChain("laser.txt","Track",0,10200);
- TCut cutFi("cutZB","");
- TCut cutFi("cutFi","abs((180*atan2(x1,x0)/pi-20)%40)<5");
+ //
+ //
+ // Analyze LASER scan
+ //
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+ TChain * chain = tool.MakeChain("laserScan.txt","Mean",0,10200);
+ chain->Lookup();
+ AliTPCcalibLaser::DumpScanInfo(chain)
+ TFile fscan("laserScan.root")
+ TTree * treeT = (TTree*)fscan.Get("Mean")
+ //
+ // Analyze laser
+ //
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+ AliXRDPROOFtoolkit tool;
+ AliXRDPROOFtoolkit::FilterList("laserDebug.list","* driftvN",1)
+ TChain * chainDrift = tool.MakeChainRandom("laser.txt.Good","driftv",0,50);
+ chainDrift->Lookup();
+ TChain * chainDriftN = tool.MakeChainRandom("laserDebug.list.Good","driftvN",0,300);
+ chainDriftN->Lookup();
+
+
+ TChain * chain = tool.MakeChain("laser.txt","Residuals",0,10200);
+ chain->Lookup();
+ TChain * chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
+ chainFit->Lookup();
+ TChain * chainTrack = tool.MakeChainRandom("laser.txt","Track",0,30);
+ chainTrack->Lookup();
+
*/
#include "TLinearFitter.h"
#include "AliTPCcalibLaser.h"
#include "AliExternalTrackParam.h"
+#include "AliESDEvent.h"
+#include "AliESDfriend.h"
#include "AliESDtrack.h"
#include "AliTPCTracklet.h"
#include "TH1D.h"
+#include "TH1F.h"
+#include "TProfile.h"
#include "TVectorD.h"
#include "TTreeStream.h"
#include "TFile.h"
#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliTracker.h"
+#include "AliLog.h"
#include "TClonesArray.h"
+#include "TPad.h"
+#include "TSystem.h"
+#include "TCut.h"
+#include "TChain.h"
+#include "TH2F.h"
+#include "TStatToolkit.h"
+#include "TROOT.h"
+#include "TDatabasePDG.h"
#include "TTreeStream.h"
#include <iostream>
#include <sstream>
#include "AliTPCLaserTrack.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
+#include "AliTPCParamSR.h"
+#include "TTimeStamp.h"
+#include "AliDCSSensorArray.h"
+#include "AliDCSSensor.h"
+#include "AliGRPObject.h"
+#include "AliTPCROC.h"
using namespace std;
ClassImp(AliTPCcalibLaser)
AliTPCcalibLaser::AliTPCcalibLaser():
- AliTPCcalibBase()
+ AliTPCcalibBase(),
+ fESD(0),
+ fESDfriend(0),
+ fNtracks(0),
+ fTracksMirror(336),
+ fTracksEsd(336),
+ fTracksEsdParam(336),
+ fTracksTPC(336),
+ fFullCalib(kTRUE),
+ fDeltaZ(336),
+ fDeltaP3(336),
+ fDeltaP4(336),
+ fDeltaPhi(336),
+ fDeltaPhiP(336),
+ fSignals(336),
+ //
+ fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
+ fHisNclIn(0), //->Number of clusters inner
+ fHisNclOut(0), //->Number of clusters outer
+ fHisNclIO(0), //->Number of cluster inner outer
+ fHisLclIn(0), //->Level arm inner
+ fHisLclOut(0), //->Level arm outer
+ fHisLclIO(0), //->Number of cluster inner outer
+ fHisdEdx(0), //->dEdx histo
+ fHisdZfit(0), //->distance to the mirror after linear fit
+ //
+ //
+ fHisChi2YIn1(0), //->chi2 y inner - line
+ fHisChi2YOut1(0), //->chi2 y inner - line
+ fHisChi2YIn2(0), //->chi2 y inner - parabola
+ fHisChi2YOut2(0), //->chi2 y inner - parabola
+ fHisChi2YIO1(0), //->chi2 y IO - common
+ fHisChi2ZIn1(0), //->chi2 z inner - line
+ fHisChi2ZOut1(0), //->chi2 z inner - line
+ fHisChi2ZIn2(0), //->chi2 z inner - parabola
+ fHisChi2ZOut2(0), //->chi2 z inner - parabola
+ fHisChi2ZIO1(0), //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
+ fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
+ fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
+ fHisPy2vP2In(0), //-> Curv P2inner - parabola
+ fHisPy2vP2Out(0), //-> Curv P2outer - parabola
+ fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
+ //
+ //
+ fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
+ fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
+ fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
+ fHisPz2vP2In(0), //-> Curv P2inner - parabola
+ fHisPz2vP2Out(0), //-> Curv P2outer - parabola
+ fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
+ //
+ fDeltaYres(336), //->2D histo of residuals
+ fDeltaZres(336), //->2D histo fo residuals
+ fDeltaYres2(336), //->2D histo of residuals
+ fDeltaZres2(336), //->2D histo fo residuals
+ fDeltaYresAbs(336), //->2D histo of residuals
+ fHisYAbsErrors(0), //-> total errors per beam in the abs res y analysis
+ fDeltaZresAbs(336), //->2D histo of residuals
+ fHisZAbsErrors(0), //-> total errors per beam in the abs res z analysis
+ //fDeltaYres3(336), //->2D histo of residuals
+ //fDeltaZres3(336), //->2D histo fo residuals
+ fFitAside(new TVectorD(5)),
+ fFitCside(new TVectorD(5)),
+ fFitACside(new TVectorD(6)),
+ fEdgeXcuts(3),
+ fEdgeYcuts(3),
+ fNClCuts(5),
+ fNcuts(0),
+ fBeamSectorOuter(336),
+ fBeamSectorInner(336),
+ fBeamOffsetYOuter(336),
+ fBeamSlopeYOuter(336),
+ fBeamOffsetYInner(336),
+ fBeamSlopeYInner(336),
+ fBeamOffsetZOuter(336),
+ fBeamSlopeZOuter(336),
+ fBeamOffsetZInner(336),
+ fBeamSlopeZInner(336),
+ fInverseSlopeZ(kTRUE),
+ fUseFixedDriftV(0),
+ fFixedFitAside0(0.0),
+ fFixedFitAside1(1.0),
+ fFixedFitCside0(0.0),
+ fFixedFitCside1(1.0)
{
//
// Constructor
//
+ fTracksEsdParam.SetOwner(kTRUE);
+ for (Int_t i=0; i<336; i++) {
+ fFitZ[i]=0;
+ fCounter[i]=0; //! counter of usage
+ fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+ fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
+ }
}
-AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title):
- AliTPCcalibBase()
+AliTPCcalibLaser::AliTPCcalibLaser(const Text_t *name, const Text_t *title, Bool_t full):
+ AliTPCcalibBase(),
+ fESD(0),
+ fESDfriend(0),
+ fNtracks(0),
+ fTracksMirror(336),
+ fTracksEsd(336),
+ fTracksEsdParam(336),
+ fTracksTPC(336),
+ fFullCalib(full),
+ //
+ fDeltaZ(336), // array of histograms of delta z for each track
+ fDeltaP3(336), // array of histograms of delta z for each track
+ fDeltaP4(336), // array of histograms of P3 for each track
+ fDeltaPhi(336), // array of histograms of P4 for each track
+ fDeltaPhiP(336), // array of histograms of delta z for each track
+ fSignals(336), // array of dedx signals
+ //
+ //
+ fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
+
+ fHisNclIn(0), //->Number of clusters inner
+ fHisNclOut(0), //->Number of clusters outer
+ fHisNclIO(0), //->Number of cluster inner outer
+ fHisLclIn(0), //->Level arm inner
+ fHisLclOut(0), //->Level arm outer
+ fHisLclIO(0), //->Number of cluster inner outer
+ fHisdEdx(0), //->dEdx histo
+ fHisdZfit(0), //->distance to the mirror after linear fit
+ //
+ //
+ fHisChi2YIn1(0), //->chi2 y inner - line
+ fHisChi2YOut1(0), //->chi2 y inner - line
+ fHisChi2YIn2(0), //->chi2 y inner - parabola
+ fHisChi2YOut2(0), //->chi2 y inner - parabola
+ fHisChi2YIO1(0), //->chi2 y IO - common
+ fHisChi2ZIn1(0), //->chi2 z inner - line
+ fHisChi2ZOut1(0), //->chi2 z inner - line
+ fHisChi2ZIn2(0), //->chi2 z inner - parabola
+ fHisChi2ZOut2(0), //->chi2 z inner - parabola
+ fHisChi2ZIO1(0), //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0(0), //-> delta y P0outer-P0inner - line
+ fHisPy2vP0(0), //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0(0), //-> delta y P0outer-P0inner - common parabola
+ fHisPy1vP1(0), //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1(0), //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1(0), //-> delta ky P1outer-P1inner - common parabola
+ fHisPy2vP2In(0), //-> Curv P2inner - parabola
+ fHisPy2vP2Out(0), //-> Curv P2outer - parabola
+ fHisPy3vP2IO(0), //-> Curv P2outerinner - common parabola
+ //
+ //
+ fHisPz1vP0(0), //-> delta z P0outer-P0inner - line
+ fHisPz2vP0(0), //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0(0), //-> delta z P0outer-P0inner - common parabola
+ fHisPz1vP1(0), //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1(0), //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1(0), //-> delta kz P1outer-P1inner - common parabola
+ fHisPz2vP2In(0), //-> Curv P2inner - parabola
+ fHisPz2vP2Out(0), //-> Curv P2outer - parabola
+ fHisPz3vP2IO(0), //-> Curv P2outerinner - common parabola
+ //
+ //
+ //
+ fDeltaYres(336),
+ fDeltaZres(336),
+ fDeltaYres2(336),
+ fDeltaZres2(336),
+ fDeltaYresAbs(336),
+ fHisYAbsErrors(0),
+ fDeltaZresAbs(336),
+ fHisZAbsErrors(0),
+ // fDeltaYres3(336),
+ //fDeltaZres3(336),
+ fFitAside(new TVectorD(5)), // drift fit - A side
+ fFitCside(new TVectorD(5)), // drift fit - C- side
+ fFitACside(new TVectorD(6)), // drift fit - AC- side
+ fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
+ fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
+ fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
+ fNcuts(0), // number of cuts
+ fBeamSectorOuter(336),
+ fBeamSectorInner(336),
+ fBeamOffsetYOuter(336),
+ fBeamSlopeYOuter(336),
+ fBeamOffsetYInner(336),
+ fBeamSlopeYInner(336),
+ fBeamOffsetZOuter(336),
+ fBeamSlopeZOuter(336),
+ fBeamOffsetZInner(336),
+ fBeamSlopeZInner(336),
+ fInverseSlopeZ(kTRUE),
+ fUseFixedDriftV(0),
+ fFixedFitAside0(0.0),
+ fFixedFitAside1(1.0),
+ fFixedFitCside0(0.0),
+ fFixedFitCside1(1.0)
{
SetName(name);
SetTitle(title);
//
// Constructor
//
-
+ fTracksEsdParam.SetOwner(kTRUE);
+ for (Int_t i=0; i<336; i++) {
+ fFitZ[i]=0;
+ fCounter[i]=0; //! counter of usage
+ fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+ fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
+ }
+}
+
+AliTPCcalibLaser::AliTPCcalibLaser(const AliTPCcalibLaser& calibLaser):
+ AliTPCcalibBase(calibLaser),
+ fESD(0),
+ fESDfriend(0),
+ fNtracks(0),
+ fTracksMirror(336),
+ fTracksEsd(336),
+ fTracksEsdParam(336),
+ fTracksTPC(336),
+ fFullCalib(calibLaser.fFullCalib),
+ //
+ fDeltaZ(calibLaser.fDeltaZ), // array of histograms of delta z for each track
+ fDeltaP3(((calibLaser.fDeltaP3))), // array of histograms of delta z for each track
+ fDeltaP4(((calibLaser.fDeltaP4))), // array of histograms of P3 for each track
+ fDeltaPhi(((calibLaser.fDeltaPhi))), // array of histograms of P4 for each track
+ fDeltaPhiP(((calibLaser.fDeltaPhiP))), // array of histograms of delta z for each track
+ fSignals(((calibLaser.fSignals))), // array of dedx signals
+ //
+ //
+ fHisLaser(0), // N dim histogram of laser
+ fHisLaserPad(0), // N dim histogram of laser
+ fHisLaserTime(0), // N dim histogram of laser
+
+ fHisNclIn(new TH2F(*(calibLaser.fHisNclIn))), //->Number of clusters inner
+ fHisNclOut(new TH2F(*(calibLaser.fHisNclOut))), //->Number of clusters outer
+ fHisNclIO(new TH2F(*(calibLaser.fHisNclIO))), //->Number of cluster inner outer
+ fHisLclIn(new TH2F(*(calibLaser.fHisLclIn))), //->Level arm inner
+ fHisLclOut(new TH2F(*(calibLaser.fHisLclOut))), //->Level arm outer
+ fHisLclIO(new TH2F(*(calibLaser.fHisLclIO))), //->Number of cluster inner outer
+ fHisdEdx(new TH2F(*(calibLaser.fHisdEdx))), //
+ fHisdZfit(new TH2F(*(calibLaser.fHisdZfit))), //->distance to the mirror after linear fit
+ //
+ //
+ fHisChi2YIn1(new TH2F(*(calibLaser.fHisChi2YIn1))), //->chi2 y inner - line
+ fHisChi2YOut1(new TH2F(*(calibLaser.fHisChi2YOut1))), //->chi2 y inner - line
+ fHisChi2YIn2(new TH2F(*(calibLaser.fHisChi2YIn2))), //->chi2 y inner - parabola
+ fHisChi2YOut2(new TH2F(*(calibLaser.fHisChi2YOut2))), //->chi2 y inner - parabola
+ fHisChi2YIO1(new TH2F(*(calibLaser.fHisChi2YIO1))), //->chi2 y IO - common
+ fHisChi2ZIn1(new TH2F(*(calibLaser.fHisChi2ZIn1))), //->chi2 z inner - line
+ fHisChi2ZOut1(new TH2F(*(calibLaser.fHisChi2ZOut1))), //->chi2 z inner - line
+ fHisChi2ZIn2(new TH2F(*(calibLaser.fHisChi2ZIn2))), //->chi2 z inner - parabola
+ fHisChi2ZOut2(new TH2F(*(calibLaser.fHisChi2ZOut2))), //->chi2 z inner - parabola
+ fHisChi2ZIO1(new TH2F(*(calibLaser.fHisChi2ZIO1))), //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0(new TH2F(*(calibLaser.fHisPy1vP0))), //-> delta y P0outer-P0inner - line
+ fHisPy2vP0(new TH2F(*(calibLaser.fHisPy2vP0))), //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0(new TH2F(*(calibLaser.fHisPy3vP0))), //-> delta y P0outer-P0inner - common parabola
+ fHisPy1vP1(new TH2F(*(calibLaser.fHisPy1vP1))), //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1(new TH2F(*(calibLaser.fHisPy2vP1))), //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1(new TH2F(*(calibLaser.fHisPy3vP1))), //-> delta ky P1outer-P1inner - common parabola
+ fHisPy2vP2In(new TH2F(*(calibLaser.fHisPy2vP2In))), //-> Curv P2inner - parabola
+ fHisPy2vP2Out(new TH2F(*(calibLaser.fHisPy2vP2Out))), //-> Curv P2outer - parabola
+ fHisPy3vP2IO(new TH2F(*(calibLaser.fHisPy3vP2IO))), //-> Curv P2outerinner - common parabola
+ //
+ //
+ fHisPz1vP0(new TH2F(*(calibLaser.fHisPz1vP0))), //-> delta z P0outer-P0inner - line
+ fHisPz2vP0(new TH2F(*(calibLaser.fHisPz2vP0))), //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0(new TH2F(*(calibLaser.fHisPz3vP0))), //-> delta z P0outer-P0inner - common parabola
+ fHisPz1vP1(new TH2F(*(calibLaser. fHisPz1vP1))), //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1(new TH2F(*(calibLaser.fHisPz2vP1))), //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1(new TH2F(*(calibLaser.fHisPz3vP1))), //-> delta kz P1outer-P1inner - common parabola
+ fHisPz2vP2In(new TH2F(*(calibLaser.fHisPz2vP2In))), //-> Curv P2inner - parabola
+ fHisPz2vP2Out(new TH2F(*(calibLaser. fHisPz2vP2Out))), //-> Curv P2outer - parabola
+ fHisPz3vP2IO(new TH2F(*(calibLaser.fHisPz3vP2IO))), //-> Curv P2outerinner - common parabola
+ //
+ //
+ fDeltaYres(((calibLaser.fDeltaYres))),
+ fDeltaZres(((calibLaser.fDeltaZres))),
+ fDeltaYres2(((calibLaser.fDeltaYres))),
+ fDeltaZres2(((calibLaser.fDeltaZres))),
+ fDeltaYresAbs(((calibLaser.fDeltaYresAbs))),
+ fHisYAbsErrors(new TH1F(*(calibLaser.fHisYAbsErrors))),
+ fDeltaZresAbs(((calibLaser.fDeltaZresAbs))),
+ fHisZAbsErrors(new TH1F(*(calibLaser.fHisZAbsErrors))),
+ // fDeltaYres3(((calibLaser.fDeltaYres))),
+ //fDeltaZres3(((calibLaser.fDeltaZres))),
+ fFitAside(new TVectorD(5)), // drift fit - A side
+ fFitCside(new TVectorD(5)), // drift fit - C- side
+ fFitACside(new TVectorD(6)), // drift fit - C- side
+ fEdgeXcuts(3), // cuts in local x direction; used in the refit of the laser tracks
+ fEdgeYcuts(3), // cuts in local y direction; used in the refit of the laser tracks
+ fNClCuts(5), // cuts on the number of clusters per tracklet; used in the refit of the laser tracks
+ fNcuts(0), // number of cuts
+ fBeamSectorOuter(336),
+ fBeamSectorInner(336),
+ fBeamOffsetYOuter(336),
+ fBeamSlopeYOuter(336),
+ fBeamOffsetYInner(336),
+ fBeamSlopeYInner(336),
+ fBeamOffsetZOuter(336),
+ fBeamSlopeZOuter(336),
+ fBeamOffsetZInner(336),
+ fBeamSlopeZInner(336),
+ fInverseSlopeZ(calibLaser.fInverseSlopeZ),
+ fUseFixedDriftV(calibLaser.fUseFixedDriftV),
+ fFixedFitAside0(calibLaser.fFixedFitAside0),
+ fFixedFitAside1(calibLaser.fFixedFitAside1),
+ fFixedFitCside0(calibLaser.fFixedFitCside0),
+ fFixedFitCside1(calibLaser.fFixedFitCside1)
+{
+ //
+ // copy constructor
+ //
+ for (Int_t i=0; i<336; i++) {
+ fFitZ[i]=0;
+ fCounter[i]=0; //! counter of usage
+ fClusterCounter[i]=0; //!couter of clusters in "sensitive are"
+ fClusterSatur[i]=0; //!couter of saturated clusters in "sensitive are"
+ }
+}
+
+
+
+AliTPCcalibLaser & AliTPCcalibLaser::operator=(const AliTPCcalibLaser& calibLaser){
+ //
+ // assgnment operator
+ //
+ if (this != &calibLaser) {
+ new (this) AliTPCcalibLaser(calibLaser);
+ }
+ return *this;
+
}
+
+
+
AliTPCcalibLaser::~AliTPCcalibLaser() {
//
// destructor
//
+ if ( fHisNclIn){
+ delete fHisLaser; //->
+ delete fHisLaserPad; //->
+ delete fHisLaserTime; //->
+
+ delete fHisNclIn; //->Number of clusters inner
+ delete fHisNclOut; //->Number of clusters outer
+ delete fHisNclIO; //->Number of cluster inner outer
+ delete fHisLclIn; //->Level arm inner
+ delete fHisLclOut; //->Level arm outer
+ delete fHisLclIO; //->Number of cluster inner outer
+ delete fHisdEdx;
+ delete fHisdZfit;
+ //
+ //
+ delete fHisChi2YIn1; //->chi2 y inner - line
+ delete fHisChi2YOut1; //->chi2 y inner - line
+ delete fHisChi2YIn2; //->chi2 y inner - parabola
+ delete fHisChi2YOut2; //->chi2 y inner - parabola
+ delete fHisChi2YIO1; //->chi2 y IO - common
+ delete fHisChi2ZIn1; //->chi2 z inner - line
+ delete fHisChi2ZOut1; //->chi2 z inner - line
+ delete fHisChi2ZIn2; //->chi2 z inner - parabola
+ delete fHisChi2ZOut2; //->chi2 z inner - parabola
+ delete fHisChi2ZIO1; //->chi2 z IO - common
+ //
+ //
+ delete fHisPy1vP0; //-> delta y P0outer-P0inner - line
+ delete fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
+ delete fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
+ delete fHisPy1vP1; //-> delta ky P1outer-P1inner - line
+ delete fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
+ delete fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
+ delete fHisPy2vP2In; //-> Curv P2inner - parabola
+ delete fHisPy2vP2Out; //-> Curv P2outer - parabola
+ delete fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
+ //
+ //
+ delete fHisPz1vP0; //-> delta z P0outer-P0inner - line
+ delete fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
+ delete fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
+ delete fHisPz1vP1; //-> delta kz P1outer-P1inner - line
+ delete fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
+ delete fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
+ delete fHisPz2vP2In; //-> Curv P2inner - parabola
+ delete fHisPz2vP2Out; //-> Curv P2outer - parabola
+ delete fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
+ }
+ //
+ //
+ //
+ fDeltaZ.SetOwner(); //-> array of histograms of delta z for each track
+ fDeltaP3.SetOwner(); //-> array of histograms of P3 for each track
+ fDeltaP4.SetOwner(); //-> array of histograms of P4 for each track
+ fDeltaPhi.SetOwner(); //-> array of histograms of delta z for each track
+ fDeltaPhiP.SetOwner(); //-> array of histograms of delta z for each track
+ fSignals.SetOwner(); //->Array of dedx signals
+
+ fDeltaZ.Delete(); //-> array of histograms of delta z for each track
+ fDeltaP3.Delete(); //-> array of histograms of P3 for each track
+ fDeltaP4.Delete(); //-> array of histograms of P4 for each track
+ fDeltaPhi.Delete(); //-> array of histograms of delta z for each track
+ fDeltaPhiP.Delete(); //-> array of histograms of delta z for each track
+ fSignals.Delete(); //->Array of dedx signals
+
+ fDeltaYres.SetOwner();
+ fDeltaYres.Delete();
+ delete fHisYAbsErrors;
+ fDeltaZres.SetOwner();
+ fDeltaZres.Delete();
+ delete fHisZAbsErrors;
+ fDeltaYres2.SetOwner();
+ fDeltaYres2.Delete();
+ fDeltaZres2.SetOwner();
+ fDeltaZres2.Delete();
+
+ fDeltaYresAbs.SetOwner();
+ fDeltaYresAbs.Delete();
+ fDeltaZresAbs.SetOwner();
+ fDeltaZresAbs.Delete();
+}
+
+
+
+void AliTPCcalibLaser::Process(AliESDEvent * event) {
+ //
+ //
+ // Loop over tracks and call Process function
+ //
+ const Int_t kMinTracks=20;
+ const Int_t kMinClusters=40;
+
+ fESD = event;
+ if (!fESD) {
+ return;
+ }
+ fESDfriend=static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
+ if (!fESDfriend) {
+ return;
+ }
+ if (fESDfriend->TestSkipBit()) return;
+ if (fESD->GetNumberOfTracks()<kMinTracks) return; //not enough tracks
+ AliDebug(4,Form("Event number in current file: %d",event->GetEventNumberInFile()));
+ //
+ // find CE background if present
+ //
+ if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks();
+ TH1D hisCE("hhisCE","hhisCE",100,-100,100);
+ for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
+ AliESDtrack *track=fESD->GetTrack(i);
+ if (!track) continue;
+ hisCE.Fill(track->GetZ());
+ hisCE.Fill(track->GetZ()+2);
+ hisCE.Fill(track->GetZ()-2);
+ }
+ //
+ //
+
+
+ fTracksTPC.Clear();
+ fTracksEsd.Clear();
+ fTracksEsdParam.Delete();
+ for (Int_t id=0; id<336;id++) {
+ fCounter[id]=0;
+ fClusterCounter[id]=0;
+ fClusterSatur[id]=0;
+ }
+ //
+ Int_t n=fESD->GetNumberOfTracks();
+ Int_t counter=0;
+ for (Int_t i=0;i<n;++i) {
+ AliESDfriendTrack *friendTrack=fESDfriend->GetTrack(i);
+ if (!friendTrack) continue;
+ AliESDtrack *track=fESD->GetTrack(i);
+ if (!track) continue;
+ Double_t binC = hisCE.GetBinContent(hisCE.FindBin(track->GetZ()));
+ if (binC>336) continue; //remove CE background
+ TObject *calibObject=0;
+ AliTPCseed *seed=0;
+ for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j)
+ if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
+ break;
+ if (track&&seed&&track->GetTPCNcls()>kMinClusters && seed->GetNumberOfClusters() >kMinClusters) {
+ //filter CE tracks
+ Int_t id = FindMirror(track,seed);
+ if (id>=0) counter++;
+ }
+ //
+ }
+ fNtracks=counter;
+ if (counter<kMinTracks) return;
+
+ //FitDriftV();
+ FitDriftV(0.2);
+ if (!fFullCalib) return;
+ static Bool_t init=kFALSE;
+ if (!init){
+ init = kTRUE; // way around for PROOF - to be investigated
+ MakeFitHistos();
+ }
+ //
+ for (Int_t id=0; id<336; id++){
+ //
+ if (!fTracksEsdParam.At(id)) continue;
+ if (fClusterSatur[id]>0.3) continue; // tracks in saturation
+ DumpLaser(id);
+ if ( AcceptLaser(id) && fFitZ[id]<0.5){
+ RefitLaserJW(id);
+ MakeDistHisto(id);
+ }
+ }
+
+}
+
+void AliTPCcalibLaser::MakeDistHisto(Int_t id){
+ //
+ //
+ //
+ // for (Int_t id=0; id<336; id++){
+ //
+ //
+ if (!fTracksEsdParam.At(id)) return;
+ if (!AcceptLaser(id)) return;
+ Double_t xhis[12]={0,0,0,0,0,0,0,0,0,0,0,0};
+ //
+ //
+ TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
+ if (!hisdz) MakeFitHistos();
+ hisdz = (TH1F*)fDeltaZ.At(id);
+ TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
+ TH1F * hisP4 = (TH1F*)fDeltaP4.At(id);
+ TH1F * hisdphi = (TH1F*)fDeltaPhi.At(id);
+ TH1F * hisdphiP = (TH1F*)fDeltaPhiP.At(id);
+ TH1F * hisSignal = (TH1F*)fSignals.At(id);
+ //
+
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ if (!param) return;
+ if (!ltrp) return;
+ if (!track) return;
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ //
+ Float_t dz = param->GetZ()-ltrp->GetZ();
+ Float_t dphi = (TMath::ATan2(xyz[1],xyz[0])- TMath::ATan2(lxyz[1],lxyz[0]))*254.;
+ Float_t dphiP = param->GetParameter()[2]-ltrp->GetParameter()[2];
+ if (hisdz) hisdz->Fill(dz);
+ if (hisP3) hisP3->Fill(param->GetParameter()[3]);
+ if (hisP4) hisP4->Fill(param->GetParameter()[4]);
+ if (hisdphi) hisdphi->Fill(dphi);
+ if (hisdphiP) hisdphiP->Fill(dphiP);
+ if (hisSignal) hisSignal->Fill(TMath::Sqrt(TMath::Abs(track->GetTPCsignal())));
+ // fill HisLaser
+ xhis[0] = ltrp->GetId();
+ xhis[1] = ltrp->GetSide();
+ xhis[2] = ltrp->GetRod();
+ xhis[3] = ltrp->GetBundle();
+ xhis[4] = ltrp->GetBeam();
+ xhis[5] = dphi;
+ xhis[6] = fFitZ[id];
+ xhis[7] = param->GetParameter()[2]-ltrp->GetParameter()[2]; //dp2
+ xhis[8] = param->GetParameter()[3]-ltrp->GetParameter()[3]; //dp3
+ xhis[9] = param->GetParameter()[4];
+ xhis[10]= track->GetTPCNcls();
+ xhis[11]= TMath::Sqrt(TMath::Abs(track->GetTPCsignal()));
+ // }
+ fHisLaser->Fill(xhis);
+ //
+
}
-void AliTPCcalibLaser::Process(AliESDtrack *track, Int_t run) {
+void AliTPCcalibLaser::FitDriftV(){
+ //
+ // Fit corrections to the drift velocity - linear approximation in the z and global y
+ //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
+ //
+ /*
+ Formulas:
+
+ z = s* (z0 - vd*(t-t0))
+
+ s - side -1 and +1
+ t0 - time 0
+ vd - nominal drift velocity
+ zs - miscalibrated position
+
+ zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
+ vr - relative change of the drift velocity
+ dzt - vd*dt
+ dr = zz0-s*z
+ ..
+ ==>
+ zs ~ z - s*vr*(z0-s*z)+s*dzt
+ --------------------------------
+ 1. Correction function vr constant:
+
+
+ dz = zs-z = -s*vr *(z0-s*z)+s*dzt
+ dzs/dl = dz/dl +s*s*vr*dz/dl
+ d(dz/dl) = vr*dz/dl
+ */
+
+ const Float_t kZCut = 240; // remove the closest laser beam
+ const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
+ const Float_t kDistCut = 3; // distance sigma cut
+ const Float_t kDistCutAbs = 0.25;
+ const Float_t kMinClusters = 60; // minimal amount of the clusters
+ const Float_t kMinSignal = 16; // minimal mean height of the signal
+ const Float_t kChi2Cut = 0.1; // chi2 cut to accept drift fit
+ static TLinearFitter fdriftA(3,"hyp2");
+ static TLinearFitter fdriftC(3,"hyp2");
+ static TLinearFitter fdriftAC(4,"hyp3");
+ TVectorD fitA(3),fitC(3),fitAC(4);
+
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCParam * tpcparam = calib->GetParameters();
+ //
+ for (Int_t id=0; id<336; id++) fFitZ[id]=0;
+
+ Float_t chi2A = 10;
+ Float_t chi2C = 10;
+ Float_t chi2AC = 10;
+ Int_t npointsA=0;
+ Int_t npointsC=0;
+ Int_t npointsAC=0;
+
+
+ for (Int_t iter=0; iter<3; iter++){
+ fdriftA.ClearPoints();
+ fdriftC.ClearPoints();
+ fdriftAC.ClearPoints();
+ //
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ if (!AcceptLaser(id)) continue;
+ if ( fClusterSatur[id]>kSaturCut) continue;
+ if ( fClusterCounter[id]<kMinClusters) continue;
+ AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ if (track->GetTPCsignal()<kMinSignal) continue;
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ if (TMath::Abs(lxyz[2])>kZCut) continue;
+ Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
+ if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
+ if (TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
+ if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
+
+ // drift distance
+ Double_t zlength = tpcparam->GetZLength(0);
+ Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+ Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+ Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ if (ltrp->GetSide()==0){
+ fdriftA.AddPoint(xxx,mdrift,1);
+ }else{
+ fdriftC.AddPoint(xxx,mdrift,1);
+ }
+ Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ fdriftAC.AddPoint(xxx2,mdrift,1);
+ }
+ //
+ if (fdriftA.GetNpoints()>10){
+ //
+ fdriftA.Eval();
+ if (iter==0) fdriftA.EvalRobust(0.7);
+ else fdriftA.EvalRobust(0.8);
+ fdriftA.GetParameters(fitA);
+ npointsA= fdriftA.GetNpoints();
+ chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
+ if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
+ if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
+ (*fFitAside)[0] = fitA[0];
+ (*fFitAside)[1] = fitA[1];
+ (*fFitAside)[2] = fitA[2];
+ (*fFitAside)[3] = fdriftA.GetNpoints();
+ (*fFitAside)[4] = chi2A;
+ }
+ }
+ if (fdriftC.GetNpoints()>10){
+ fdriftC.Eval();
+ if (iter==0) fdriftC.EvalRobust(0.7);
+ else fdriftC.EvalRobust(0.8);
+ //
+ fdriftC.GetParameters(fitC);
+ npointsC= fdriftC.GetNpoints();
+ chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
+ if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
+ if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
+ (*fFitCside)[0] = fitC[0];
+ (*fFitCside)[1] = fitC[1];
+ (*fFitCside)[2] = fitC[2];
+ (*fFitCside)[3] = fdriftC.GetNpoints();
+ (*fFitCside)[4] = chi2C;
+ }
+ }
+
+ if (fdriftAC.GetNpoints()>10&&fdriftC.GetNpoints()>10&&fdriftA.GetNpoints()>10){
+ fdriftAC.Eval();
+ if (iter==0) fdriftAC.EvalRobust(0.7);
+ else fdriftAC.EvalRobust(0.8);
+ //
+ fdriftAC.GetParameters(fitAC);
+ npointsAC= fdriftAC.GetNpoints();
+ chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
+ if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
+ }
+
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ //
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ Double_t zlength = tpcparam->GetZLength(0);
+ Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+ Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+
+ Float_t fz =0;
+ if (ltrp->GetSide()==0){
+ fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }else{
+ fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ if (npointsAC>10){
+ fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ fFitZ[id]=mdrift-fz;
+ }
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ TTimeStamp tstamp(fTime);
+ Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
+ Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
+ Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
+ Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
+ Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
+ Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
+ TVectorD vecGoofie(20);
+ AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
+ if (goofieArray)
+ for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
+ AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
+ if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
+ }
+
+ if (cstream){
+ (*cstream)<<"driftv"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ // Environment values
+ "press0="<<valuePressure0<<
+ "press1="<<valuePressure1<<
+ "pt0="<<ptrelative0<<
+ "pt1="<<ptrelative1<<
+ "temp0="<<temp0<<
+ "temp1="<<temp1<<
+ "vecGoofie.="<<&vecGoofie<<
+ //
+ //
+ "iter="<<iter<<
+ "driftA.="<<fFitAside<<
+ "driftC.="<<fFitCside<<
+ "driftAC.="<<fFitACside<<
+ "chi2A="<<chi2A<<
+ "chi2C="<<chi2C<<
+ "chi2AC="<<chi2AC<<
+ "nA="<<npointsA<<
+ "nC="<<npointsC<<
+ "nAC="<<npointsAC<<
+ "\n";
+ }
+ }
+ }
+}
+Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
//
+ // Fit corrections to the drift velocity - linear approximation in the z and global y
+ //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
//
+ // Source of outlyers :
+ // 0. Track in the saturation - postpeak
+ // 1. gating grid close the part of the signal for first bundle
+
+ // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
+ // 1. Robust fit is used in the itteration number 0
+ // only fraction of laser used
+ // 2. Only the tracks close to the fit used in the second itteration
+ /*
+ Formulas:
+
+ z = s* (z0 - vd*(t-t0))
+
+ s - side -1 and +1
+ t0 - time 0
+ vd - nominal drift velocity
+ zs - miscalibrated position
+
+ zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
+ vr - relative change of the drift velocity
+ dzt - vd*dt
+ dr = zz0-s*z
+ ..
+ ==>
+ zs ~ z - s*vr*(z0-s*z)+s*dzt
+ --------------------------------
+ 1. Correction function vr constant:
+
+
+ dz = zs-z = -s*vr *(z0-s*z)+s*dzt
+ dzs/dl = dz/dl +s*s*vr*dz/dl
+ d(dz/dl) = vr*dz/dl
+ */
+ const Int_t knLaser = 336; //n laser tracks
+ const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
+
+ const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
+ const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
+ const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
+ const Float_t kMinClusters = 40.; // minimal amount of the clusters
+ const Float_t kMinSignal = 2.5; // minimal mean height of the signal
+ const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
+ //
+ static TLinearFitter fdriftA(3,"hyp2");
+ static TLinearFitter fdriftC(3,"hyp2");
+ static TLinearFitter fdriftAC(4,"hyp3");
+ TVectorD fitA(3),fitC(3),fitAC(4);
+
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCParam * tpcparam = calib->GetParameters();
+ //
+ // reset old data
+ //
+ for (Int_t id=0; id<336; id++) fFitZ[id]=0;
+ if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
+ if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
+ for (Int_t i=0;i<5; i++){
+ (*fFitCside)[i]=0;
+ (*fFitAside)[i]=0;
+ }
+ //
+ //
+ Float_t chi2A = 10;
+ Float_t chi2C = 10;
+ Float_t chi2AC = 10;
+ Int_t npointsA=0;
+ Int_t npointsC=0;
+ Int_t npointsAC=0;
+ Int_t nbA[4]={0,0,0,0};
+ Int_t nbC[4]={0,0,0,0};
+ TVectorD vecZM(336); // measured z potion of laser
+ TVectorD vecA(336); // accepted laser
+ TVectorD vecZF(336); // fitted position
+ TVectorD vecDz(336); // deltaZ
+ TVectorD vecZS(336); // surveyed position of laser
+ // additional variable to cut
+ TVectorD vecdEdx(336); // dEdx
+ TVectorD vecSy(336); // shape y
+ TVectorD vecSz(336); // shape z
+ //
+ //
+ for (Int_t id=0; id<336; id++){
+ Int_t reject=0;
+ AliTPCLaserTrack *ltrp =
+ (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
+ vecZM(id)= (param==0) ? 0:param->GetZ();
+ vecZS(id)= ltrp->GetZ();
+ vecDz(id)= 0;
+ vecA(id)=1;
+ vecdEdx(id)=(seed)?seed->GetdEdx():0;
+ vecSy(id) =(seed)?seed->CookShape(1):0;
+ vecSz(id) =(seed)?seed->CookShape(2):0;
+ //
+ fFitZ[id]=0;
+ if (!fTracksEsdParam.At(id)) reject|=1;
+ if (!param) continue;
+ if (!AcceptLaser(id)) reject|=2;
+ if ( fClusterSatur[id]>kSaturCut) reject|=4;
+ if ( fClusterCounter[id]<kMinClusters) reject|=8;
+ vecA(id)=reject;
+ if (reject>0) continue;
+ if (ltrp->GetSide()==0){
+ npointsA++;
+ nbA[ltrp->GetBundle()]++;
+ }
+ if (ltrp->GetSide()==1){
+ npointsC++;
+ nbC[ltrp->GetBundle()]++;
+ }
+ }
+ //
+ // reject "bad events"
+ //
+ Bool_t isOK=kTRUE;
+ Int_t naA=0;
+ Int_t naC=0;
+ if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
+ isOK=kFALSE;
+ for (Int_t i=0;i<4;i++){
+ //count accepted for all layers
+ if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
+ if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
+ }
+ if (naA<3 &&naC<3) isOK=kFALSE;
+ if (isOK==kFALSE) return kFALSE;
+
+ //
+ //
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fdriftA.ClearPoints();
+ fdriftC.ClearPoints();
+ fdriftAC.ClearPoints();
+ npointsA=0; npointsC=0; npointsAC=0;
+ //
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ if (!AcceptLaser(id)) continue;
+ if ( fClusterSatur[id]>kSaturCut) continue;
+ if ( fClusterCounter[id]<kMinClusters) continue;
+ AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ if (track->GetTPCsignal()<kMinSignal) continue;
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
+ //if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
+ if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
+ if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
+
+// // drift distance
+// Double_t zlength = tpcparam->GetZLength(0);
+// Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+// Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+ //
+ Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+ Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+ Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
+
+ Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ if (iter==0 &<rp->GetBundle()==0) continue;
+ // skip bundle 0 - can be wrong in the 0 iteration
+ if (ltrp->GetSide()==0){
+ fdriftA.AddPoint(xxx,mdrift,1);
+ }else{
+ fdriftC.AddPoint(xxx,mdrift,1);
+ }
+ Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ fdriftAC.AddPoint(xxx2,mdrift,1);
+ }
+ //
+ if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
+ //
+ fdriftA.Eval();
+ //if (iter==0) fdriftA.FixParameter(2,0); //fix global y gradient
+ npointsA= fdriftA.GetNpoints();
+ chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
+ fdriftA.EvalRobust(kFraction[iter]);
+ fdriftA.GetParameters(fitA);
+ if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
+ if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
+ (*fFitAside)[0] = fitA[0];
+ (*fFitAside)[1] = fitA[1];
+ (*fFitAside)[2] = fitA[2];
+ (*fFitAside)[3] = fdriftA.GetNpoints();
+ (*fFitAside)[4] = chi2A;
+ }
+ }
+ if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
+ fdriftC.Eval();
+ //if (iter==0) fdriftC.FixParameter(2,0); //fix global y gradient
+ npointsC= fdriftC.GetNpoints();
+ chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
+ fdriftC.EvalRobust(kFraction[iter]);
+ fdriftC.GetParameters(fitC);
+ if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
+ if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
+ (*fFitCside)[0] = fitC[0];
+ (*fFitCside)[1] = fitC[1];
+ (*fFitCside)[2] = fitC[2];
+ (*fFitCside)[3] = fdriftC.GetNpoints();
+ (*fFitCside)[4] = chi2C;
+ }
+ }
+
+ if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
+ fdriftAC.Eval();
+ //if (iter==0) fdriftAC.FixParameter(2,0); //fix global y gradient
+ npointsAC= fdriftAC.GetNpoints();
+ chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
+ fdriftAC.EvalRobust(kFraction[iter]);
+ fdriftAC.GetParameters(fitAC);
+ if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
+ (*fFitACside)[0] = fitAC[0];
+ (*fFitACside)[1] = fitAC[1];
+ (*fFitACside)[2] = fitAC[2];
+ (*fFitACside)[3] = fdriftAC.GetNpoints();
+ (*fFitACside)[4] = chi2AC;
+ }
+
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ //
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ //Double_t zlength = tpcparam->GetZLength(0);
+ //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+ //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+ Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+ Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+ Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
+
+ Float_t fz =0;
+ if (ltrp->GetSide()==0){
+ fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }else{
+ fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ if (npointsAC>10){
+ //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ fFitZ[id]=mdrift-fz;
+ vecZF[id]=fz;
+ vecDz[id]=mdrift-fz;
+ }
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ TTimeStamp tstamp(fTime);
+ Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
+ Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
+ Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
+ Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
+ Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
+ Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
+ TVectorD vecGoofie(20);
+ AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
+ if (goofieArray)
+ for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
+ AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
+ if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
+ }
+
+ if (cstream){
+ (*cstream)<<"driftvN"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ // Environment values
+ "press0="<<valuePressure0<<
+ "press1="<<valuePressure1<<
+ "pt0="<<ptrelative0<<
+ "pt1="<<ptrelative1<<
+ "temp0="<<temp0<<
+ "temp1="<<temp1<<
+ "vecGoofie.="<<&vecGoofie<<
+ //
+ //
+ "vecZM.="<<&vecZM<< // measured z position
+ "vecZS.="<<&vecZS<< // surveyed z position
+ "vecZF.="<<&vecZF<< // fitted z position
+ "vecDz.="<<&vecDz<< // fitted z position
+ "vecA.="<<&vecA<< // accept laser flag
+ "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
+ "vecSy.="<<&vecSy<< // shape y - to cut on
+ "vecSz.="<<&vecSz<< // shape z - to cut on
+ //
+ "iter="<<iter<<
+ "driftA.="<<fFitAside<<
+ "driftC.="<<fFitCside<<
+ "driftAC.="<<fFitACside<<
+ "chi2A="<<chi2A<<
+ "chi2C="<<chi2C<<
+ "chi2AC="<<chi2AC<<
+ "nA="<<npointsA<<
+ "nC="<<npointsC<<
+ "nAC="<<npointsAC<<
+ "\n";
+ /*
+ //
+ variables to check in debug mode:
+ //
+ chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
+ chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
+ chainDriftN->SetAlias("driftF","vecZF.fElements");
+ chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
+ TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
+ TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
+
+ */
+ }
+ }
+ }
+ return kTRUE;
+}
+
+Float_t AliTPCcalibLaser::GetDistance(AliExternalTrackParam *param, AliTPCLaserTrack *ltrp){
+ //
+ // get distance between mirror and laser track
+ //
+ //
+ Double_t xyz[3];
+ Double_t lxyz[3];
+ param->GetXYZ(xyz);
+ ltrp->GetXYZ(lxyz);
+ //
+ //
+ Double_t dist = 0;
+ //radial distance
+ dist+=TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX());
+ //
+ // z distance
+ // apply drift correction if already exist
+ //
+ Float_t dz = 0;
+ if (ltrp->GetSide()==0){
+ if ((*fFitAside)[1]>0.) dz = ((*fFitAside)[0]+(*fFitAside)[1]*lxyz[2]+(*fFitAside)[2]*lxyz[1])-xyz[2];
+ }else{
+ if ((*fFitCside)[1]>0.) dz = ((*fFitCside)[0]+(*fFitCside)[1]*lxyz[2]+(*fFitCside)[2]*lxyz[1])-xyz[2];
+ }
+ if (TMath::Abs(dz)>3*(TMath::Abs(lxyz[2]-xyz[2])+1)) dz= TMath::Abs(lxyz[2]-xyz[2]);
+ dist+=TMath::Abs(dz);
+ //
+ // phi dist - divergence on 50 cm
+ //
+ dist = TMath::Abs((param->GetParameter()[2]-ltrp->GetParameter()[2])*50);
+ return dist;
+}
+
+
+Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){
+ //
+ //
+ //
+ /*
+ TCut cutP0("cutP0","abs((atan2(x1,x0)-atan2(lx1,lx0))*254)<1.5");
+ TCut cutP1("cutP1","abs(LTr.fP[1]-Tr.fP[1])<30");
+ TCut cutP2("cutP2","abs(LTr.fP[2]-Tr.fP[2])<0.03");
+ TCut cutP3("cutP3","abs(Tr.fP[3])<0.05");
+
+ TCut cutA = cutP0+cutP1+cutP2+cutP3+cutP4;
+ */
+ AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ //AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ Double_t xyz[3];
+ Double_t lxyz[3];
+ param->GetXYZ(xyz);
+ ltrp->GetXYZ(lxyz);
+ if (TMath::Abs((TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*254)>1.5) return kFALSE; //cut y- P0
+ if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>30) return kFALSE; // cutZ -P1
+ if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.03) return kFALSE; // cut -P2
+ if (TMath::Abs(param->GetParameter()[3])>0.05) return kFALSE; // cut Tl -P3
+ //
+ //
+
+ return kTRUE;
+}
+
+Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){
//
- // 1. Propagate track to the mirror radius
- Float_t kRadius = 271.6;
- if (!track->GetOuterParam()) return;
+ // Find corresponding mirror
+ // add the corresponding tracks
+
+ if (!track->GetOuterParam()) return -1;
+
+ Float_t kRadius0 = 252;
+ Float_t kRadius = 254.2;
+ Int_t countercl=0;
+ Float_t counterSatur=0;
+ Int_t csideA =0;
+ Int_t csideC =0;
+ for (Int_t irow=158;irow>-1;--irow) {
+ AliTPCclusterMI *c=seed->GetClusterPointer(irow);
+ if (!c) continue;
+ Double_t pedgeY = c->GetX()*TMath::DegToRad()*(10)-TMath::Abs(c->GetY());
+ Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
+ if (pedgeY<3) continue;
+ if (pedgeX<3) continue;
+ countercl++;
+ if (c->GetDetector()%36<18) csideA++;
+ if (c->GetDetector()%36>=18) csideC++;
+ if (c->GetMax()>900) counterSatur++;
+ }
+ counterSatur/=(countercl+1);
+ //
+ //
+ //
+ if (csideA<0.9*seed->GetNumberOfClusters() && csideC<0.9*seed->GetNumberOfClusters()) return 0; // cross laser track can not happen
+
+ Int_t side= 0;
+ if (csideC>0.5*seed->GetNumberOfClusters()) side=1;
+
+
AliExternalTrackParam param(*(track->GetOuterParam()));
- AliTracker::PropagateTrackTo(¶m,270.,0.0005,3,kTRUE);
- AliTracker::PropagateTrackTo(¶m,kRadius,0.0005,0.1,kTRUE);
+ AliTracker::PropagateTrackTo(¶m,kRadius0,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),3,kTRUE);
+ AliTracker::PropagateTrackTo(¶m,kRadius,TDatabasePDG::Instance()->GetParticle("mu+")->Mass(),0.1,kTRUE);
AliTPCLaserTrack ltr;
AliTPCLaserTrack *ltrp=0x0;
- Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m);
- if (id!=-1) ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
- else ltrp=<r;
+ // AliTPCLaserTrack *ltrpjw=0x0;
+ //
+ Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m,side);
+ // Int_t idjw = AliTPCLaserTrack::IdentifyTrackJW(¶m);
+ //AliDebug(4,Form("Identified Track: %03d (%03d)",id,idjw));
+
+ if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id)))
+ ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ else
+ ltrp=<r;
+
+ if (id<0) return -1;
+ if (ltrp->GetSide()!=side) return -1;
+ fCounter[id]++;
+ //
+ //
+ //
+ //
+ if (counterSatur>fClusterSatur[id]) fClusterSatur[id]=counterSatur;
+ //
//
+ Float_t radius=TMath::Abs(ltrp->GetX());
+ param.Rotate(ltrp->GetAlpha());
+ AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kFALSE);
+ //
+ if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id);
+ Bool_t accept=kTRUE;
+ //
+ // choose closer track
+ //
+ AliExternalTrackParam * param0 = (AliExternalTrackParam *)fTracksEsdParam.At(id);
+ if (param0){
+ Float_t dist0=GetDistance(param0,ltrp);
+ Float_t dist1=GetDistance(¶m,ltrp);
+ if (dist0<dist1) accept=kFALSE;
+ }
+
+ if (accept){
+ fClusterCounter[id]=countercl;
+ fTracksEsdParam.AddAt(param.Clone(),id);
+ fTracksEsd.AddAt(track,id);
+ fTracksTPC.AddAt(seed,id);
+ }
+ return id;
+}
+
+
+
+void AliTPCcalibLaser::DumpLaser(Int_t id) {
+ //
+ // Dump Laser info to the tree
+ //
+ AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ //
+ // Fast laser ID
+ //
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ Float_t dist3D = GetDistance(param,ltrp);
+ Float_t dist0 = (TMath::ATan2(xyz[1],xyz[0])-TMath::ATan2(lxyz[1],lxyz[0]))*param->GetX();
+ Float_t distphi = (param->GetParameter()[2]-ltrp->GetParameter()[2])*50;
+
+
if (fStreamLevel>0){
TTreeSRedirector *cstream = GetDebugStreamer();
- Double_t xyz[3];
- Double_t pxyz[3];
- param.GetXYZ(xyz);
- param.GetPxPyPz(pxyz);
- Int_t side = (param.GetZ()>0) ? 0:1;
- Int_t beam = 0;
- if (side==0) beam = TMath::Nint((180*atan2(xyz[1],xyz[0])/TMath::Pi()+20)/60.);
- if (side==1) beam = TMath::Nint((180*atan2(xyz[1],xyz[0])/TMath::Pi()-20)/60.);
- //Int_t id(180*atan2(x1,x0)/pi+20)/60.;
- Int_t bundle=TMath::Nint(param.GetZ()/80.);
+ Int_t time = fESD->GetTimeStamp();
+ Bool_t accept = AcceptLaser(id);
if (cstream){
(*cstream)<<"Track"<<
- "run="<<run<<
+ //
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+
"id="<<id<<
- "fSide="<<side<< // side A-C
- "fBeam="<<beam<< // phi id
- "fBundle="<<bundle<< // laser Z
+ "accept="<<accept<<
+ "driftA.="<<fFitAside<<
+ "driftC.="<<fFitCside<<
+ "time="<<time<<
+ "dist3D="<<dist3D<<
+ "dist0="<<dist0<<
+ "distphi="<<distphi<<
+ //
+ //
+ "counter="<<fCounter[id]<<
+ "clcounter="<<fClusterCounter[id]<<
+ "satur="<<fClusterSatur[id]<<
+ "fz="<<fFitZ[id]<<
//
"LTr.="<<ltrp<<
"Esd.="<<track<<
- "Tr.="<<¶m<<
+ "Tr.="<<param<<
"x0="<<xyz[0]<<
"x1="<<xyz[1]<<
"x2="<<xyz[2]<<
"px0="<<pxyz[0]<<
"px1="<<pxyz[1]<<
"px2="<<pxyz[2]<<
+ //
+ "lx0="<<lxyz[0]<<
+ "lx1="<<lxyz[1]<<
+ "lx2="<<lxyz[2]<<
+ "lpx0="<<lpxyz[0]<<
+ "lpx1="<<lpxyz[1]<<
+ "lpx2="<<lpxyz[2]<<
"\n";
}
}
}
-void AliTPCcalibLaser::Analyze(){
+void AliTPCcalibLaser::RefitLaserJW(Int_t id){
//
+ // Refit the track with different tracklet models:
+ // 1. Per ROC using the kalman filter, different edge cuts
+ // 2. Per ROC linear in y and z
+ // 3. Per ROC quadratic in y and z
+ // 4. Per track offset for each sector, linear for each sector, common quadratic
+ // store x, y, z information for all models and the cluster to calculate the residuals
//
+
+ // The clusters which do not fulfill given criteria are skipped
//
-}
+ // Cluters removed from fit
+ // 1. Extended shape > kShapeCut
+ // 2. In saturation Max > kMax
+ // 3. Distance to edge < cutEdge
+ //
+ // Clusters not used for the calibration:
+ //
+ // 1. Extended shape > kShapeCut
+ // 2. In saturation Max > kMax
+
+
+ AliTPCseed *track = (AliTPCseed*)fTracksTPC.At(id);
+ AliExternalTrackParam *extparam=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
+
+ AliTPCclusterMI dummyCl;
+
+ //two tracklets
+ Int_t kMaxTracklets=2;
+ Float_t kShapeCut = 1.3;
+ Float_t kRatioCut = 2.;
+
+ Float_t kMax = 900.;
+
+
+ //=============================================//
+ // Linear Fitters for the Different Approaches //
+ //=============================================//
+ //linear fit model in y and z; inner - outer sector, combined with offset
+ static TLinearFitter fy1I(2,"hyp1");
+ static TLinearFitter fy1O(2,"hyp1");
+ static TLinearFitter fz1I(2,"hyp1");
+ static TLinearFitter fz1O(2,"hyp1");
+ static TLinearFitter fy1IO(3,"hyp2");
+ static TLinearFitter fz1IO(3,"hyp2");
+ //quadratic fit model in y and z; inner - sector
+ static TLinearFitter fy2I(3,"hyp2");
+ static TLinearFitter fy2O(3,"hyp2");
+ static TLinearFitter fz2I(3,"hyp2");
+ static TLinearFitter fz2O(3,"hyp2");
+ //common quadratic fit for IROC and OROC in y and z
+ static TLinearFitter fy4(5,"hyp4");
+ static TLinearFitter fz4(5,"hyp4");
+
+
+ //set standard cuts
+ if ( fNcuts==0 ){
+ fNcuts=1;
+ fEdgeXcuts[0]=4;
+ fEdgeYcuts[0]=3;
+ fNClCuts[0]=20;
+ }
+ //=============================//
+ // Loop over all Tracklet Cuts //
+ //=============================//
+ for (Int_t icut=0; icut<fNcuts; icut++){
+ Float_t xinMin = 2500, xinMax=-90;
+ Float_t xoutMin=2500, xoutMax=-90;
+ Float_t msigmaYIn=0;
+ Float_t msigmaYOut=0;
+ Float_t mqratioIn=0;
+ Float_t mqratioOut=0;
+
+
+ AliDebug(4,Form("Processing cut %d for track with ID %d",icut,id));
+ //cut parameters
+ Double_t edgeCutX = fEdgeXcuts[icut];
+ Double_t edgeCutY = fEdgeYcuts[icut];
+ Int_t nclCut = (Int_t)fNClCuts[icut];
+ //=========================//
+ // Parameters to calculate //
+ //=========================//
+ //fit parameter inner, outer tracklet and combined fit
+ TVectorD vecy1resInner(2),vecz1resInner(2); //pol1 fit parameters inner
+ TVectorD vecy2resInner(3),vecz2resInner(3); //pol2 fit parameters inner
+ //
+ TVectorD vecy1resOuter(2),vecz1resOuter(2); //pol1 fit parameters outer
+ TVectorD vecy2resOuter(3),vecz2resOuter(3); //pol2 fit parameters outer
+ TVectorD vecy4res(5),vecz4res(5);
+ TVectorD vecy1resIO(3),vecz1resIO(3);
+ // cluster and track positions for each row - used for residuals
+ TVectorD vecgX(159); // global X
+ TVectorD vecgY(159); // global Y
+ TVectorD vecgZ(159); // global Z
+
+ TVectorD vecX(159); // x is the same for all (row center)
+ TVectorD vecYkalman(159); // y from kalman fit
+ TVectorD vecZkalman(159); // z from kalman fit
+ TVectorD vecY1(159); // y from pol1 fit per ROC
+ TVectorD vecZ1(159); // z from pol1 fit per ROC
+ TVectorD vecY1IO(159); // y from pol1 fit per ROC
+ TVectorD vecZ1IO(159); // z from pol1 fit per ROC
+ TVectorD vecY2(159); // y from pol2 fit per ROC
+ TVectorD vecZ2(159); // z from pol2 fit per ROC
+ TVectorD vecY4(159); // y from sector fit
+ TVectorD vecZ4(159); // z from sector fit
+ TVectorD vecClY(159); // y cluster position
+ TVectorD vecClZ(159); // z cluster position
+ TVectorD vecSec(159); // sector for each row
+ TVectorD isReject(159); // flag - cluster to be rejected
+ //chi2 of fits
+ Double_t chi2I1z=0; // chi2 of pol1 fit in z (inner)
+ Double_t chi2I1y=0; // chi2 of pol1 fit in y (inner)
+ Double_t chi2O1z=0; // chi2 of pol1 fit in z (outer)
+ Double_t chi2O1y=0; // chi2 of pol1 fit in y (outer)
+ Double_t chi2IO1z=0; // chi2 of pol1 fit in z (outer)
+ Double_t chi2IO1y=0; // chi2 of pol1 fit in y (outer)
+ Double_t chi2I2z=0; // chi2 of pol2 fit in z (inner)
+ Double_t chi2I2y=0; // chi2 of pol2 fit in y (inner)
+ Double_t chi2O2z=0; // chi2 of pol2 fit in z (outer)
+ Double_t chi2O2y=0; // chi2 of pol2 fit in y (outer)
+ Double_t chi2IOz=0; // chi2 of hyp4 fit in z (inner+outer)
+ Double_t chi2IOy=0; // chi2 of hyp4 fit in y (inner+outer)
+ //more
+ Int_t innerSector = -1; // number of inner sector
+ Int_t outerSector = -1; // number of outer sector
+ Int_t nclI=0; // number of clusters (inner)
+ Int_t nclO=0; // number of clusters (outer)
+ //=================================================//
+ // Perform the Kalman Fit using the Tracklet Class //
+ //=================================================//
+ AliTPCTracklet::SetEdgeCut(edgeCutX,edgeCutY);
+ TObjArray tracklets=
+ AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,
+ kFALSE,nclCut,kMaxTracklets);
+ // tracklet pointers
+ AliTPCTracklet *trInner = (AliTPCTracklet*)tracklets.At(0);
+ AliTPCTracklet *trOuter = (AliTPCTracklet*)tracklets.At(1);
+ AliTPCTracklet *tr=0x0;
+ AliTPCTracklet dummy;
+ //continue if we didn't find a tracklet
+ if ( !trInner && !trOuter ) continue;
+ //================================================================================//
+ // Swap Inner and Outer if necessary (inner sector is definde by smaller local x) //
+ //================================================================================//
+ if ( trInner && trOuter ){
+ if ( !trInner->GetInner() || !trOuter->GetInner() ) continue;
+ if ( trInner->GetInner()->GetX() > trOuter->GetInner()->GetX() ){
+ tr = trInner;
+ trInner=trOuter;
+ trOuter=tr;
+ AliDebug(5,Form("Swapped Sectors: %02d (%f) <-> %02d (%f)", trOuter->GetSector(), trOuter->GetInner()->GetX(), trInner->GetSector(), trInner->GetInner()->GetX()));
+ }
+ } else {
+ if ( trInner ){
+ if ( !trInner->GetInner() ) continue;
+ trOuter=&dummy;
+ if ( trInner->GetSector()>35 ){
+ trOuter=trInner;
+ trInner=&dummy;
+ }
+ } else { //trOuter
+ if ( !trOuter->GetInner() ) continue;
+ trInner=&dummy;
+ if ( trOuter->GetSector()<36 ){
+ trInner=trOuter;
+ trOuter=&dummy;
+ }
+ }
+ }
+ innerSector = trInner->GetSector();
+ if ( innerSector>=0 ) AliDebug(5,Form("Found inner Sector %02d at X %.2f", innerSector, trInner->GetInner()->GetX()));
+ outerSector = trOuter->GetSector();
+ if ( outerSector>=0 ) AliDebug(5,Form("Found outer Sector %02d at X %.2f", outerSector, trOuter->GetInner()->GetX()));
+
+ // array of clusters
+ TClonesArray arrCl("AliTPCclusterMI",159);
+ arrCl.ExpandCreateFast(159);
+ //=======================================//
+ // fill fitters with cluster information //
+ //=======================================//
+ AliDebug(3,"Filing Cluster Information");
+
+ //
+ //
+
+ for (Int_t irow=158;irow>-1;--irow) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ AliTPCclusterMI &cl = (AliTPCclusterMI&) (*arrCl[irow]);
+ cl=dummyCl;
+ vecX[irow] = 0;
+ vecClY[irow] = 0;
+ vecClZ[irow] = 0;
+ Float_t meanY=0, sumY=0;
+ for (Int_t drow=-1;drow<=1;drow++) {
+ if (irow+drow<0) continue;
+ if (irow+drow>158) continue;
+ AliTPCclusterMI *ccurrent=track->GetClusterPointer(irow);
+ if (!ccurrent) continue;
+ Int_t roc = static_cast<Int_t>(ccurrent->GetDetector());
+ if ( roc!=innerSector && roc!=outerSector ) continue;
+ if (ccurrent->GetX()<10) continue;
+ if (ccurrent->GetY()==0) continue;
+ meanY+=ccurrent->GetY();
+ sumY++;
+ }
+ if (sumY>0) meanY/=sumY;
+
+ //
+ vecSec[irow]=-1;
+ if (!c) continue;
+ Double_t pedgeY = c->GetX()*TMath::Tan(TMath::DegToRad()*(10))-TMath::Abs(meanY);
+ Double_t pedgeX = TMath::Min((irow)*0.75, (159.-irow)*1.5);
+
+ //
+ Int_t roc = static_cast<Int_t>(c->GetDetector());
+ if ( roc!=innerSector && roc!=outerSector ) continue;
+ vecSec[irow]=roc;
+ //store clusters in clones array
+ cl=*c;
+ //
+ if (c->GetMax()<4) continue; // noise cluster?
+ if (TMath::Abs(c->GetY())<0.0001) continue; // noise cluster?
+ //cluster position
+ vecX[irow] = c->GetX();
+ vecClY[irow] = c->GetY();
+ vecClZ[irow] = c->GetZ();
+ //
+// Float_t gxyz[3];
+// c->GetGlobalXYZ(gxyz);
+// vecgX[irow] = gxyz[0];
+// vecgY[irow] = gxyz[1];
+// vecgZ[irow] = gxyz[2];
+ //
+ Double_t x = vecX[irow]-133.4; //reference is between IROC and OROC
+ Double_t y = vecClY[irow];
+ Double_t z = vecClZ[irow];
+ //
+ Double_t x2[2]={x,x*x}; //linear and parabolic parameters
+ Double_t x4[4]={0,0,0,0}; //hyp4 parameters
+ Double_t xIO[2]={0,x}; //common linear + offset IROC-OROC
+ if ( roc == innerSector ) {
+ x4[0]=1; //offset inner - outer sector
+ x4[1]=x; //slope parameter inner sector
+ xIO[0]=1;
+ } else {
+ x4[2]=x; //slope parameter outer sector
+ }
+ x4[3]=x*x; //common parabolic shape
+ if (pedgeX < fEdgeXcuts[icut]) continue;
+ if (pedgeY < fEdgeYcuts[icut]) continue;
+ if (c->GetMax()>900) continue; // cluster in saturation
+ //
+ if ( roc==innerSector ){
+ fy1I.AddPoint(x2,y);
+ fz1I.AddPoint(x2,z);
+ fy2I.AddPoint(x2,y);
+ fz2I.AddPoint(x2,z);
+ ++nclI;
+ if (c->GetX()<xinMin) xinMin=c->GetX();
+ if (c->GetX()>xinMax) xinMax=c->GetX();
+
+ msigmaYIn +=TMath::Sqrt(c->GetSigmaY2());
+ mqratioIn +=c->GetMax()/c->GetQ();
+
+ }
+ if ( roc==outerSector ){
+ fy1O.AddPoint(x2,y);
+ fz1O.AddPoint(x2,z);
+ fy2O.AddPoint(x2,y);
+ fz2O.AddPoint(x2,z);
+ ++nclO;
+ if (c->GetX()<xoutMin) xoutMin=c->GetX();
+ if (c->GetX()>xoutMax) xoutMax=c->GetX();
+ msigmaYOut +=TMath::Sqrt(c->GetSigmaY2());
+ mqratioOut +=c->GetMax()/c->GetQ();
+
+ }
+ fy4.AddPoint(x4,y);
+ fz4.AddPoint(x4,z);
+ fy1IO.AddPoint(xIO,y);
+ fz1IO.AddPoint(xIO,z);
+ }
+ if (nclI>0) {
+ msigmaYIn/=nclI;
+ mqratioIn/=nclI;
+ }
+ if (nclO>0) {
+ msigmaYOut/=nclO;
+ mqratioOut/=nclO;
+ }
+ //======================================//
+ // Evaluate and retrieve fit parameters //
+ //======================================//
+ AliDebug(5,Form("Evaluating fits with inner (outer) Sec: %02d (%02d)",innerSector,outerSector));
+ //inner sector
+ if ( (innerSector>-1) && (fy1I.GetNpoints()>0) ){
+ fy1I.Eval();
+ fz1I.Eval();
+ fy2I.Eval();
+ fz2I.Eval();
+ fy1I.GetParameters(vecy1resInner);
+ fz1I.GetParameters(vecz1resInner);
+ fy2I.GetParameters(vecy2resInner);
+ fz2I.GetParameters(vecz2resInner);
+ chi2I1y=fy1I.GetChisquare()/(fy1I.GetNpoints()-2);
+ chi2I1z=fz1I.GetChisquare()/(fz1I.GetNpoints()-2);
+ chi2I2y=fy2I.GetChisquare()/(fy2I.GetNpoints()-3);
+ chi2I2z=fz2I.GetChisquare()/(fz2I.GetNpoints()-3);
+ }
+ //outer sector
+ if ( (outerSector>-1) && (fy1O.GetNpoints()>0) ){
+ fy1O.Eval();
+ fz1O.Eval();
+ fy2O.Eval();
+ fz2O.Eval();
+ fy1O.GetParameters(vecy1resOuter);
+ fz1O.GetParameters(vecz1resOuter);
+ fy2O.GetParameters(vecy2resOuter);
+ fz2O.GetParameters(vecz2resOuter);
+ chi2O1y=fy1O.GetChisquare()/(fy1O.GetNpoints()-2);
+ chi2O1z=fz1O.GetChisquare()/(fz1O.GetNpoints()-2);
+ chi2O2y=fy2O.GetChisquare()/(fy2O.GetNpoints()-3);
+ chi2O2z=fz2O.GetChisquare()/(fz2O.GetNpoints()-3);
+ }
+ //combined hyp4 fit
+ if ( innerSector>0 && outerSector>0 ){
+ if (fy4.GetNpoints()>0) {
+ fy4.Eval();
+ fy4.GetParameters(vecy4res);
+ chi2IOy=fy4.GetChisquare()/(fy4.GetNpoints()-5);
+ }
+ if (fz4.GetNpoints()>0) {
+ fz4.Eval();
+ fz4.GetParameters(vecz4res);
+ chi2IOz=fz4.GetChisquare()/(fz4.GetNpoints()-5);
+ }
+ if (fy1IO.GetNpoints()>0) {
+ fy1IO.Eval();
+ fy1IO.GetParameters(vecy1resIO);
+ chi2IO1y=fy1IO.GetChisquare()/(fy1IO.GetNpoints()-3);
+ }
+ if (fz1IO.GetNpoints()>0) {
+ fz1IO.Eval();
+ fz1IO.GetParameters(vecz1resIO);
+ chi2IO1z=fz1IO.GetChisquare()/(fz1IO.GetNpoints()-3);
+ }
+ }
+ //clear points
+ fy4.ClearPoints(); fz4.ClearPoints();
+ fy1I.ClearPoints(); fy1O.ClearPoints();
+ fz1I.ClearPoints(); fz1O.ClearPoints();
+ fy2I.ClearPoints(); fy2O.ClearPoints();
+ fz2I.ClearPoints(); fz2O.ClearPoints();
+ fy1IO.ClearPoints(); fz1IO.ClearPoints();
+ //==============================//
+ // calculate tracklet positions //
+ //==============================//
+ AliDebug(4,"Calculate tracklet positions");
+ for (Int_t irow=158;irow>-1;--irow) {
+ isReject[irow]=0;
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) { // no cluster in given sectors
+ isReject[irow]+=1;
+ }
+
+ if (!c) { //no cluster
+ isReject[irow]+=2;
+ }else{
+ if (c->GetMax()>kMax){ //saturation
+ isReject[irow]+=4;
+ }
+ if ( vecSec[irow] == outerSector ) { //extended shape
+ if (c->GetMax()/c->GetQ()> mqratioOut*kRatioCut) isReject[irow]+=8;
+ if (TMath::Sqrt(c->GetSigmaY2())>msigmaYOut*kShapeCut) isReject[irow]+=16;
+ }else{
+ if (c->GetMax()/c->GetQ()> mqratioIn*kRatioCut) isReject[irow]+=8;
+ if (TMath::Sqrt(c->GetSigmaY2())>msigmaYIn*kShapeCut) isReject[irow]+=16;
+ }
+ }
+ if ( vecSec[irow]==-1 ) continue; //no cluster info
+ if ( vecSec[irow]!=innerSector && vecSec[irow]!=outerSector ) continue;
+ tr=&dummy;
+ Double_t x=vecX[irow];
+ Double_t xref=x-133.4;
+ //
+ Double_t yoffInner=0;
+ Double_t zoffInner=0;
+ Double_t yoffInner1=0;
+ Double_t zoffInner1=0;
+ Double_t yslopeInner=0;
+ Double_t yslopeOuter=0;
+ Double_t zslopeInner=0;
+ Double_t zslopeOuter=0;
+ //positions of hyperplane fits
+ if ( vecSec[irow] == outerSector ) {
+ tr=trOuter;
+ vecY1[irow]=vecy1resOuter[0]+vecy1resOuter[1]*xref;
+ vecZ1[irow]=vecz1resOuter[0]+vecz1resOuter[1]*xref;
+ vecY2[irow]=vecy2resOuter[0]+vecy2resOuter[1]*xref+vecy2resOuter[2]*xref*xref;
+ vecZ2[irow]=vecz2resOuter[0]+vecz2resOuter[1]*xref+vecz2resOuter[2]*xref*xref;
+ yslopeOuter=vecy4res[3];
+ zslopeOuter=vecz4res[3];
+ } else {
+ tr=trInner;
+ vecY1[irow]=vecy1resInner[0]+vecy1resInner[1]*xref;
+ vecZ1[irow]=vecz1resInner[0]+vecz1resInner[1]*xref;
+ vecY2[irow]=vecy2resInner[0]+vecy2resInner[1]*xref+vecy2resInner[2]*xref*xref;
+ vecZ2[irow]=vecz2resInner[0]+vecz2resInner[1]*xref+vecz2resInner[2]*xref*xref;
+ yoffInner=vecy4res[1];
+ zoffInner=vecz4res[1];
+ yoffInner1=vecy1resIO[1];
+ zoffInner1=vecz1resIO[1];
+ yslopeInner=vecy4res[2];
+ zslopeInner=vecz4res[2];
+ }
+ vecY1IO[irow]=vecy1resIO[0]+yoffInner1+vecy1resIO[2]*xref;
+ vecZ1IO[irow]=vecz1resIO[0]+zoffInner1+vecz1resIO[2]*xref;
+ vecY4[irow]=vecy4res[0]+yoffInner+yslopeInner*xref+yslopeOuter*xref+vecy4res[4]*xref*xref;
+ vecZ4[irow]=vecz4res[0]+zoffInner+zslopeInner*xref+zslopeOuter*xref+vecz4res[4]*xref*xref;
+ //positions of kalman fits
+ Double_t gxyz[3],xyz[3];
+ AliExternalTrackParam *param = 0x0;
+ //
+ param=tr->GetInner();
+ if (param){
+ param->GetXYZ(gxyz);
+ Float_t bz = AliTracker::GetBz(gxyz);
+ param->GetYAt(x, bz, xyz[1]);
+ param->GetZAt(x, bz, xyz[2]);
+ vecYkalman[irow]=xyz[1];
+ vecZkalman[irow]=xyz[2];
+ }
+ //
+ //
+ //
+
+ }
+ //=====================================================================//
+ // write results from the different tracklet fits with debug streamers //
+ //=====================================================================//
+ if (fStreamLevel>4){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (cstream){
+ Float_t dedx = track->GetdEdx();
+ (*cstream)<<"FitModels"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ //
+ "cutNr=" << icut <<
+ "edgeCutX=" << edgeCutX <<
+ "edgeCutY=" << edgeCutY <<
+ "nclCut=" << nclCut <<
+ "innerSector="<< innerSector <<
+ "outerSector="<< outerSector <<
+ "dEdx=" << dedx <<
+ "LTr.=" << ltrp <<
+ "Tr.=" << extparam <<
+ "yPol1In.=" << &vecy1resInner <<
+ "zPol1In.=" << &vecz1resInner <<
+ "yPol1InOut.="<< &vecy1resIO <<
+ "zPol1InOut.="<< &vecz1resIO <<
+ "yPol2In.=" << &vecy2resInner <<
+ "zPol2In.=" << &vecz2resInner <<
+ "yPol1Out.=" << &vecy1resOuter <<
+ "zPol1Out.=" << &vecz1resOuter <<
+ "yPol2Out.=" << &vecy2resOuter <<
+ "zPol2Out.=" << &vecz2resOuter <<
+ "yInOut.=" << &vecy4res <<
+ "zInOut.=" << &vecz4res <<
+ "chi2y1In=" << chi2I1y <<
+ "chi2z1In=" << chi2I1z <<
+ "chi2y1InOut="<< chi2IO1y <<
+ "chi2z1InOut="<< chi2IO1z <<
+ "chi2y1Out=" << chi2O1y <<
+ "chi2z1Out=" << chi2O1z <<
+ "chi2y2In=" << chi2I2y <<
+ "chi2z2In=" << chi2I2z <<
+ "chi2y2Out=" << chi2O2y <<
+ "chi2z2Out=" << chi2O2z <<
+ "chi2yInOut=" << chi2IOy <<
+ "chi2zInOut=" << chi2IOz <<
+ "trletIn.=" << trInner <<
+ "trletOut.=" << trOuter <<
+ "nclI=" << nclI <<
+ "nclO=" << nclO <<
+ "xinMin=" << xinMin<<
+ "xinMax=" << xinMax<<
+ "xoutMin=" << xoutMin<<
+ "xoutMax=" << xoutMax<<
+ "msigmaYIn=" <<msigmaYIn<<
+ "msigmaYOut=" <<msigmaYOut<<
+ "mqratioIn=" <<mqratioIn<<
+ "mqratioOut=" << mqratioOut <<
+ "\n";
+ }
+ }
+ // wirte residuals information
+ if (fStreamLevel>5){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (cstream){
+ Float_t dedx = track->GetdEdx();
+ (*cstream)<<"Residuals"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ //
+ "cutNr=" << icut <<
+ "edgeCutX=" << edgeCutX <<
+ "edgeCutY=" << edgeCutY <<
+ "nclCut=" << nclCut <<
+ "LTr.=" << ltrp <<
+ "Tr.=" << extparam<<
+ "dEdx=" << dedx <<
+ "Cl.=" << &arrCl <<
+ "vX.=" << &vecgX<< // global x
+ "vY.=" << &vecgY<< // global y
+ "vZ.=" << &vecgZ<< // global z
+ "TrX.=" << &vecX <<
+ "TrYpol1.=" << &vecY1 <<
+ "TrZpol1.=" << &vecZ1 <<
+ "TrYpol2.=" << &vecY2 <<
+ "TrZpol2.=" << &vecZ2 <<
+ "TrYpol1InOut.="<< &vecY1IO <<
+ "TrZpol1InOut.="<< &vecZ1IO <<
+ "TrYInOut.=" << &vecY4 <<
+ "TrZInOut.=" << &vecZ4 <<
+ "ClY.=" << &vecClY <<
+ "ClZ.=" << &vecClZ <<
+ "isReject.=" << &isReject<<
+ "sec.=" << &vecSec <<
+ "nclI=" << nclI <<
+ "nclO=" << nclO <<
+ "xinMin=" << xinMin<<
+ "xinMax=" << xinMax<<
+ "xoutMin=" << xoutMin<<
+ "xoutMax=" << xoutMax<<
+ "msigmaYIn=" <<msigmaYIn<<
+ "msigmaYOut=" <<msigmaYOut<<
+ "mqratioIn=" <<mqratioIn<<
+ "mqratioOut=" << mqratioOut <<
+ "yInOut.=" << &vecy4res <<
+ "zInOut.=" << &vecz4res <<
+ //chi2s
+ "chi2y1In=" << chi2I1y << //
+ "chi2z1In=" << chi2I1z <<
+ "chi2y1Out=" << chi2O1y <<
+ "chi2z1Out=" << chi2O1z <<
+ "chi2y1InOut="<< chi2IO1y <<
+ "chi2z1InOut="<< chi2IO1z <<
+ "chi2y2In=" << chi2I2y <<
+ "chi2z2In=" << chi2I2z <<
+ "chi2y2Out=" << chi2O2y <<
+ "chi2z2Out=" << chi2O2z <<
+ "chi2yInOut=" << chi2IOy <<
+ "chi2zInOut=" << chi2IOz <<
+ // fit parameters
+ "yPol1In.=" << &vecy1resInner <<
+ "zPol1In.=" << &vecz1resInner <<
+ "yPol2In.=" << &vecy2resInner <<
+ "zPol2In.=" << &vecz2resInner <<
+ "yPol1Out.=" << &vecy1resOuter <<
+ "zPol1Out.=" << &vecz1resOuter <<
+ "yPol1InOut.="<< &vecy1resIO <<
+ "zPol1InOut.="<< &vecz1resIO <<
+ "yPol2Out.=" << &vecy2resOuter <<
+ "zPol2Out.=" << &vecz2resOuter <<
+
+ "\n";
+
+ }
+ }
+ //==========================//
+ // Fill Residual Histograms //
+ //==========================//
+ if (!fHisNclIn) MakeFitHistos();
+
+ TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
+ TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
+ TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
+ TH2F *profz2 = (TH2F*)fDeltaZres2.UncheckedAt(id);
+ TH2F *profyabs = (TH2F*)fDeltaYresAbs.UncheckedAt(id);
+ TH2F *profzabs = (TH2F*)fDeltaZresAbs.UncheckedAt(id);
+ // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
+ //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
+ //
+ for (Int_t irow=158;irow>-1;--irow) {
+ if (vecSec[irow]==-1)continue; // no cluster info
+ if (isReject[irow]>0.5) continue; //
+ Double_t ycl = vecClY[irow];
+ Double_t yfit = vecY1[irow];
+ Double_t yfit2 = vecY2[irow];
+ Double_t x = vecX[irow];
+ Double_t yabsbeam = -1000;
+ if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
+ yabsbeam = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id];
+ else if(innerSector==fBeamSectorInner[id])
+ yabsbeam = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id];
+
+ // Double_t yfit3 = vecY2[irow];
+ Double_t zcl = vecClZ[irow];
+ Double_t zfit = vecZ1[irow];
+ Double_t zfit2 = vecZ2[irow];
+ //Double_t zfit3 = vecZ2[irow];
+
+ // dz abs
+ // The expressions for zcorrected has been obtained by
+ // inverting the fits in the FitDriftV() method (ignoring the
+ // global y dependence for now):
+ // A side:
+ // 250 - zmeasured = [0] + [1]*(250 - zreal) + .... (yglobal)
+ // =>
+ // zreal = (zmeasured + [0] - (1-[1])*250.0)/[1]
+ //
+ // C side:
+ // 250 + zmeasured = [0] + [1]*(250+zreal) + .... (yglobal)
+ // =>
+ // zreal = (zmeasured - [0] + (1 - [1])*250.0)/[1]
+
+ Double_t dzabs = -1000;
+ Double_t zcorrected = -1000;
+ if (ltrp->GetSide()==0){
+ if ((*fFitAside)[1]>0. || fUseFixedDriftV) {
+ // ignore global y dependence for now
+ zcorrected = 0;
+ if(!fUseFixedDriftV)
+ zcorrected = (zcl + (*fFitAside)[0] -
+ (1.0-(*fFitAside)[1])*250.0)/(*fFitAside)[1];
+ else
+ zcorrected = (zcl + fFixedFitAside0 -
+ (1.0-fFixedFitAside1)*250.0)/fFixedFitAside1;
+ // zcorrected = zcl;
+ if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
+ dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
+ else if(innerSector==fBeamSectorInner[id])
+ dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
+ }
+ } else {
+ if ((*fFitCside)[1]>0. || fUseFixedDriftV) {
+
+ if(!fUseFixedDriftV)
+ zcorrected = (zcl - (*fFitCside)[0] +
+ (1.0-(*fFitCside)[1])*250.0)/(*fFitCside)[1];
+ else
+ zcorrected = (zcl - fFixedFitCside0 +
+ (1.0-fFixedFitCside1)*250.0)/fFixedFitCside1;
+
+ // zcorrected = zcl;
+ if(vecSec[irow]==outerSector && outerSector==fBeamSectorOuter[id])
+ dzabs = zcorrected -fBeamSlopeZOuter[id]*x -fBeamOffsetZOuter[id];
+ else if(innerSector==fBeamSectorInner[id])
+ dzabs = zcorrected -fBeamSlopeZInner[id]*x -fBeamOffsetZInner[id];
+ }
+ }
+
+ if (TMath::Abs(yfit-ycl)<2&&TMath::Abs(zfit-zcl)<2){
+ if (profy){
+ profy->Fill(irow,ycl-yfit);
+ profy2->Fill(irow,ycl-yfit2);
+ if(yabsbeam<-100) {
+ fHisYAbsErrors->Fill(id);
+ // profyabs->Fill(irow,-0.99);
+ } else
+ profyabs->Fill(irow,ycl-yabsbeam);
+
+ // profy3->Fill(irow,ycl-yfit3);
+ }
+ if (profz) {
+ profz->Fill(irow,zcl-zfit);
+ profz2->Fill(irow,zcl-zfit2);
+ //profz3->Fill(irow,zcl-zfit3);
+ if(dzabs<-100) {
+
+ fHisZAbsErrors->Fill(id);
+ }else
+ profzabs->Fill(irow,dzabs);
+ }
+ }
+ }
+ //
+ //
+ // Fill laser fit histograms
+ //
+ Float_t dedx = track->GetdEdx();
+ if (nclI>20&&nclO>20){
+ fHisNclIn->Fill(id,nclI); //->Number of clusters inner
+ fHisNclOut->Fill(id,nclO); //->Number of clusters outer
+ fHisNclIO->Fill(id,nclI+nclO); //->Number of cluster inner outer
+ //
+ fHisLclIn->Fill(id,xinMax-xinMin); //->Level arm inner
+ fHisLclOut->Fill(id,xoutMax-xoutMin); //->Level arm outer
+ fHisLclIO->Fill(id,xoutMax-xinMin); //->Number of cluster inner outer
+ //
+ fHisdEdx->Fill(id,TMath::Sqrt(TMath::Abs(dedx)));
+ fHisdZfit->Fill(id,fFitZ[id]);
+ //
+ //
+ fHisChi2YIn1->Fill(id,TMath::Sqrt(chi2I1y)); //->chi2 y inner - line
+ fHisChi2YOut1->Fill(id,TMath::Sqrt(chi2O1y)); //->chi2 y inner - line
+ fHisChi2YIn2->Fill(id,TMath::Sqrt(chi2I2y)); //->chi2 y inner - parabola
+ fHisChi2YOut2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 y inner - parabola
+ fHisChi2YIO1->Fill(id,TMath::Sqrt(chi2IOy)); //->chi2 y IO - common
+
+
+ fHisChi2ZIn1->Fill(id,TMath::Sqrt(chi2I1z)); //->chi2 z inner - line
+ fHisChi2ZOut1->Fill(id,TMath::Sqrt(chi2O1z)); //->chi2 z inner - line
+ fHisChi2ZIn2->Fill(id,TMath::Sqrt(chi2O2y)); //->chi2 z inner - parabola
+ fHisChi2ZOut2->Fill(id,TMath::Sqrt(chi2O2z)); //->chi2 z inner - parabola
+ fHisChi2ZIO1->Fill(id,TMath::Sqrt(chi2IOz)); //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0->Fill(id,vecy1resOuter[0]-vecy1resInner[0]); //-> delta y P0outer-P0inner - line
+ fHisPy2vP0->Fill(id,vecy2resOuter[0]-vecy2resInner[0]); //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0->Fill(id,vecy4res[1]); //-> delta y P0outer-P0inner - common parabola
+ //
+ fHisPy1vP1->Fill(id,vecy1resOuter[1]-vecy1resInner[1]); //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1->Fill(id,vecy2resOuter[1]-vecy2resInner[1]); //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1->Fill(id,vecy4res[3]-vecy4res[2]); //-> delta ky P1outer-P1inner - common parabola
+ //
+ fHisPy3vP2IO->Fill(id,vecy4res[4]); //-> Curv P2outerinner - common parabola
+ fHisPz1vP0->Fill(id,vecz1resOuter[0]-vecz1resInner[0]); //-> delta z P0outer-P0inner - line
+ fHisPz2vP0->Fill(id,vecz2resOuter[0]-vecz2resInner[0]); //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0->Fill(id,vecz4res[1]); //-> delta z P0outer-P0inner - common parabola
+ //
+ fHisPz1vP1->Fill(id,vecz1resOuter[1]-vecz1resInner[1]); //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1->Fill(id,vecz2resOuter[1]-vecz2resInner[1]); //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1->Fill(id,vecz4res[3]-vecz4res[2]); //-> delta kz P1outer-P1inner - common parabola
+ fHisPz3vP2IO->Fill(id,vecz4res[4]); //-> Curv P2outerinner - common parabola
+ }
+ if(nclI>20){
+ fHisPy2vP2In->Fill(id,vecy2resInner[2]); //-> Curv P2inner - parabola
+ fHisPz2vP2In->Fill(id,vecz2resInner[2]); //-> Curv P2inner - parabola
+ }
+ //
+ if (nclO>20){
+ fHisPz2vP2Out->Fill(id,vecz2resOuter[2]); //-> Curv P2outer - parabola
+ fHisPy2vP2Out->Fill(id,vecy2resOuter[2]); //-> Curv P2outer - parabola
+ }
+
+ }
+ //
+ // Fill raw THnSparses
+ //
+ for (Int_t irow=0;irow<159;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetMax()>800) continue; // saturation cut
+ //if (TMath::Sqrt(TMath::Abs(c->GetSigmaY2()))>1) continue; // saturation cut
+ //
+ Double_t deltaY=c->GetY()-(*ltrp->GetVecLY())[irow];
+ Double_t deltaZ=c->GetZ()-(*ltrp->GetVecLZ())[irow];
+ //TString axisName[6]={"Delta","bin", "rms shape", "Q", "row","trackID"}
+ Double_t xyz[6]={0, 0, 0,TMath::Sqrt(c->GetMax()),irow,id};
+ xyz[0]=deltaY;
+ xyz[1]=c->GetPad();
+ xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaY2()));
+ fHisLaserPad->Fill(xyz);
+ xyz[0]=deltaZ;
+ xyz[1]=c->GetTimeBin();
+ xyz[2]=TMath::Sqrt(TMath::Abs(c->GetSigmaZ2()));
+ fHisLaserTime->Fill(xyz);
+ }
+}
+
+
+
+void AliTPCcalibLaser::DumpMeanInfo(Int_t run){
+ //
+ // Dump information about laser beams
+ // isOK variable indicates usability of the beam
+ // Beam is not usable if:
+ // a. No entries in range (krmsCut0)
+ // b. Big sperad (krmscut1)
+ // c. RMSto fit sigma bigger then (kmultiCut)
+ // d. Too big angular spread
+ //
+
+ const Float_t krmsCut0=0.001;
+ const Float_t krmsCut1=0.16;
+ const Float_t kmultiCut=2;
+ const Float_t kcutP0=0.002;
+ AliMagF* magF= dynamic_cast<AliMagF*> (TGeoGlobalMagField::Instance()->GetField());
+ Double_t xyz[3]={90,0,10}; // tmp. global position
+ Double_t bxyz[3]={90,0,10}; // tmp. mag field integral - cylindrical
+ Double_t bgxyz[3]={90,0,10}; // tmp. mag field integral - local
+ //
+ AliTPCcalibLaser *laser = this;
+ TTreeSRedirector *pcstream = new TTreeSRedirector("laserMean.root");
+ TF1 fg("fg","gaus");
+ AliTPCParam * tpcparam = 0;
+ // start set up for absolute residuals analysis
+ //
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ tpcparam = calib->GetParameters();
+ if (!tpcparam) tpcparam = new AliTPCParamSR;
+ tpcparam->Update();
+ AliGRPObject *grp = AliTPCcalibDB::GetGRP(run);
+ Float_t current=0;
+ Float_t bfield = 0, bz=0;
+
+ if (grp){
+ Float_t polarity = (grp->GetL3Polarity()>0) ? -1.:1;
+ current = grp->GetL3Current((AliGRPObject::Stats)0);
+ bfield = polarity*5*current/30000.;
+ bz = polarity*5*current/30000.;
+ printf("Run%d\tL3 current%f\tBz\t%f\n",run,current,bz);
+ }
+
+ SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
+ SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
+ TLinearFitter lfabsyInner(2);
+ lfabsyInner.SetFormula("1 ++ x");
+ TLinearFitter lfabszInner(2);
+ lfabszInner.SetFormula("1 ++ x");
+
+ TLinearFitter lfabsyOuter(2);
+ lfabsyOuter.SetFormula("1 ++ x");
+ TLinearFitter lfabszOuter(2);
+ lfabszOuter.SetFormula("1 ++ x");
+ // end set up for absolute residuals analysis
+
+ //
+ //
+ for (Int_t id=0; id<336; id++){
+ Bool_t isOK=kTRUE;
+ TH1F * hisphi = (TH1F*)laser->fDeltaPhi.At(id);
+ TH1F * hisphiP = (TH1F*)laser->fDeltaPhiP.At(id);
+ TH1F * hisZ = (TH1F*)laser->fDeltaZ.At(id);
+ TH1F * hisP3 = (TH1F*)laser->fDeltaP3.At(id);
+ TH1F * hisP4 = (TH1F*)laser->fDeltaP4.At(id);
+ TH1F * hisS = (TH1F*)laser->fSignals.At(id);
+ //if (!hisphi) continue;
+ Double_t entries = (hisphi==0)? 0: hisphi->GetEntries();
+ //if (entries<minEntries) continue;
+ //
+ AliTPCLaserTrack *ltrp = (AliTPCLaserTrack*)fTracksMirror.At(id);
+ if (!ltrp) {
+ AliTPCLaserTrack::LoadTracks();
+ ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ }
+ ltrp->UpdatePoints();
+ pcstream->GetFile()->cd();
+ if (hisphi) hisphi->Write();
+ if (hisphiP) hisphiP->Write();
+ if (hisZ) hisZ->Write();
+ if (hisP3) hisP3->Write();
+ if (hisP4) hisP4->Write();
+
+ Float_t meanphi = hisphi->GetMean();
+ Float_t rmsphi = hisphi->GetRMS();
+ //
+ Float_t meanphiP = hisphiP->GetMean();
+ Float_t rmsphiP = hisphiP->GetRMS();
+ Float_t meanZ = hisZ->GetMean();
+ Float_t rmsZ = hisZ->GetRMS();
+ if (hisphi->GetRMS()>hisphi->GetBinWidth(1))
+ hisphi->Fit(&fg,"","",hisphi->GetMean()-4*hisphi->GetRMS(),hisphi->GetMean()+4*hisphi->GetRMS());
+ Double_t gphi1 = fg.GetParameter(1);
+ Double_t gphi2 = fg.GetParameter(2);
+ if (hisphiP->GetRMS()>0)
+ hisphiP->Fit(&fg,"","",hisphiP->GetMean()-4*hisphiP->GetRMS(),hisphiP->GetMean()+4*hisphiP->GetRMS());
+ Double_t gphiP1 = fg.GetParameter(1);
+ Double_t gphiP2 = fg.GetParameter(2);
+ //
+ if (hisZ->GetRMS()>hisZ->GetBinWidth(1))
+ hisZ->Fit(&fg,"","");
+ Double_t gz1 = fg.GetParameter(1);
+ Double_t gz2 = fg.GetParameter(2);
+ //
+ if (hisP3->GetRMS()>hisP3->GetBinWidth(1))
+ hisP3->Fit(&fg,"","",hisP3->GetMean()-4*hisP3->GetRMS(),hisP3->GetMean()+4*hisP3->GetRMS());
+ Double_t gp31 = fg.GetParameter(1);
+ Double_t gp32 = fg.GetParameter(2);
+ Double_t meanp3 = hisP3->GetMean();
+ Double_t rmsp3 = hisP3->GetRMS();
+ //
+ if (hisP4->GetRMS()>hisP4->GetBinWidth(1))
+ hisP4->Fit(&fg,"","",hisP4->GetMean()-4*hisP4->GetRMS(),hisP4->GetMean()+4*hisP4->GetRMS());
+ Double_t gp41 = fg.GetParameter(1);
+ Double_t gp42 = fg.GetParameter(2);
+ Double_t meanp4 = hisP4->GetMean();
+ Double_t rmsp4 = hisP4->GetRMS();
+ //
+ Float_t meanS=hisS->GetMean();
+ //
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+
+ if (rmsphi<krmsCut0) isOK=kFALSE; // empty in range - not entries inside
+ if (rmsphi>krmsCut1) isOK=kFALSE; // empty in range - not entries inside
+ if (rmsphi>krmsCut0) if (gphi2/rmsphi>kmultiCut) isOK=kFALSE; // multi peak structure
+ if (gphiP2>kcutP0) isOK=kFALSE;
+ //
+ //
+ //
+ TH1 *his =0;
+ //
+ his = fHisNclIn->ProjectionY("aaa",id+1,id+1);
+ Float_t mnclIn = his->GetMean();
+ delete his;
+ his = fHisNclOut->ProjectionY("aaa",id+1,id+1);
+ Float_t mnclOut = his->GetMean();
+ delete his;
+ his = fHisNclIO->ProjectionY("aaa",id+1,id+1);
+ Float_t mnclIO = his->GetMean();
+ delete his;
+ his = fHisLclIn->ProjectionY("aaa",id+1,id+1);
+ Float_t mLclIn = his->GetMean();
+ delete his;
+ his = fHisLclOut->ProjectionY("aaa",id+1,id+1);
+ Float_t mLclOut = his->GetMean();
+ delete his;
+ his = fHisLclIO->ProjectionY("aaa",id+1,id+1);
+ Float_t mLclIO = his->GetMean();
+ delete his;
+ //
+ his = fHisdEdx->ProjectionY("aaa",id+1,id+1);
+ Float_t mdEdx = his->GetMean();
+ delete his;
+ //
+ //
+ //
+ //
+ his = fHisChi2YIn1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
+ Float_t mChi2YIn1= his->GetMean();
+ delete his;
+ his = fHisChi2YOut1->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - line
+ Float_t mChi2YOut1 = his->GetMean();
+ delete his;
+ his = fHisChi2YIn2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
+ Float_t mChi2YIn2 = his->GetMean();
+ delete his;
+ his = fHisChi2YOut2->ProjectionY("aaa",id+1,id+1); //->chi2 y inner - parabola
+ Float_t mChi2YOut2 = his->GetMean();
+ delete his;
+ his = fHisChi2YIO1->ProjectionY("aaa",id+1,id+1); //->chi2 y IO - common
+ Float_t mChi2YIO1 = his->GetMean();
+ delete his;
+ his = fHisChi2ZIn1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
+ Float_t mChi2ZIn1 = his->GetMean();
+ delete his;
+ his = fHisChi2ZOut1->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - line
+ Float_t mChi2ZOut1 = his->GetMean();
+ delete his;
+ his = fHisChi2ZIn2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
+ Float_t mChi2ZIn2 = his->GetMean();
+ delete his;
+ his = fHisChi2ZOut2->ProjectionY("aaa",id+1,id+1); //->chi2 z inner - parabola
+ Float_t mChi2ZOut2 = his->GetMean();
+ delete his;
+ his = fHisChi2ZIO1->ProjectionY("aaa",id+1,id+1); //->chi2 z IO - common
+ Float_t mChi2ZIO1 = his->GetMean();
+ delete his;
+ //
+ // fit res. histos
+ //
+ his = fHisdZfit->ProjectionY("aaa",id+1,id+1);
+ Float_t edZfit = his->GetEntries();
+ Float_t mdZfit = his->GetMean();
+ Float_t rdZfit = his->GetRMS();
+ delete his;
+
+ his = fHisPy1vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - line
+ Float_t ePy1vP0 = his->GetEntries();
+ Float_t mPy1vP0 = his->GetMean();
+ Float_t rPy1vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPy2vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - parabola
+ Float_t ePy2vP0 = his->GetEntries();
+ Float_t mPy2vP0 = his->GetMean();
+ Float_t rPy2vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPy3vP0->ProjectionY("aaa",id+1,id+1); //-> delta y P0outer-P0inner - common parabola
+ Float_t ePy3vP0 = his->GetEntries();
+ Float_t mPy3vP0 = his->GetMean();
+ Float_t rPy3vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPy1vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - line
+ Float_t ePy1vP1 = his->GetEntries();
+ Float_t mPy1vP1 = his->GetMean();
+ Float_t rPy1vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPy2vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - parabola
+ Float_t ePy2vP1 = his->GetEntries();
+ Float_t mPy2vP1 = his->GetMean();
+ Float_t rPy2vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPy3vP1->ProjectionY("aaa",id+1,id+1); //-> delta ky P1outer-P1inner - common parabola
+ Float_t ePy3vP1 = his->GetEntries();
+ Float_t mPy3vP1 = his->GetMean();
+ Float_t rPy3vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPy2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
+ Float_t ePy2vP2In = his->GetEntries();
+ Float_t mPy2vP2In = his->GetMean();
+ Float_t rPy2vP2In = his->GetRMS();
+ delete his;
+
+ his = fHisPy2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
+ Float_t ePy2vP2Out = his->GetEntries();
+ Float_t mPy2vP2Out = his->GetMean();
+ Float_t rPy2vP2Out = his->GetRMS();
+ delete his;
+
+ his = fHisPy3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
+ Float_t ePy3vP2IO = his->GetEntries();
+ Float_t mPy3vP2IO = his->GetMean();
+ Float_t rPy3vP2IO = his->GetRMS();
+ delete his;
+
+ //
+ //
+ his = fHisPz1vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - line
+ Float_t ePz1vP0 = his->GetEntries();
+ Float_t mPz1vP0 = his->GetMean();
+ Float_t rPz1vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPz2vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - parabola
+ Float_t ePz2vP0 = his->GetEntries();
+ Float_t mPz2vP0 = his->GetMean();
+ Float_t rPz2vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPz3vP0->ProjectionY("aaa",id+1,id+1); //-> delta z P0outer-P0inner - common parabola
+ Float_t ePz3vP0 = his->GetEntries();
+ Float_t mPz3vP0 = his->GetMean();
+ Float_t rPz3vP0 = his->GetRMS();
+ delete his;
+
+ his = fHisPz1vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - line
+ Float_t ePz1vP1 = his->GetEntries();
+ Float_t mPz1vP1 = his->GetMean();
+ Float_t rPz1vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPz2vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - parabola
+ Float_t ePz2vP1 = his->GetEntries();
+ Float_t mPz2vP1 = his->GetMean();
+ Float_t rPz2vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPz3vP1->ProjectionY("aaa",id+1,id+1); //-> delta kz P1outer-P1inner - common parabola
+ Float_t ePz3vP1 = his->GetEntries();
+ Float_t mPz3vP1 = his->GetMean();
+ Float_t rPz3vP1 = his->GetRMS();
+ delete his;
+
+ his = fHisPz2vP2In->ProjectionY("aaa",id+1,id+1); //-> Curv P2inner - parabola
+ Float_t ePz2vP2In = his->GetEntries();
+ Float_t mPz2vP2In = his->GetMean();
+ Float_t rPz2vP2In = his->GetRMS();
+ delete his;
+
+ his = fHisPz2vP2Out->ProjectionY("aaa",id+1,id+1); //-> Curv P2outer - parabola
+ Float_t ePz2vP2Out = his->GetEntries();
+ Float_t mPz2vP2Out = his->GetMean();
+ Float_t rPz2vP2Out = his->GetRMS();
+ delete his;
+
+ his = fHisPz3vP2IO->ProjectionY("aaa",id+1,id+1); //-> Curv P2outerinner - common parabola
+ Float_t ePz3vP2IO = his->GetEntries();
+ Float_t mPz3vP2IO = his->GetMean();
+ Float_t rPz3vP2IO = his->GetRMS();
+ delete his;
+
+ // Fit absolute laser residuals
+ TH2F* histAbsY = (TH2F*)laser->fDeltaYresAbs[id];
+ TH2F* histAbsZ = (TH2F*)laser->fDeltaZresAbs[id];
+
+ Int_t secInner = TMath::Nint(fBeamSectorInner[id]);
+ Int_t secOuter = TMath::Nint(fBeamSectorOuter[id]);
+
+ TVectorD vecX(159); // X
+ TVectorD vecY(159); // Y
+ TVectorD vecR(159); // R
+ TVectorD vecDY(159); // absolute residuals in Y
+ TVectorD vecDZ(159); // absolute residuals in Z
+ TVectorD vecN(159); // number of clusters
+ TVectorD vecEy(159); //error y
+ TVectorD vecEz(159); //error z
+ TVectorD vecPhi(159); // local tangent
+ TVectorD vecPhiR(159); // local tangent
+ // magnetic field integrals
+ TVectorD vecIBR(159); // radial
+ TVectorD vecIBRPhi(159); // r-phi
+ TVectorD vecIBLX(159); // local x
+ TVectorD vecIBLY(159); // local y
+ TVectorD vecIBGX(159); // local x
+ TVectorD vecIBGY(159); // local y
+ TVectorD vecIBZ(159); // z
+ //
+ for (Int_t irow=0;irow<159;irow++){
+ vecIBR[irow]=0;
+ vecIBRPhi[irow]=0;
+ vecIBLX[irow]=0;
+ vecIBLY[irow]=0;
+ vecIBGX[irow]=0;
+ vecIBGY[irow]=0;
+ vecIBZ[irow]=0;
+ Double_t gx =(*(ltrp->fVecGX))[irow];
+ Double_t gy =(*(ltrp->fVecGY))[irow];
+ Int_t lsec =TMath::Nint((*(ltrp->fVecSec))[irow]);
+ Double_t ca =TMath::Cos(TMath::Pi()*(lsec+0.5)/9.);
+ Double_t sa =TMath::Sin(TMath::Pi()*(lsec+0.5)/9.);
+ xyz[2]=(*(ltrp->fVecGZ))[irow];
+ xyz[0]=TMath::Sqrt(gx*gx+gy*gy);
+ xyz[1]=TMath::ATan2(gy,gx);
+ Double_t gxyz[3]={gx,gy,(*(ltrp->fVecGZ))[irow]};
+ if (magF){
+ magF->GetTPCIntCyl(xyz,bxyz);
+ magF->GetTPCInt(gxyz,bgxyz);
+ vecIBR[irow]=bxyz[0];
+ vecIBRPhi[irow]=bxyz[1];
+ //
+ vecIBGX[irow]=bgxyz[0];
+ vecIBGY[irow]=bgxyz[1];
+ //
+ vecIBLX[irow]= bgxyz[0]*ca+bgxyz[1]*sa;
+ vecIBLY[irow]= -bgxyz[0]*sa+bgxyz[1]*ca;
+ //
+
+ vecIBZ[irow]=bxyz[2];
+ }
+ }
+
+
+ lfabsyInner.ClearPoints();
+ lfabszInner.ClearPoints();
+ lfabsyOuter.ClearPoints();
+ lfabszOuter.ClearPoints();
+ // dummy fit values
+ Int_t nClInY = 0;
+ Float_t yAbsInOffset = -100;
+ Float_t yAbsInSlope = -100;
+ Float_t yAbsInDeltay = -100;
+ Int_t nClInZ = 0;
+ Float_t zAbsInOffset = -100;
+ Float_t zAbsInSlope = -100;
+ Float_t zAbsInDeltaz = -100;
+ Int_t nClOutY = 0;
+ Float_t yAbsOutOffset = -100;
+ Float_t yAbsOutSlope = -100;
+ Float_t yAbsOutDeltay = -100;
+ Int_t nClOutZ = 0;
+ Float_t zAbsOutOffset = -100;
+ Float_t zAbsOutSlope = -100;
+ Float_t zAbsOutDeltaz = -100;
+
+ Float_t lasTanPhiLocIn = -100;
+ Float_t lasTanPhiLocOut = -100;
+
+ if(histAbsY && histAbsY->GetEntries()>0) {
+
+ Double_t rotAngOut = 10;
+ Double_t rotAngIn = 10;
+ if((secInner%36)!=(secOuter%36))
+ rotAngIn += 20; // 30 degrees
+
+ // Calculate laser mirror X position in local frame
+ Double_t laserposOut =
+ TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngOut*TMath::DegToRad()));
+ Double_t laserposIn =
+ TMath::Abs(ltrp->GetX()*TMath::Cos(rotAngIn*TMath::DegToRad()));
+
+ // Calculate laser tan phi in local frame
+ lasTanPhiLocOut = TMath::ASin(ltrp->GetSnp());
+ if(lasTanPhiLocOut<0) {
+ lasTanPhiLocIn = lasTanPhiLocOut - rotAngIn*TMath::DegToRad();
+ lasTanPhiLocOut -= rotAngOut*TMath::DegToRad();
+ } else {
+
+ lasTanPhiLocIn = lasTanPhiLocOut + rotAngIn*TMath::DegToRad();
+ lasTanPhiLocOut += rotAngOut*TMath::DegToRad();
+ }
+
+ lasTanPhiLocIn = TMath::Tan(lasTanPhiLocIn);
+ lasTanPhiLocOut = TMath::Tan(lasTanPhiLocOut);
+
+ TProfile* yprof = histAbsY->ProfileX("yprof");
+ TProfile* zprof = histAbsZ->ProfileX("zprof");
+
+ for(Int_t bin = 1; bin<=159; bin++) {
+
+ if(yprof->GetBinEntries(bin)<5&&
+ zprof->GetBinEntries(bin)<5) {
+ continue;
+ }
+
+ // There is a problem in defining inner and outer sectors for
+ // the outer beams (0 and 6) where both sectors are OROCs. To
+ // make sure there is no overlap row 94 to 99 are cutted.
+ if(((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin>=95&&bin<=100)
+ continue;
+
+ Int_t row = (bin-1);
+ if(row>62)
+ row -= 63;
+
+ Bool_t isOuter = kTRUE;
+ Int_t sector = TMath::Nint(fBeamSectorOuter[id]);
+
+ if(bin<=62 ||
+ (((ltrp->GetBeam()==0)||(ltrp->GetBeam()==6))&&bin<=95)) {
+
+ isOuter = kFALSE;
+ sector = TMath::Nint(fBeamSectorInner[id]);
+ }
+
+
+ Double_t x = tpcparam->GetPadRowRadii(sector, row); // slope
+ vecN[bin-1] =yprof->GetBinEntries(bin);
+ vecEy[bin-1]=yprof->GetBinError(bin);
+ vecEz[bin-1]=zprof->GetBinError(bin);
+ vecX[bin-1] = x;
+ vecDY[bin-1] = yprof->GetBinContent(bin);
+ vecDZ[bin-1] = zprof->GetBinContent(bin);
+ if (bin>0&&bin<159){
+ //
+ //truncated mean - skip first and the last pad row
+ //
+ Int_t firstBin=TMath::Max(bin-5,0);
+ Int_t lastBin =TMath::Min(bin+5,158);
+ histAbsY->GetXaxis()->SetRangeUser(firstBin,lastBin);
+ histAbsY->GetYaxis()->SetRangeUser(-2,2);
+ vecEy[bin-1]=histAbsY->GetRMS(2);
+ vecDY[bin-1]=histAbsY->GetMean(2);
+ histAbsY->GetXaxis()->SetRangeUser(firstBin+2,lastBin-2);//use+-2 bins
+ histAbsY->GetYaxis()->SetRangeUser(vecDY[bin-1]-4*vecEy[bin-1],vecDY[bin-1]+4*vecEy[bin-1]);
+ if (yprof->GetBinEntries(bin-1)>0) vecEy[bin-1]=histAbsY->GetRMS(2)/TMath::Sqrt(yprof->GetBinEntries(bin-1));
+ vecDY[bin-1]=histAbsY->GetMean(2);
+ }
+
+ if(!isOuter) { // inner
+ vecPhi[bin-1]=lasTanPhiLocIn;
+ // Calculate local y from residual and database
+ Double_t y = fBeamSlopeYInner[id]*x + fBeamOffsetYInner[id]
+ + vecDY[bin-1];
+ vecY[bin-1] = y;
+ Double_t r = TMath::Sqrt(x*x + y*y);
+ vecR[bin-1] = r;
+ // Find angle between laser vector and R vector
+ // cos angle = (x,y).(1,fBeamSlopeYInner)/R/sqrt(1+fBeamSlopeYInner**2)
+ Double_t cosPhi = x + y*fBeamSlopeYInner[id];
+ cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYInner[id]*fBeamSlopeYInner[id]);
+ vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
+ if(lasTanPhiLocIn<0)
+ vecPhiR[bin-1]*=-1; // to have the same sign
+
+ if(yprof->GetBinEntries(bin)>=10) {
+ lfabsyInner.AddPoint(&x, yprof->GetBinContent(bin),
+ TMath::Max(yprof->GetBinError(bin), 0.001));
+ }
+ if(zprof->GetBinEntries(bin)>=10) {
+ lfabszInner.AddPoint(&x, zprof->GetBinContent(bin),
+ TMath::Max(zprof->GetBinError(bin), 0.001));
+ }
+ } else { // outer
+ vecPhi[bin-1]=lasTanPhiLocOut;
+ // Calculate local y from residual and database
+ Double_t y = fBeamSlopeYOuter[id]*x + fBeamOffsetYOuter[id]
+ + vecDY[bin-1];
+ vecY[bin-1] = y;
+ Double_t r = TMath::Sqrt(x*x + y*y);
+ vecR[bin-1] = r;
+
+ Double_t cosPhi = x + y*fBeamSlopeYOuter[id];
+ cosPhi /= r*TMath::Sqrt(1+fBeamSlopeYOuter[id]*fBeamSlopeYOuter[id]);
+ vecPhiR[bin-1] = TMath::Tan(TMath::ACos(cosPhi));
+ if(lasTanPhiLocOut<0)
+ vecPhiR[bin-1]*=-1; // to have the same sign
+
+ if(yprof->GetBinEntries(bin)>=10) {
+ lfabsyOuter.AddPoint(&x, yprof->GetBinContent(bin),
+ TMath::Max(yprof->GetBinError(bin), 0.001));
+ }
+ if(zprof->GetBinEntries(bin)>=10) {
+ lfabszOuter.AddPoint(&x, zprof->GetBinContent(bin),
+ TMath::Max(zprof->GetBinError(bin), 0.001));
+ }
+ }
+ // global position
+
+ }
+
+ delete yprof; delete zprof;
+
+
+ // Fit laser abs residuals with linear robust fit (exclude 5% outliers)
+ nClInY = lfabsyInner.GetNpoints();
+ if(lfabsyInner.GetNpoints()>10) {
+ lfabsyInner.EvalRobust(0.95);
+
+ TVectorD result(2);
+ lfabsyInner.GetParameters(result);
+ yAbsInOffset = result[0];
+ yAbsInSlope = result[1];
+ yAbsInDeltay = yAbsInSlope*laserposIn + yAbsInOffset;
+ }
+ nClInZ = lfabszInner.GetNpoints();
+ if(lfabszInner.GetNpoints()>10) {
+ lfabszInner.EvalRobust(0.95);
+
+ TVectorD result(2);
+ lfabszInner.GetParameters(result);
+ zAbsInOffset = result[0];
+ zAbsInSlope = result[1];
+ zAbsInDeltaz = zAbsInSlope*laserposIn + zAbsInOffset;
+ }
+ nClOutY = lfabsyOuter.GetNpoints();
+ if(lfabsyOuter.GetNpoints()>10) {
+ lfabsyOuter.EvalRobust(0.95);
+
+ TVectorD result(2);
+ lfabsyOuter.GetParameters(result);
+ yAbsOutOffset = result[0];
+ yAbsOutSlope = result[1];
+ yAbsOutDeltay = yAbsOutSlope*laserposOut + yAbsOutOffset;
+ }
+ nClOutZ = lfabszOuter.GetNpoints();
+ if(lfabszOuter.GetNpoints()>10) {
+ lfabszOuter.EvalRobust(0.95);
+
+ TVectorD result(2);
+ lfabszOuter.GetParameters(result);
+ zAbsOutOffset = result[0];
+ zAbsOutSlope = result[1];
+ zAbsOutDeltaz = zAbsOutSlope*laserposOut + zAbsOutOffset;
+ }
+ }
+
+
+ Int_t itime=-1;
+ Float_t coverIA=AliTPCcalibDB::GetCoverVoltage(run,0,itime);
+ Float_t coverIC=AliTPCcalibDB::GetCoverVoltage(run,18,itime);
+ Float_t coverOA=AliTPCcalibDB::GetCoverVoltage(run,36,itime);
+ Float_t coverOC=AliTPCcalibDB::GetCoverVoltage(run,54,itime);
+ //
+ Float_t skirtA=AliTPCcalibDB::GetSkirtVoltage(run,0,itime);
+ Float_t skirtC=AliTPCcalibDB::GetSkirtVoltage(run,18,itime);
+ //
+ Float_t ggOffA=AliTPCcalibDB::GetGGoffsetVoltage(run,0,itime);
+ Float_t ggOffC=AliTPCcalibDB::GetGGoffsetVoltage(run,18,itime);
+
+ //
+ (*pcstream)<<"Mean"<<
+ "run="<<run<< //
+ //voltages
+ "UcIA="<<coverIA<<
+ "UcIC="<<coverIC<<
+ "UcOA="<<coverOA<<
+ "UcOC="<<coverOC<<
+ "UsA="<<skirtA<<
+ "UsC="<<skirtC<<
+ "UggA="<<ggOffA<<
+ "UggC="<<ggOffC<<
+ //
+ "isOK="<<isOK<< //
+ "id="<<id<< // track id
+ "entries="<<entries<< // number of entries
+ "bz="<<bfield<< // bfield
+ "LTr.="<<ltrp<< // refernece track
+ //
+ "mnclIn="<<mnclIn<< // mean number of clusters in inner
+ "mnclOut="<<mnclOut<< // mean number of clusters in outer
+ "mnclIO="<<mnclIO<< // mean number of clusters in inner+outer
+ "mLclIn="<<mLclIn<< // mean number of clusters in inner
+ "mLclOut="<<mLclOut<< // mean number of clusters in outer
+ "mLclIO="<<mLclIO<< // mean number of clusters in inner+outer
+ "mdEdx="<<mdEdx<< // mean dEdx
+ "edZfit="<<edZfit<< // entries z fit
+ "mdZfit="<<mdZfit<< // mean z fit
+ "rdZfit="<<rdZfit<< // RMS z fit
+ //
+ //
+ "mChi2YIn1="<<mChi2YIn1<< //->chi2 y inner - line
+ "mChi2YOut1="<<mChi2YOut1<< //->chi2 y inner - line
+ "mChi2YIn2="<<mChi2YIn2<< //->chi2 y inner - parabola
+ "mChi2YOut2="<<mChi2YOut2<< //->chi2 y inner - parabola
+ "mChi2YIO1="<<mChi2YIO1<< //->chi2 y IO - common
+ "mChi2ZIn1="<<mChi2ZIn1<< //->chi2 z inner - line
+ "mChi2ZOut1="<<mChi2ZOut1<< //->chi2 z inner - line
+ "mChi2ZIn2="<<mChi2ZIn2<< //->chi2 z inner - parabola
+ "mChi2ZOut2="<<mChi2ZOut2<< //->chi2 z inner - parabola
+ "mChi2ZIO1="<<mChi2ZIO1<< //->chi2 z IO - common
+ //
+ //
+ "ePy1vP0="<<ePy1vP0<<
+ "mPy1vP0="<<mPy1vP0<<
+ "rPy1vP0="<<rPy1vP0<<
+ "ePy2vP0="<<ePy2vP0<<
+ "mPy2vP0="<<mPy2vP0<<
+ "rPy2vP0="<<rPy2vP0<<
+ "ePy3vP0="<<ePy3vP0<<
+ "mPy3vP0="<<mPy3vP0<<
+ "rPy3vP0="<<rPy3vP0<<
+ "ePy1vP1="<<ePy1vP1<<
+ "mPy1vP1="<<mPy1vP1<<
+ "rPy1vP1="<<rPy1vP1<<
+ "ePy2vP1="<<ePy2vP1<<
+ "mPy2vP1="<<mPy2vP1<<
+ "rPy2vP1="<<rPy2vP1<<
+ "ePy3vP1="<<ePy3vP1<<
+ "mPy3vP1="<<mPy3vP1<<
+ "rPy3vP1="<<rPy3vP1<<
+ "ePy2vP2In="<<ePy2vP2In<<
+ "mPy2vP2In="<<mPy2vP2In<<
+ "rPy2vP2In="<<rPy2vP2In<<
+ "ePy2vP2Out="<<ePy2vP2Out<<
+ "mPy2vP2Out="<<mPy2vP2Out<<
+ "rPy2vP2Out="<<rPy2vP2Out<<
+ "ePy3vP2IO="<<ePy3vP2IO<<
+ "mPy3vP2IO="<<mPy3vP2IO<<
+ "rPy3vP2IO="<<rPy3vP2IO<<
+ //
+ //
+ "ePz1vP0="<<ePz1vP0<<
+ "mPz1vP0="<<mPz1vP0<<
+ "rPz1vP0="<<rPz1vP0<<
+ "ePz2vP0="<<ePz2vP0<<
+ "mPz2vP0="<<mPz2vP0<<
+ "rPz2vP0="<<rPz2vP0<<
+ "ePz3vP0="<<ePz3vP0<<
+ "mPz3vP0="<<mPz3vP0<<
+ "rPz3vP0="<<rPz3vP0<<
+ "ePz1vP1="<<ePz1vP1<<
+ "mPz1vP1="<<mPz1vP1<<
+ "rPz1vP1="<<rPz1vP1<<
+ "ePz2vP1="<<ePz2vP1<<
+ "mPz2vP1="<<mPz2vP1<<
+ "rPz2vP1="<<rPz2vP1<<
+ "ePz3vP1="<<ePz3vP1<<
+ "mPz3vP1="<<mPz3vP1<<
+ "rPz3vP1="<<rPz3vP1<<
+ "ePz2vP2In="<<ePz2vP2In<<
+ "mPz2vP2In="<<mPz2vP2In<<
+ "rPz2vP2In="<<rPz2vP2In<<
+ "ePz2vP2Out="<<ePz2vP2Out<<
+ "mPz2vP2Out="<<mPz2vP2Out<<
+ "rPz2vP2Out="<<rPz2vP2Out<<
+ "ePz3vP2IO="<<ePz3vP2IO<<
+ "mPz3vP2IO="<<mPz3vP2IO<<
+ "rPz3vP2IO="<<rPz3vP2IO<<
+ //
+ //
+ //
+ "lx0="<<lxyz[0]<< // reference x
+ "lx1="<<lxyz[1]<< // reference y
+ "lx2="<<lxyz[2]<< // refernece z
+ "lpx0="<<lpxyz[0]<< // reference x
+ "lpx1="<<lpxyz[1]<< // reference y
+ "lpx2="<<lpxyz[2]<< // refernece z
+ //
+ "msig="<<meanS<<
+ //
+ "mphi="<<meanphi<< //
+ "rmsphi="<<rmsphi<< //
+ "gphi1="<<gphi1<<
+ "gphi2="<<gphi2<<
+ //
+ "mphiP="<<meanphiP<< //
+ "rmsphiP="<<rmsphiP<< //
+ "gphiP1="<<gphiP1<<
+ "gphiP2="<<gphiP2<<
+ //
+ "meanZ="<<meanZ<<
+ "rmsZ="<<rmsZ<<
+ "gz1="<<gz1<<
+ "gz2="<<gz2<<
+ //
+ "gp31="<<gp31<< //gaus mean - tgl
+ "gp32="<<gp32<< //gaus rms -tgl
+ "meanp3="<<meanp3<< //mean - tgl
+ "rmsp3="<<rmsp3<< //rms -tgl
+ "gp41="<<gp41<< //gaus mean - P4
+ "gp42="<<gp42<< //gaus rms - P4
+ "meanp4="<<meanp4<< //mean - P4
+ "rmsp4="<<rmsp4<< //rms - P4
+ // Parameters from abs res analysis
+ "SecIn="<<secInner<< // inner sector
+ "SecOut="<<secOuter<< // outer sector
+ "lasTanPhiLocIn="<<lasTanPhiLocIn<< // laser tan phi in local frame (inner)
+ "lasTanPhiLocOut="<<lasTanPhiLocOut<<// laser tan phi in local frame (outer)
+ "ibr.="<<&vecIBR<< // radial filed integral
+ "ibrphi.="<<&vecIBRPhi<< // r=phifiled integral
+ "ibr.="<<&vecIBR<< // radial filed integral
+ "ibz.="<<&vecIBZ<< // z filed integral
+ //
+ "iblx.="<<&vecIBLX<< // local bx integral
+ "ibly.="<<&vecIBLY<< // local by integral
+ "ibgx.="<<&vecIBGX<< // global bx integral
+ "ibgy.="<<&vecIBGY<< // global by integral
+ //
+ "X.="<<&vecX<< // local x
+ "Y.="<<&vecY<< // local y
+ "R.="<<&vecR<< // radius
+ "dY.="<<&vecDY<< // abs y residuals
+ "dZ.="<<&vecDZ<< // abs z residuals
+ "eY.="<<&vecEy<< // error of y residuals
+ "eZ.="<<&vecEz<< // error of z residuals
+ "kY.="<<&vecPhi<< // local tangent y (fixed for sector)
+ "kYR.="<<&vecPhiR<< // tangent between laser and R vector (varies inside sector)
+ "nCl.="<<&vecN<< // number of clusters
+ //
+ "nClInY="<<nClInY<< // Number of clusters for inner y
+ "yAbsInOffset="<<yAbsInOffset<< // fitted offset absres (inner y)
+ "yAbsInSlope="<<yAbsInSlope << // fitted slope absres (inner y)
+ "yAbsInDeltay="<<yAbsInDeltay<< // fitted laser offset absres (inner y)
+ "nClInZ="<<nClInZ<< // Number of clusters for inner z
+ "zAbsInOffset="<<zAbsInOffset<< // fitted offset absres (inner z)
+ "zAbsInSlope="<<zAbsInSlope << // fitted slope absres (inner z)
+ "zAbsInDeltaz="<<zAbsInDeltaz<< // fitted laser offset absres (inner z)
+ //
+ "nClOutY="<<nClOutY<< // Number of clusters for outer y
+ "yAbsOutOffset="<<yAbsOutOffset<< // fitted offset absres (outer y)
+ "yAbsOutSlope="<<yAbsOutSlope << // fitted slope absres (outer y)
+ "yAbsOutDeltay="<<yAbsOutDeltay<< // fitted laser offset absres (outer y)
+ "nClOutZ="<<nClOutZ<< // Number of clusters for outer z
+ "zAbsOutOffset="<<zAbsOutOffset<< // fitted offset absres (outer z)
+ "zAbsOutSlope="<<zAbsOutSlope << // fitted slope absres (outer z)
+ "zAbsOutDeltaz="<<zAbsOutDeltaz<< // fitted laser offset absres (outer z)
+ //
+ "\n";
+ }
+ delete pcstream;
+ /*
+ Browse the content
+ TFile fmean("laserMean.root");
+
+
+ */
+
+
+}
+
+
+
+void AliTPCcalibLaser::DumpScanInfo(TTree * chain, const char * cutUser){
+ //
+ //
+ //
+ TTreeSRedirector *pcstream = new TTreeSRedirector("laserScan.root");
+ TFile * f = pcstream->GetFile();
+ f->mkdir("dirphi");
+ f->mkdir("dirphiP");
+ f->mkdir("dirZ");
+ TF1 fp("p1","pol1");
+ //
+ //
+ char cut[1000];
+ char grname[1000];
+ char grnamefull[1000];
+
+ Double_t mphi[100];
+ Double_t mphiP[100];
+ Double_t smphi[100];
+ Double_t smphiP[100];
+ Double_t mZ[100];
+ Double_t smZ[100];
+ Double_t bz[100];
+ Double_t sbz[100];
+ // fit parameters
+ Double_t pphi[3];
+ Double_t pphiP[3];
+ Double_t pmZ[3];
+
+ //
+ for (Int_t id=0; id<336; id++){
+ // id =205;
+ snprintf(cut,1000,"fId==%d&&%s",id,cutUser);
+ Int_t entries = chain->Draw("bz",cut,"goff");
+ if (entries<3) continue;
+ AliTPCLaserTrack *ltrp = 0;
+ if (!AliTPCLaserTrack::GetTracks()) AliTPCLaserTrack::LoadTracks();
+ ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+
+ chain->Draw("bz",cut,"goff");
+ memcpy(bz, chain->GetV1(), entries*sizeof(Double_t));
+ chain->Draw("0.01*abs(bz)+0.02",cut,"goff");
+ memcpy(sbz, chain->GetV1(), entries*sizeof(Double_t));
+ //
+ chain->Draw("gphi1",cut,"goff");
+ memcpy(mphi, chain->GetV1(), entries*sizeof(Double_t));
+ chain->Draw("0.05*abs(mphi)+abs(gphi2)*0.5+0.05",cut,"goff");
+ memcpy(smphi, chain->GetV1(), entries*sizeof(Double_t));
+ //
+ chain->Draw("gphiP1",cut,"goff");
+ memcpy(mphiP, chain->GetV1(), entries*sizeof(Double_t));
+ chain->Draw("0.05*abs(mphiP)+abs(gphiP2)*0.5+0.001",cut,"goff");
+ memcpy(smphiP, chain->GetV1(), entries*sizeof(Double_t));
+ //
+ chain->Draw("gz1",cut,"goff");
+ memcpy(mZ, chain->GetV1(), entries*sizeof(Double_t));
+ chain->Draw("0.01*abs(meanZ)+abs(gz2)*0.5+0.1",cut,"goff");
+ memcpy(smZ, chain->GetV1(), entries*sizeof(Double_t));
+ //
+ //
+ snprintf(grnamefull,1000,"Side_%d_Bundle_%d_Rod_%d_Beam_%d",
+ ltrp->GetSide(), ltrp->GetBundle(), ltrp->GetRod(), ltrp->GetBeam());
+ // store data
+ // phi
+ f->cd("dirphi");
+ Float_t phiB0 =0;
+ Float_t phiPB0=0;
+ Float_t zB0=0;
+ for (Int_t ientry=0; ientry<entries; ientry++){
+ if (TMath::Abs(bz[ientry])<0.05){
+ phiB0 = mphi[ientry];
+ phiPB0 = mphiP[ientry];
+ zB0 = mZ[ientry];
+ }
+ }
+ TGraphErrors *grphi = new TGraphErrors(entries,bz,mphi,sbz,smphi);
+ grphi->Draw("a*");
+ grphi->Fit(&fp);
+ pphi[0] = fp.GetParameter(0); // offset
+ pphi[1] = fp.GetParameter(1); // slope
+ pphi[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
+ snprintf(grname,1000,"phi_id%d",id);
+ grphi->SetName(grname); grphi->SetTitle(grnamefull);
+ grphi->GetXaxis()->SetTitle("b_{z} (T)");
+ grphi->GetYaxis()->SetTitle("#Delta r#phi (cm)");
+ grphi->SetMaximum(1.2);
+ grphi->SetMinimum(-1.2);
+ grphi->Draw("a*");
+
+ grphi->Write();
+ gPad->SaveAs(Form("pic/phi/phi_%s.gif",grnamefull));
+ // phiP
+ f->cd("dirphiP)");
+ TGraphErrors *grphiP = new TGraphErrors(entries,bz,mphiP,sbz,smphiP);
+ grphiP->Draw("a*");
+ grphiP->Fit(&fp);
+ pphiP[0] = fp.GetParameter(0); // offset
+ pphiP[1] = fp.GetParameter(1); // slope
+ pphiP[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
+ snprintf(grname,1000,"phiP_id%d",id);
+ grphiP->SetName(grname); grphiP->SetTitle(grnamefull);
+ grphiP->GetXaxis()->SetTitle("b_{z} (T)");
+ grphiP->GetYaxis()->SetTitle("#Delta #phi (rad)");
+ grphiP->SetMaximum(pphiP[0]+0.005);
+ grphiP->SetMinimum(pphiP[0]-0.005);
+
+ gPad->SaveAs(Form("pic/phiP/phiP_%s.gif",grnamefull));
+ grphiP->Write();
+ //
+ //Z
+ f->cd("dirZ");
+ TGraphErrors *grmZ = new TGraphErrors(entries,bz,mZ,sbz,smZ);
+ grmZ->Draw("a*");
+ grmZ->Fit(&fp);
+ pmZ[0] = fp.GetParameter(0); // offset
+ pmZ[1] = fp.GetParameter(1); // slope
+ pmZ[2] = TMath::Sqrt(fp.GetChisquare()/(entries-2.)); // normalized chi2
+ snprintf(grname,1000,"mZ_id%d",id);
+ grmZ->SetName(grname); grmZ->SetTitle(grnamefull);
+ grmZ->GetXaxis()->SetTitle("b_{z} (T)");
+ grmZ->GetYaxis()->SetTitle("#Delta z (cm)");
+
+ gPad->SaveAs(Form("pic/z/z_%s.gif",grnamefull));
+ grmZ->Write();
+ //
+ // P4
+ //
+
+ for (Int_t ientry=0; ientry<entries; ientry++){
+ (*pcstream)<<"Mean"<<
+ "id="<<id<<
+ "LTr.="<<ltrp<<
+ "entries="<<entries<<
+ "bz="<<bz[ientry]<<
+ "lx0="<<lxyz[0]<< // reference x
+ "lx1="<<lxyz[1]<< // reference y
+ "lx2="<<lxyz[2]<< // refernece z
+ "lpx0="<<lpxyz[0]<< // reference x
+ "lpx1="<<lpxyz[1]<< // reference y
+ "lpx2="<<lpxyz[2]<< // refernece z
+ //values
+ "phiB0="<<phiB0<< // position shift at 0 field
+ "phiPB0="<<phiPB0<< // angular shift at 0 field
+ "zB0="<<zB0<< // z shift for 0 field
+ //
+ "gphi1="<<mphi[ientry]<< // mean - from gaus fit
+ "pphi0="<<pphi[0]<< // offset
+ "pphi1="<<pphi[1]<< // slope
+ "pphi2="<<pphi[2]<< // norm chi2
+ //
+ "gphiP1="<<mphiP[ientry]<< // mean - from gaus fit
+ "pphiP0="<<pphiP[0]<< // offset
+ "pphiP1="<<pphiP[1]<< // slope
+ "pphiP2="<<pphiP[2]<< // norm chi2
+ //
+ "gz1="<<mZ[ientry]<<
+ "pmZ0="<<pmZ[0]<< // offset
+ "pmZ1="<<pmZ[1]<< // slope
+ "pmZ2="<<pmZ[2]<< // norm chi2
+ "\n";
+ }
+ }
+
+ delete pcstream;
+
+}
+
+
+void AliTPCcalibLaser::Analyze(){
+ //
+ //
+ //
+}
+
+
+Long64_t AliTPCcalibLaser::Merge(TCollection *li) {
+
+ TIterator* iter = li->MakeIterator();
+ AliTPCcalibLaser* cal = 0;
+ static Int_t counter0=0;
+ while ((cal = (AliTPCcalibLaser*)iter->Next())) {
+ if (!cal->InheritsFrom(AliTPCcalibLaser::Class())) {
+ // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
+ return -1;
+ }
+ printf("Marging number %d\n", counter0);
+ counter0++;
+ //
+ MergeFitHistos(cal);
+ TH1F *h=0x0;
+ TH1F *hm=0x0;
+ TH2F *h2=0x0;
+ TH2F *h2m=0x0;
+ // TProfile *hp=0x0;
+ //TProfile *hpm=0x0;
+
+ for (Int_t id=0; id<336; id++){
+ // merge fDeltaZ histograms
+ hm = (TH1F*)cal->fDeltaZ.At(id);
+ h = (TH1F*)fDeltaZ.At(id);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fDeltaZ.AddAt(h,id);
+ }
+ if (h && hm) h->Add(hm);
+
+ // merge fP3 histograms
+ hm = (TH1F*)cal->fDeltaP3.At(id);
+ h = (TH1F*)fDeltaP3.At(id);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fDeltaP3.AddAt(h,id);
+ }
+ if (h && hm) h->Add(hm);
+ // merge fP4 histograms
+ hm = (TH1F*)cal->fDeltaP4.At(id);
+ h = (TH1F*)fDeltaP4.At(id);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fDeltaP4.AddAt(h,id);
+ }
+ if (h&&hm) h->Add(hm);
+
+ //
+ // merge fDeltaPhi histograms
+ hm = (TH1F*)cal->fDeltaPhi.At(id);
+ h = (TH1F*)fDeltaPhi.At(id);
+ if (!h &&hm &&hm->GetEntries()>0) {
+ h= (TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fDeltaPhi.AddAt(h,id);
+ }
+ if (h&&hm) h->Add(hm);
+ // merge fDeltaPhiP histograms
+ hm = (TH1F*)cal->fDeltaPhiP.At(id);
+ h = (TH1F*)fDeltaPhiP.At(id);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fDeltaPhiP.AddAt(h,id);
+ }
+ if (h&&hm) h->Add(hm);
+ // merge fSignals histograms
+ hm = (TH1F*)cal->fSignals.At(id);
+ h = (TH1F*)fSignals.At(id);
+ if (!h&&hm &&hm->GetEntries()>0) {
+ h=(TH1F*)hm->Clone();
+ h->SetDirectory(0);
+ fSignals.AddAt(h,id);
+ }
+ if (h&&hm) h->Add(hm);
+ //
+ //
+ // merge ProfileY histograms -0
+ h2m = (TH2F*)cal->fDeltaYres.At(id);
+ h2 = (TH2F*)fDeltaYres.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+ //
+ h2m = (TH2F*)cal->fDeltaZres.At(id);
+ h2 = (TH2F*)fDeltaZres.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+ // merge ProfileY histograms - 2
+ h2m = (TH2F*)cal->fDeltaYres2.At(id);
+ h2 = (TH2F*)fDeltaYres2.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+ //
+ h2m = (TH2F*)cal->fDeltaZres2.At(id);
+ h2 = (TH2F*)fDeltaZres2.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+
+ // merge ProfileY histograms - abs
+ h2m = (TH2F*)cal->fDeltaYresAbs.At(id);
+ h2 = (TH2F*)fDeltaYresAbs.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+ if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaYresAbs.AddAt(h2,id);}
+ h2m = (TH2F*)cal->fDeltaZresAbs.At(id);
+ h2 = (TH2F*)fDeltaZresAbs.At(id);
+ if (h2m&&h2) h2->Add(h2m);
+ if (h2m&&!h2) { h2=(TH2F*)h2m->Clone(); h2->SetDirectory(0); fDeltaZresAbs.AddAt(h2,id);}
+ // merge ProfileY histograms - 3
+ //h2m = (TH2F*)cal->fDeltaYres3.At(id);
+ //h2 = (TH2F*)fDeltaYres3.At(id);
+ //if (h2m) h2->Add(h2m);
+ //
+ //h2m = (TH2F*)cal->fDeltaZres3.At(id);
+ //h2 = (TH2F*)fDeltaZres3.At(id);
+ //if (h2m) h->Add(h2m);
+ //
+ //
+ }
+ }
+ return 0;
+}
+
+void AliTPCcalibLaser::MakeFitHistos(){
+ //
+ // Make a fit histograms
+ //
+ // Number of clusters
+ //
+ //TH2F *fHisNclIn; //->Number of clusters inner
+ //TH2F *fHisNclOut; //->Number of clusters outer
+ //TH2F *fHisNclIO; //->Number of cluster inner outer
+ // TH2F *fHisdEdx; //->dEdx histo
+ fHisNclIn = new TH2F("HisNclIn","HisNclIn",336,0,336,64,10,64);
+ fHisNclOut = new TH2F("HisNclOut","HisNclOut",336,0,336,100,10,100);
+ fHisNclIO = new TH2F("HisNclIO","HisNclIO",336,0,336,160,10,160);
+ //
+ fHisLclIn = new TH2F("HisLclIn","HisLclIn",336,0,336,64,10,64);
+ fHisLclOut = new TH2F("HisLclOut","HisLclOut",336,0,336,100,10,150);
+ fHisLclIO = new TH2F("HisLclIO","HisLclIO",336,0,336,160,10,160);
+ //
+ fHisdEdx = new TH2F("HisdEdx","HisdEdx",336,0,336,160,1,50);
+ fHisdZfit = new TH2F("HisdZfit","HisdZfit",336,0,336,300,-0.3,0.3);
+
+ //
+ // Chi2
+ //
+ // TH2F *fHisChi2YIn1; //->chi2 y inner - line
+ // TH2F *fHisChi2YOut1; //->chi2 y inner - line
+ // TH2F *fHisChi2YIn2; //->chi2 y inner - parabola
+ // TH2F *fHisChi2YOut2; //->chi2 y inner - parabola
+ // TH2F *fHisChi2YIO1; //->chi2 y IO - common
+ fHisChi2YIn1 = new TH2F("Chi2YIn1","Chi2YIn1",336,0,336,500,0.001,0.5);
+ fHisChi2YOut1 = new TH2F("Chi2YOut1","Chi2YOut1",336,0,336,500,0.001,0.5);
+ fHisChi2YIn2 = new TH2F("Chi2YIn2","Chi2YIn2",336,0,336,500,0.001,0.5);
+ fHisChi2YOut2 = new TH2F("Chi2YOut2","Chi2YOut2",336,0,336,500,0.001,0.5);
+ fHisChi2YIO1 = new TH2F("Chi2YIO1","Chi2YIO1",336,0,336,500,0.001,0.5);
+ // TH2F *fHisChi2ZIn1; //->chi2 z inner - line
+ // TH2F *fHisChi2ZOut1; //->chi2 z inner - line
+ // TH2F *fHisChi2ZIn2; //->chi2 z inner - parabola
+ // TH2F *fHisChi2ZOut2; //->chi2 z inner - parabola
+ // TH2F *fHisChi2ZIO1; //->chi2 z IO - common
+ fHisChi2ZIn1 = new TH2F("Chi2ZIn1","Chi2ZIn1",336,0,336,500,0.001,0.5);
+ fHisChi2ZOut1 = new TH2F("Chi2ZOut1","Chi2ZOut1",336,0,336,500,0.001,0.5);
+ fHisChi2ZIn2 = new TH2F("Chi2ZIn2","Chi2ZIn2",336,0,336,500,0.001,0.5);
+ fHisChi2ZOut2 = new TH2F("Chi2ZOut2","Chi2ZOut2",336,0,336,500,0.001,0.5);
+ fHisChi2ZIO1 = new TH2F("Chi2ZIO1","Chi2ZIO1",336,0,336,500,0.001,0.5);
+ //
+ // Fit
+ //
+ //
+ // TH2F *fHisPy1vP0; //-> delta y P0outer-P0inner - line
+ // TH2F *fHisPy2vP0; //-> delta y P0outer-P0inner - parabola
+ // TH2F *fHisPy3vP0; //-> delta y P0outer-P0inner - common parabola
+ // TH2F *fHisPy1vP1; //-> delta ky P1outer-P1inner - line
+ // TH2F *fHisPy2vP1; //-> delta ky P1outer-P1inner - parabola
+ // TH2F *fHisPy3vP1; //-> delta ky P1outer-P1inner - common parabola
+ // TH2F *fHisPy2vP2In; //-> Curv P2inner - parabola
+ // TH2F *fHisPy2vP2Out; //-> Curv P2outer - parabola
+ // TH2F *fHisPy3vP2IO; //-> Curv P2outerinner - common parabola
+ fHisPy1vP0 = new TH2F("HisPy1vP0", "HisPy1vP0",336,0,336,500,-0.3,0.3);
+ fHisPy2vP0 = new TH2F("HisPy2vP0", "HisPy2vP0",336,0,336,500,-0.3,0.3);
+ fHisPy3vP0 = new TH2F("HisPy3vP0", "HisPy3vP0",336,0,336,500,-0.3,0.3);
+ fHisPy1vP1 = new TH2F("HisPy1vP1", "HisPy1vP1",336,0,336,500,-0.01,0.01);
+ fHisPy2vP1 = new TH2F("HisPy2vP1", "HisPy2vP1",336,0,336,500,-0.01,0.01);
+ fHisPy3vP1 = new TH2F("HisPy3vP1", "HisPy3vP1",336,0,336,500,-0.01,0.01);
+ fHisPy2vP2In = new TH2F("HisPy2vP2In", "HisPy2vP2In",336,0,336,500,-0.0003,0.0003);
+ fHisPy2vP2Out= new TH2F("HisPy2vP2Out","HisPy2vP2Out",336,0,336,500,-0.0003,0.0003);
+ fHisPy3vP2IO = new TH2F("HisPy3vP2IO", "HisPy3vP2IO",336,0,336,500,-0.0003,0.0003);
+ //
+ //
+ // TH2F *fHisPz1vP0; //-> delta z P0outer-P0inner - line
+ // TH2F *fHisPz2vP0; //-> delta z P0outer-P0inner - parabola
+ // TH2F *fHisPz3vP0; //-> delta z P0outer-P0inner - common parabola
+ // TH2F *fHisPz1vP1; //-> delta kz P1outer-P1inner - line
+ // TH2F *fHisPz2vP1; //-> delta kz P1outer-P1inner - parabola
+ // TH2F *fHisPz3vP1; //-> delta kz P1outer-P1inner - common parabola
+ // TH2F *fHisPz2vP2In; //-> Curv P2inner - parabola
+ // TH2F *fHisPz2vP2Out; //-> Curv P2outer - parabola
+ // TH2F *fHisPz3vP2IO; //-> Curv P2outerinner - common parabola
+ fHisPz1vP0 = new TH2F("HisPz1vP0", "HisPz1vP0",336,0,336,500,-0.3,0.3);
+ fHisPz2vP0 = new TH2F("HisPz2vP0", "HisPz2vP0",336,0,336,500,-0.3,0.3);
+ fHisPz3vP0 = new TH2F("HisPz3vP0", "HisPz3vP0",336,0,336,500,-0.3,0.3);
+ fHisPz1vP1 = new TH2F("HisPz1vP1", "HisPz1vP1",336,0,336,500,-0.01,0.01);
+ fHisPz2vP1 = new TH2F("HisPz2vP1", "HisPz2vP1",336,0,336,500,-0.01,0.01);
+ fHisPz3vP1 = new TH2F("HisPz3vP1", "HisPz3vP1",336,0,336,500,-0.01,0.01);
+ fHisPz2vP2In = new TH2F("HisPz2vP2In", "HisPz2vP2In",336,0,336,500,-0.0003,0.0003);
+ fHisPz2vP2Out= new TH2F("HisPz2vP2Out","HisPz2vP2Out",336,0,336,500,-0.0002,0.0002);
+ fHisPz3vP2IO = new TH2F("HisPz3vP2IO", "HisPz3vP2IO",336,0,336,500,-0.0002,0.0002);
+
+ fHisYAbsErrors = new TH1F("HisYabsErrors", "Errors per beam (y)", 336, 0, 336);
+ fHisZAbsErrors = new TH1F("HisZabsErrors", "Errors per beam (z)", 336, 0, 336);
+
+ fHisNclIn->SetDirectory(0); //->Number of clusters inner
+ fHisNclOut->SetDirectory(0); //->Number of clusters outer
+ fHisNclIO->SetDirectory(0); //->Number of cluster inner outer
+ fHisLclIn->SetDirectory(0); //->Level arm inner
+ fHisLclOut->SetDirectory(0); //->Level arm outer
+ fHisLclIO->SetDirectory(0); //->Number of cluster inner outer
+ fHisdEdx->SetDirectory(0); //->Number of cluster inner outer
+ fHisdZfit->SetDirectory(0); //->Number of cluster inner outer
+ //
+ //
+ fHisChi2YIn1->SetDirectory(0); //->chi2 y inner - line
+ fHisChi2YOut1->SetDirectory(0); //->chi2 y inner - line
+ fHisChi2YIn2->SetDirectory(0); //->chi2 y inner - parabola
+ fHisChi2YOut2->SetDirectory(0); //->chi2 y inner - parabola
+ fHisChi2YIO1->SetDirectory(0); //->chi2 y IO - common
+ fHisChi2ZIn1->SetDirectory(0); //->chi2 z inner - line
+ fHisChi2ZOut1->SetDirectory(0); //->chi2 z inner - line
+ fHisChi2ZIn2->SetDirectory(0); //->chi2 z inner - parabola
+ fHisChi2ZOut2->SetDirectory(0); //->chi2 z inner - parabola
+ fHisChi2ZIO1->SetDirectory(0); //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0->SetDirectory(0); //-> delta y P0outer-P0inner - line
+ fHisPy2vP0->SetDirectory(0); //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0->SetDirectory(0); //-> delta y P0outer-P0inner - common parabola
+ fHisPy1vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1->SetDirectory(0); //-> delta ky P1outer-P1inner - common parabola
+ fHisPy2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
+ fHisPy2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
+ fHisPy3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
+ //
+ //
+ fHisPz1vP0->SetDirectory(0); //-> delta z P0outer-P0inner - line
+ fHisPz2vP0->SetDirectory(0); //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0->SetDirectory(0); //-> delta z P0outer-P0inner - common parabola
+ fHisPz1vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1->SetDirectory(0); //-> delta kz P1outer-P1inner - common parabola
+ fHisPz2vP2In->SetDirectory(0); //-> Curv P2inner - parabola
+ fHisPz2vP2Out->SetDirectory(0); //-> Curv P2outer - parabola
+ fHisPz3vP2IO->SetDirectory(0); //-> Curv P2outerinner - common parabola
+
+ fHisYAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res y analysis
+ fHisZAbsErrors->SetDirectory(0); //-> total errors per beam in the abs res z analysis
+
+
+
+ //
+ //
+ //
+ for (Int_t id=0; id<336;id++){
+ TH2F *profy = (TH2F*)fDeltaYres.UncheckedAt(id);
+ TH2F *profz = (TH2F*)fDeltaZres.UncheckedAt(id);
+ //TH2F *profy2 = (TH2F*)fDeltaYres2.UncheckedAt(id);
+ TH2F *profy2 = 0;
+ TH2F *profz2 = 0;//(TH2F*)fDeltaZres2.UncheckedAt(id);
+ TH2F *profyabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
+ TH2F *profzabs = 0; //(TH2F*)fDeltaYresAbs.UncheckedAt(id);
+ // TH2F *profy3 = (TH2F*)fDeltaYres3.UncheckedAt(id);
+ //TH2F *profz3 = (TH2F*)fDeltaZres3.UncheckedAt(id);
+ if (!profy){
+ profy=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Linear",id),160,0,160,50,-0.5,0.5);
+ profy->SetDirectory(0);
+ fDeltaYres.AddAt(profy,id);
+ profy2=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d -Parabolic",id),160,0,160,50,-0.5,0.5);
+ profy2->SetDirectory(0);
+ fDeltaYres2.AddAt(profy2,id);
+ if(!fUseFixedDriftV)
+ profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+ else
+ profyabs=new TH2F(Form("pryabs%03d",id),Form("Y Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
+ profyabs->SetDirectory(0);
+ fDeltaYresAbs.AddAt(profyabs,id);
+ //profy3=new TH2F(Form("pry%03d",id),Form("Y Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
+ //profy3->SetDirectory(0);
+ //fDeltaYres3.AddAt(profy3,id);
+ }
+ if (!profz){
+ profz=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d Linear",id),160,0,160,50,-0.5,0.5);
+ profz->SetDirectory(0);
+ fDeltaZres.AddAt(profz,id);
+ profz2=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d - Parabolic",id),160,0,160,50,-0.5,0.5);
+ profz2->SetDirectory(0);
+ fDeltaZres2.AddAt(profz2,id);
+ if(!fUseFixedDriftV)
+ profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,100,-1.0,1.0); // has to be bigger based on earlier studies
+ else
+ profzabs=new TH2F(Form("przabs%03d",id),Form("Z Residuals for Laser Beam %03d -Absolute",id),160,0,160,200,-2.0,2.0); // has to be bigger based on earlier studies
+ profzabs->SetDirectory(0);
+ fDeltaZresAbs.AddAt(profzabs,id);
+ //profz3=new TH2F(Form("prz%03d",id),Form("Z Residuals for Laser Beam %03d- Parabolic2",id),160,0,160,100,-0.5,0.5);
+ //profz3->SetDirectory(0);
+ //fDeltaZres3.AddAt(profz3,id);
+ }
+ }
+ //
+ //
+ for (Int_t id=0; id<336;id++){
+ TH1F * hisdz = (TH1F*)fDeltaZ.At(id);
+ //TH1F * hisP3 = (TH1F*)fDeltaP3.At(id);
+ TH1F * hisP3 = 0;
+ TH1F * hisP4 = 0;
+
+ TH1F * hisdphi = 0;//(TH1F*)fDeltaPhi.At(id);
+ TH1F * hisdphiP = 0;//(TH1F*)fDeltaPhiP.At(id);
+ TH1F * hisSignal = 0; //(TH1F*)fSignals.At(id);
+
+ if (!hisdz){
+ hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10);
+ hisdz->SetDirectory(0);
+ fDeltaZ.AddAt(hisdz,id);
+
+ hisP3 = new TH1F(Form("hisPar3v%d",id),Form("hisPar3v%d",id),400,-0.06,0.06);
+ hisP3->SetDirectory(0);
+ fDeltaP3.AddAt(hisP3,id);
+ //
+ hisP4 = new TH1F(Form("hisPar4v%d",id),Form("hisPar4v%d",id),200,-0.06,0.06);
+ hisP4->SetDirectory(0);
+ fDeltaP4.AddAt(hisP4,id);
+
+ //
+ hisdphi = new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1);
+ hisdphi->SetDirectory(0);
+ fDeltaPhi.AddAt(hisdphi,id);
+ //
+ hisdphiP = new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01);
+ hisdphiP->SetDirectory(0);
+ fDeltaPhiP.AddAt(hisdphiP,id);
+ hisSignal = new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),100,0,300);
+ hisSignal->SetDirectory(0);
+ fSignals.AddAt(hisSignal,id);
+ }
+ }
+
+ SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2);
+ SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3);
+ SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0);
+ SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1);
+ //
+ // Make THnSparse
+ //
+ // id side rod bundle beam dP0 dP1 dP2 dP3 dP4 ncl dEdx
+ Int_t binsLaser[12]= {336, //id
+ 2, //side
+ 6, //rod
+ 4, //bundle
+ 7, //beam
+ 300, //dP0
+ 300, //dP1
+ 300, //dP2
+ 300, //dP3
+ 300, //dP4
+ 80, //ncl
+ 50}; //dEdx
+ Double_t xminLaser[12]= {0, //id
+ 0, //side
+ 0, //rod
+ 0, //bundle
+ 0, //beam
+ -1, //dP0
+ -1, //dP1
+ -0.01, //dP2
+ -0.01, //dP3
+ -0.1, //dP4
+ 0, //ncl
+ 0}; //sqrt dEdx
+ Double_t xmaxLaser[12]= {336, //id
+ 2, //side
+ 6, //rod
+ 4, //bundle
+ 7, //beam
+ 1, //dP0
+ 1, //dP1
+ 0.01, //dP2
+ 0.01, //dP3
+ 0.1, //dP4
+ 160, //ncl
+ 40}; //sqrt dEdx
+
+ TString nameLaser[12]= {"id",
+ "side",
+ "rod",
+ "bundle",
+ "beam",
+ "dP0",
+ "dP1",
+ "dP2",
+ "dP3",
+ "dP4",
+ "ncl",
+ "sqrt dEdx"};
+ TString titleLaser[12]= {"id",
+ "side",
+ "rod",
+ "bundle",
+ "beam",
+ "#Delta_{P0}",
+ "#Delta_{P1}",
+ "#Delta_{P2}",
+ "#Delta_{P3}",
+ "#Delta_{P4}",
+ "N_{cl}",
+ "#sqrt{dEdx}"};
+ fHisLaser = new THnSparseS("dLaser","#Delta_{Laser}", 12, binsLaser,xminLaser, xmaxLaser);
+ for (Int_t iaxis=1; iaxis<12; iaxis++){
+ fHisLaser->GetAxis(iaxis)->SetName(nameLaser[iaxis]);
+ fHisLaser->GetAxis(iaxis)->SetTitle(titleLaser[iaxis]);
+ }
+ //
+ // Delta Time bin
+ // Pad SigmaShape Q charge pad row trackID
+ Int_t binsRow[6]={200, 10000, 20, 30, 159, 336};
+ Double_t axisMin[6]={-1, 0, 0, 1, 0 , 0};
+ Double_t axisMax[6]={ 1, 1000, 1, 30, 159, 336};
+ TString axisName[6]={"Delta","bin", "rms shape", "sqrt(Q)", "row","trackID"};
+
+ binsRow[1]=2000;
+ axisMin[1]=0;
+ axisMax[1]=200;
+ fHisLaserPad = new THnSparseS("laserPad","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
+ //
+ binsRow[0]=1000;
+ axisMin[0]=-20;
+ axisMax[0]=20;
+ binsRow[1]=10000;
+ axisMin[1]=0;
+ axisMax[1]=1000;
+ //
+ fHisLaserTime= new THnSparseS("laserTime","#Delta_{Laser}", 6, binsRow,axisMin, axisMax);
+ //
+ for (Int_t iaxis=0; iaxis<6; iaxis++){
+ fHisLaserPad->GetAxis(iaxis)->SetName(axisName[iaxis]);
+ fHisLaserTime->GetAxis(iaxis)->SetTitle(axisName[iaxis]);
+ }
+}
+
+void AliTPCcalibLaser::MergeFitHistos(AliTPCcalibLaser * laser){
+ //
+ // Merge content of histograms
+ //
+ // Only first histogram is checked - all other should be the same
+ if (fHisLaser &&laser->fHisLaser) fHisLaser->Add(laser->fHisLaser);
+ if (fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad->Add(laser->fHisLaserPad);
+ if (!fHisLaserPad &&laser->fHisLaserPad) fHisLaserPad=(THnSparseS*)laser->fHisLaserPad->Clone();
+ if (fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime->Add(laser->fHisLaserTime);
+ if (!fHisLaserTime &&laser->fHisLaserTime) fHisLaserTime=(THnSparseS*)laser->fHisLaserTime->Clone();
+
+ if (!laser->fHisNclIn) laser->MakeFitHistos(); // empty histograms
+ if (!fHisNclIn) MakeFitHistos();
+ if (fHisNclIn->GetEntries()<1) MakeFitHistos();
+ //
+ fHisNclIn->Add(laser->fHisNclIn ); //->Number of clusters inner
+ fHisNclOut->Add(laser->fHisNclOut ); //->Number of clusters outer
+ fHisNclIO->Add(laser->fHisNclIO ); //->Number of cluster inner outer
+ fHisLclIn->Add(laser->fHisLclIn ); //->Level arm inner
+ fHisLclOut->Add(laser->fHisLclOut ); //->Level arm outer
+ fHisLclIO->Add(laser->fHisLclIO ); //->Number of cluster inner outer
+ fHisdEdx->Add(laser->fHisdEdx ); //->dedx
+ fHisdZfit->Add(laser->fHisdZfit ); //->dedx
+ //
+ //
+ fHisChi2YIn1->Add(laser->fHisChi2YIn1 ); //->chi2 y inner - line
+ fHisChi2YOut1->Add(laser->fHisChi2YOut1 ); //->chi2 y inner - line
+ fHisChi2YIn2->Add(laser->fHisChi2YIn2 ); //->chi2 y inner - parabola
+ fHisChi2YOut2->Add(laser->fHisChi2YOut2 ); //->chi2 y inner - parabola
+ fHisChi2YIO1->Add(laser->fHisChi2YIO1 ); //->chi2 y IO - common
+ fHisChi2ZIn1->Add(laser->fHisChi2ZIn1 ); //->chi2 z inner - line
+ fHisChi2ZOut1->Add(laser->fHisChi2ZOut1 ); //->chi2 z inner - line
+ fHisChi2ZIn2->Add(laser->fHisChi2ZIn2 ); //->chi2 z inner - parabola
+ fHisChi2ZOut2->Add(laser->fHisChi2ZOut2 ); //->chi2 z inner - parabola
+ fHisChi2ZIO1->Add(laser->fHisChi2ZIO1 ); //->chi2 z IO - common
+ //
+ //
+ fHisPy1vP0->Add(laser->fHisPy1vP0 ); //-> delta y P0outer-P0inner - line
+ fHisPy2vP0->Add(laser->fHisPy2vP0 ); //-> delta y P0outer-P0inner - parabola
+ fHisPy3vP0->Add(laser->fHisPy3vP0 ); //-> delta y P0outer-P0inner - common parabola
+ fHisPy1vP1->Add(laser->fHisPy1vP1 ); //-> delta ky P1outer-P1inner - line
+ fHisPy2vP1->Add(laser->fHisPy2vP1 ); //-> delta ky P1outer-P1inner - parabola
+ fHisPy3vP1->Add(laser->fHisPy3vP1 ); //-> delta ky P1outer-P1inner - common parabola
+ fHisPy2vP2In->Add(laser-> fHisPy2vP2In ); //-> Curv P2inner - parabola
+ fHisPy2vP2Out->Add(laser->fHisPy2vP2Out ); //-> Curv P2outer - parabola
+ fHisPy3vP2IO->Add(laser->fHisPy3vP2IO ); //-> Curv P2outerinner - common parabola
+ //
+ //
+ fHisPz1vP0->Add(laser->fHisPz1vP0 ); //-> delta z P0outer-P0inner - line
+ fHisPz2vP0->Add(laser->fHisPz2vP0 ); //-> delta z P0outer-P0inner - parabola
+ fHisPz3vP0->Add(laser->fHisPz3vP0 ); //-> delta z P0outer-P0inner - common parabola
+ fHisPz1vP1->Add(laser->fHisPz1vP1 ); //-> delta kz P1outer-P1inner - line
+ fHisPz2vP1->Add(laser->fHisPz2vP1 ); //-> delta kz P1outer-P1inner - parabola
+ fHisPz3vP1->Add(laser->fHisPz3vP1 ); //-> delta kz P1outer-P1inner - common parabola
+ fHisPz2vP2In->Add(laser->fHisPz2vP2In ); //-> Curv P2inner - parabola
+ fHisPz2vP2Out->Add(laser->fHisPz2vP2Out ); //-> Curv P2outer - parabola
+ fHisPz3vP2IO->Add(laser->fHisPz3vP2IO ); //-> Curv P2outerinner - common parabola
+ fHisYAbsErrors->Add(laser->fHisYAbsErrors); //-> total errors per beam in the abs res y analysis
+ fHisZAbsErrors->Add(laser->fHisZAbsErrors); //-> total errors per beam in the abs res z analysis
+ //
+ //
+ //
+
+
+
+
+}
+
+
+
+
+void AliTPCcalibLaser::DumpFitInfo(TTree * chainFit,Int_t id){
+ //
+ // Dump fit information - collect information from the streamers
+ //
+ /*
+ TChain * chainFit=0;
+ TChain * chainTrack=0;
+ TChain * chain=0;
+ //
+ gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+ gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
+ AliXRDPROOFtoolkit tool;
+ chainTrack = tool.MakeChain("laser.txt","Track",0,10200);
+ chainTrack->Lookup();
+ chainTrack->SetProof(kTRUE);
+ chainDrift = tool.MakeChain("laser.txt","driftv",0,10200);
+ chainDrift->Lookup();
+ chainDrift->SetProof(kTRUE);
+
+ chain = tool.MakeChain("laser.txt","Residuals",0,10200);
+ chain->Lookup();
+ chainFit = tool.MakeChain("laser.txt","FitModels",0,10200);
+ chainFit->Lookup();
+ chainFit->SetProof(kTRUE);
+ chain->SetProof(kTRUE);
+ AliTPCLaserTrack::LoadTracks();
+ //AliTPCcalibLaser::DumpFitInfo(chainFit,0);
+
+ */
+ //
+ // Fit cuts
+ //
+ TCut cutP3("abs(Tr.fP[3])<0.1");
+ TCut cutP4("abs(Tr.fP[4])<0.5");
+ TCut cutPx = cutP3+cutP4;
+ TCut cutChi2YOut("sqrt(chi2y2Out*dEdx)<5&&chi2y2Out>0");
+ TCut cutChi2ZOut("sqrt(chi2z2Out*dEdx)<5&&chi2z2Out>0");
+ TCut cutChi2YIn("sqrt(chi2y2In*dEdx)<5&&chi2y2In>0");
+ TCut cutChi2ZIn("sqrt(chi2z2In*dEdx)<5&&chi2z2In>0");
+ //
+ TCut cutdEdx("sqrt(dEdx)>3");
+ TCut cutDY("abs(yPol2In.fElements[2]*nclI*nclI/4.)<3");
+ TCut cutN("nclO>20&&nclI>20");
+ TCut cutA = cutChi2YOut+cutChi2ZOut+cutChi2YIn+cutChi2ZIn+cutN+cutdEdx+cutPx+"accept";
+ //
+ // Cluster cuts
+ //
+ TCut cutClY("abs(Cl[].fY-TrYpol2.fElements)<0.15");
+ TCut cutClZ("abs(Cl[].fZ-TrZpol2.fElements)<0.15");
+ TCut cutClX("abs(Cl[].fX)>10");
+ TCut cutE("abs(Cl[].fY/Cl[].fX)<0.14");
+ TCut cutSY("sqrt(Cl[].fSigmaY2)>0.05");
+ TCut cutSZ("sqrt(Cl[].fSigmaZ2)>0.05");
+ TCut cutQ("sqrt(Cl[].fMax)>4");
+ TCut cutCl=cutClY+cutClZ+cutClX+cutE+cutSY+cutSZ+cutQ;
+
+
+ TH1F * phisAl = 0;
+ TH1F * phisAccept = 0;
+ TH1F * phisOut = 0;
+ TProfile * pdEdx = 0;
+
+ TProfile * pP0 = 0;
+ TProfile * pP1 = 0;
+ TProfile * pP2 = 0;
+ TProfile * pP3 = 0;
+ TProfile * pP4 = 0;
+ //
+ TProfile * pNclI = 0;
+ TProfile * pNclO = 0;
+ //
+ TProfile * pchi2YIn =0;
+ TProfile * pchi2ZIn =0;
+ TProfile * pchi2YOut =0;
+ TProfile * pchi2ZOut =0;
+ TProfile * pchi2YInOut =0;
+ TProfile * pchi2ZInOut =0;;
+ // laser counters
+ chainFit->Draw("LTr.fId>>hisAl(350,0,350)","LTr.fId<350");
+ phisAl = (TH1F*)gROOT->FindObject("hisAl");
+ chainFit->Draw("LTr.fId>>hisAccept(350,0,350)","LTr.fId<350"+cutA);
+ phisAccept = (TH1F*)gROOT->FindObject("hisAccept");
+ chainFit->Draw("LTr.fId>>hisOut(350,0,350)","LTr.fId<350"+!cutA);
+ phisOut = (TH1F*)gROOT->FindObject("hisOut");
+ //
+ chainFit->Draw("sqrt(dEdx):LTr.fId>>hdEdx(350,0,350)","","prof");
+ pdEdx = (TProfile*)gROOT->FindObject("hdEdx");
+ // track param
+ //
+ chainFit->Draw("Tr.fP[0]:LTr.fId>>hP0(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
+ pP0 = (TProfile*)gROOT->FindObject("hP0");
+ chainFit->Draw("Tr.fP[1]:LTr.fId>>hP1(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
+ pP1 = (TProfile*)gROOT->FindObject("hP1");
+ chainFit->Draw("Tr.fP[2]:LTr.fId>>hP2(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
+ pP2 = (TProfile*)gROOT->FindObject("hP2");
+ chainFit->Draw("Tr.fP[3]:LTr.fId>>hP3(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
+ pP3 = (TProfile*)gROOT->FindObject("hP3");
+ chainFit->Draw("Tr.fP[4]:LTr.fId>>hP4(350,0,350)","Tr.fP[4]/sqrt(Tr.fC[14])<20"+cutA,"prof");
+ pP4 = (TProfile*)gROOT->FindObject("hP4");
+ //
+ chainFit->Draw("nclI:LTr.fId>>hNclI(350,0,350)","","prof");
+ pNclI = (TProfile*)gROOT->FindObject("hNclI");
+ chainFit->Draw("nclO:LTr.fId>>hNclO(350,0,350)","","prof");
+ pNclO = (TProfile*)gROOT->FindObject("hNclO");
+ //
+ //
+ chainFit->Draw("sqrt(chi2y2In):LTr.fId>>hChi2YIn(350,0,350)",cutA+"","prof");
+ pchi2YIn = (TProfile*)gROOT->FindObject("hChi2YIn");
+ chainFit->Draw("sqrt(chi2y2Out):LTr.fId>>hChi2YOut(350,0,350)",cutA+"","prof");
+ pchi2YOut = (TProfile*)gROOT->FindObject("hChi2YOut");
+ chainFit->Draw("sqrt(chi2yInOut):LTr.fId>>hChi2YInOut(350,0,350)",cutA+"","prof");
+ pchi2YInOut = (TProfile*)gROOT->FindObject("hChi2YInOut");
+ chainFit->Draw("sqrt(chi2z2In):LTr.fId>>hChi2ZIn(350,0,350)",cutA+"","prof");
+ pchi2ZIn = (TProfile*)gROOT->FindObject("hChi2ZIn");
+ chainFit->Draw("sqrt(chi2z2Out):LTr.fId>>hChi2ZOut(350,0,350)",cutA+"","prof");
+ pchi2ZOut = (TProfile*)gROOT->FindObject("hChi2ZOut");
+ chainFit->Draw("sqrt(chi2zInOut):LTr.fId>>hChi2ZInOut(350,0,350)",cutA+"","prof");
+ pchi2ZInOut = (TProfile*)gROOT->FindObject("hChi2ZInOut");
+ //
+ // second derivatives
+ //
+ TH2F * phisPy2In = new TH2F("Py2Inner","Py2Inner",350,0,350,100,-0.001,0.001);
+ chainFit->Draw("yPol2In.fElements[2]:LTr.fId>>Py2Inner",cutA,"");
+ TH2F * phisPy2Out = new TH2F("Py2Outer","Py2Outer",350,0,350,200,-0.0005,0.0005);
+ chainFit->Draw("yPol2Out.fElements[2]:LTr.fId>>Py2Outer",cutA,"");
+ TH2F * phisPy2InOut = new TH2F("Py2InOut","Py2InOut",350,0,350,200,-0.0005,0.0005);
+ chainFit->Draw("yInOut.fElements[4]:LTr.fId>>Py2InOut",cutA,"");
+ //
+ phisPy2In->FitSlicesY();
+ TH1D * phisPy2InEntries = (TH1D*)gROOT->FindObject("Py2Inner_0");
+ TH1D * phisPy2InMean = (TH1D*)gROOT->FindObject("Py2Inner_1");
+ TH1D * phisPy2InSigma = (TH1D*)gROOT->FindObject("Py2Inner_2");
+ //
+ phisPy2Out->FitSlicesY();
+ TH1D * phisPy2OutEntries = (TH1D*)gROOT->FindObject("Py2Outer_0");
+ TH1D * phisPy2OutMean = (TH1D*)gROOT->FindObject("Py2Outer_1");
+ TH1D * phisPy2OutSigma = (TH1D*)gROOT->FindObject("Py2Outer_2");
+ //
+ phisPy2InOut->FitSlicesY();
+ TH1D * phisPy2InOutEntries = (TH1D*)gROOT->FindObject("Py2InOut_0");
+ TH1D * phisPy2InOutMean = (TH1D*)gROOT->FindObject("Py2InOut_1");
+ TH1D * phisPy2InOutSigma = (TH1D*)gROOT->FindObject("Py2InOut_2");
+ //
+ TH2F * phisPz2In = new TH2F("Pz2Inner","Pz2Inner",350,0,350,100,-0.001,0.001);
+ chainFit->Draw("zPol2In.fElements[2]:LTr.fId>>Pz2Inner",cutA,"");
+ TH2F * phisPz2Out = new TH2F("Pz2Outer","Pz2Outer",350,0,350,200,-0.0005,0.0005);
+ chainFit->Draw("zPol2Out.fElements[2]:LTr.fId>>Pz2Outer",cutA,"");
+ TH2F * phisPz2InOut = new TH2F("Pz2InOut","Pz2InOut",350,0,350,200,-0.0005,0.0005);
+ chainFit->Draw("zInOut.fElements[4]:LTr.fId>>Pz2InOut",cutA,"");
+ //
+ phisPz2In->FitSlicesY();
+ TH1D * phisPz2InEntries = (TH1D*)gROOT->FindObject("Pz2Inner_0");
+ TH1D * phisPz2InMean = (TH1D*)gROOT->FindObject("Pz2Inner_1");
+ TH1D * phisPz2InSigma = (TH1D*)gROOT->FindObject("Pz2Inner_2");
+ //
+ phisPz2Out->FitSlicesY();
+ TH1D * phisPz2OutEntries = (TH1D*)gROOT->FindObject("Pz2Outer_0");
+ TH1D * phisPz2OutMean = (TH1D*)gROOT->FindObject("Pz2Outer_1");
+ TH1D * phisPz2OutSigma = (TH1D*)gROOT->FindObject("Pz2Outer_2");
+ //
+ phisPz2InOut->FitSlicesY();
+ TH1D * phisPz2InOutEntries = (TH1D*)gROOT->FindObject("Pz2InOut_0");
+ TH1D * phisPz2InOutMean = (TH1D*)gROOT->FindObject("Pz2InOut_1");
+ TH1D * phisPz2InOutSigma = (TH1D*)gROOT->FindObject("Pz2InOut_2");
+ //
+ //
+ //
+
+
+ {
+ TTreeSRedirector *pcstream = new TTreeSRedirector("vscan.root");
+ for (Int_t ilaser=0; ilaser<336; ilaser++){
+ Float_t all=phisAl->GetBinContent(ilaser+1);
+ Float_t accept=phisAccept->GetBinContent(ilaser+1);
+ Float_t out=phisOut->GetBinContent(ilaser+1);
+ Float_t sdedx = pdEdx->GetBinContent(ilaser+1);
+ Float_t mP0 = pP0->GetBinContent(ilaser+1);
+ Float_t mP1 = pP1->GetBinContent(ilaser+1);
+ Float_t mP2 = pP2->GetBinContent(ilaser+1);
+ Float_t mP3 = pP3->GetBinContent(ilaser+1);
+ Float_t mP4 = pP4->GetBinContent(ilaser+1);
+
+
+ Float_t nclI = pNclI->GetBinContent(ilaser+1);
+ Float_t nclO = pNclO->GetBinContent(ilaser+1);
+ //
+ Float_t chi2YIn=pchi2YIn->GetBinContent(ilaser+1);
+ Float_t chi2YOut=pchi2YOut->GetBinContent(ilaser+1);
+ Float_t chi2YInOut=pchi2YInOut->GetBinContent(ilaser+1);
+ Float_t chi2ZIn=pchi2ZIn->GetBinContent(ilaser+1);
+ Float_t chi2ZOut=pchi2ZOut->GetBinContent(ilaser+1);
+ Float_t chi2ZInOut=pchi2ZInOut->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPy2In = phisPy2InEntries->GetBinContent(ilaser+1);
+ Float_t meanPy2In = phisPy2InMean->GetBinContent(ilaser+1);
+ Float_t sigmaPy2In = phisPy2InSigma->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPy2Out = phisPy2OutEntries->GetBinContent(ilaser+1);
+ Float_t meanPy2Out = phisPy2OutMean->GetBinContent(ilaser+1);
+ Float_t sigmaPy2Out = phisPy2OutSigma->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPy2InOut = phisPy2InOutEntries->GetBinContent(ilaser+1);
+ Float_t meanPy2InOut = phisPy2InOutMean->GetBinContent(ilaser+1);
+ Float_t sigmaPy2InOut = phisPy2InOutSigma->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPz2In = phisPz2InEntries->GetBinContent(ilaser+1);
+ Float_t meanPz2In = phisPz2InMean->GetBinContent(ilaser+1);
+ Float_t sigmaPz2In = phisPz2InSigma->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPz2Out = phisPz2OutEntries->GetBinContent(ilaser+1);
+ Float_t meanPz2Out = phisPz2OutMean->GetBinContent(ilaser+1);
+ Float_t sigmaPz2Out = phisPz2OutSigma->GetBinContent(ilaser+1);
+ //
+ Float_t entriesPz2InOut = phisPz2InOutEntries->GetBinContent(ilaser+1);
+ Float_t meanPz2InOut = phisPz2InOutMean->GetBinContent(ilaser+1);
+ Float_t sigmaPz2InOut = phisPz2InOutSigma->GetBinContent(ilaser+1);
+
+ AliTPCLaserTrack* ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(ilaser);
+ (*pcstream)<<"Scan"<<
+ "Run="<<id<<
+ "LTr.="<<ltrp<<
+ "all="<<all<<
+ "accept="<<accept<<
+ "out="<<out<<
+ "sdedx="<<sdedx<<
+ "mP0="<<mP0<<
+ "mP1="<<mP1<<
+ "mP2="<<mP2<<
+ "mP3="<<mP3<<
+ "mP4="<<mP4<<
+ "nclI="<<nclI<<
+ "nclO="<<nclO<<
+ "chi2YIn="<<chi2YIn<<
+ "chi2YOut="<<chi2YOut<<
+ "chi2YInOut="<<chi2YInOut<<
+ "chi2ZIn="<<chi2ZIn<<
+ "chi2ZOut="<<chi2ZOut<<
+ "chi2ZInOut="<<chi2ZInOut<<
+ //
+ "nPy2In="<<entriesPy2In<<
+ "mPy2In="<<meanPy2In<<
+ "sPy2In="<<sigmaPy2In<<
+ //
+ "nPy2Out="<<entriesPy2Out<<
+ "mPy2Out="<<meanPy2Out<<
+ "sPy2Out="<<sigmaPy2Out<<
+ //
+ "nPy2InOut="<<entriesPy2InOut<<
+ "mPy2InOut="<<meanPy2InOut<<
+ "sPy2InOut="<<sigmaPy2InOut<<
+ //
+ "nPz2In="<<entriesPz2In<<
+ "mPz2In="<<meanPz2In<<
+ "sPz2In="<<sigmaPz2In<<
+ //
+ "nPz2Out="<<entriesPz2Out<<
+ "mPz2Out="<<meanPz2Out<<
+ "sPz2Out="<<sigmaPz2Out<<
+ //
+ "nPz2InOut="<<entriesPz2InOut<<
+ "mPz2InOut="<<meanPz2InOut<<
+ "sPz2InOut="<<sigmaPz2InOut<<
+ "\n";
+ }
+
+ delete pcstream;
+ }
+}
+
+void AliTPCcalibLaser::SetBeamParameters(TVectorD& meanOffset,
+ TVectorD& meanSlope,
+ TVectorD& sectorArray,
+ Int_t option)
+{
+ // This method should ideally go in AliTPCLaser
+ // option == 0 (pads outer - closest to beam)
+ // option == 1 (pads inner)
+ // option == 2 (time outer)
+ // option == 3 (time inner)
+ Int_t nFailures = 0;
+
+ for(Int_t id = 0; id < 336; id++) {
+
+ if (!AliTPCLaserTrack::GetTracks())
+ AliTPCLaserTrack::LoadTracks();
+ AliTPCLaserTrack *ltrp =
+ (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+
+ AliExternalTrackParam trackParam(*ltrp);
+
+ Double_t deltaangle = 10.0; // sector is 1/2 a sector away from mirror
+ if((option==1 || option==3)&& (ltrp->GetBeam()<=1 || ltrp->GetBeam()>=5))
+ deltaangle = 30.0; // inner sector is 1 sector further away from mirror
+
+ Double_t angle = trackParam.GetAlpha();
+ if(angle<0)
+ angle += 2*TMath::Pi();
+ if(trackParam.GetSnp()>0) // track points to sector "before"
+ angle -= deltaangle*TMath::DegToRad();
+ else // track points to sector "after"
+ angle += deltaangle*TMath::DegToRad();
+
+ Bool_t success = trackParam.Rotate(angle);
+
+ if(!success) {
+ // cout << "WARNING: Rotate failed for ID: " << id << endl;
+ nFailures++;
+ }
+
+ angle *= TMath::RadToDeg();
+
+ Int_t sector = TMath::Nint((angle-10.0)/20.0);
+ if(sector<0)
+ sector += 18;
+ else if(sector>=18)
+ sector -= 18;
+ if(ltrp->GetSide()==1) // C side
+ sector += 18;
+ if(option==0 || option==2)
+ sector += 36;
+ if((option==1||option==3)&&(ltrp->GetBeam()==0||ltrp->GetBeam()==6))
+ sector += 36;
+
+ sectorArray[id] = sector;
+
+ const Double_t x0 = 0;
+
+ Double_t slopey = TMath::Tan(TMath::ASin(trackParam.GetSnp()));
+ Double_t slopez = trackParam.GetTgl();
+ // One needs a factor sqrt(1+slopey**2) to take into account the
+ // longer path length
+ slopez *= TMath::Sqrt(1.0 + slopey*slopey);
+ if(fInverseSlopeZ) // wrong sign in database, should be fixed there
+ slopez *= -1;
+ // Double_t offsetz = trackParam.GetZ();
+ Double_t offsety = trackParam.GetY() + slopey*(x0-trackParam.GetX());
+ Double_t offsetz = trackParam.GetZ() + slopez*(x0-trackParam.GetX());
+ if(option==2 || option==3) {
+ meanOffset[id] = offsetz; meanSlope[id] = slopez;
+ } else {
+ meanOffset[id] = offsety; meanSlope[id] = slopey;
+ }
+ }
+
+ if(nFailures>0)
+ AliWarning(Form("Rotate method failed %d times", nFailures));
+}
+
+
+
+void AliTPCcalibLaser::DumpLaser(const char *finput, Int_t run){
+ //
+ //
+ //input="TPCLaserObjects.root"
+ //
+ // 0. OBJ: TAxis Delta
+ // 1. OBJ: TAxis bin
+ // 2. OBJ: TAxis rms shape
+ // 3. OBJ: TAxis sqrt(Q)
+ // 4. OBJ: TAxis row
+ // 5. OBJ: TAxis trackID
+
+ const Double_t kSigma=4.;
+ TFile f(finput);
+ AliTPCcalibLaser *laserTPC = (AliTPCcalibLaser*) f.Get("laserTPC");
+ THnSparse * hisPadInput = laserTPC->fHisLaserPad;
+ THnSparse * hisTimeInput = laserTPC->fHisLaserTime;
+ TTreeSRedirector *pcstream= new TTreeSRedirector("hisLasers.root");
+ TVectorD meanY(159), sigmaY(159);
+ TVectorD meanZ(159), sigmaZ(159);
+ TVectorD meanPad(159), sigmaPad(159);
+ TVectorD meanTime(159), sigmaTime(159);
+ TVectorD meanDPad(159), sigmaDPad(159);
+ TVectorD meanDTime(159), sigmaDTime(159);
+ TVectorD meandEdx(159), sigmadEdx(159);
+ TVectorD meanSTime(159), sigmaSTime(159);
+ TVectorD meanSPad(159), sigmaSPad(159);
+ TVectorD entries(159);
+ //
+ Int_t indexes[10]={0,1,2,3,4,5,6};
+ TH1 *his=0;
+ AliTPCLaserTrack::LoadTracks();
+ //
+ for (Int_t id=0; id<336; id++){ // llop over laser beams
+ printf("id=\t%d\n",id);
+ //
+ AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ //
+ hisPadInput->GetAxis(5)->SetRange(id+1,id+1);
+ hisTimeInput->GetAxis(5)->SetRange(id+1,id+1);
+ //
+ his=hisTimeInput->Projection(3);
+ Int_t firstBindEdx=his->FindFirstBinAbove(0);
+ Int_t lastBindEdx=his->FindLastBinAbove(0);
+ hisPadInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
+ hisTimeInput->GetAxis(3)->SetRange(firstBindEdx, lastBindEdx);
+ delete his;
+ //
+ his=hisTimeInput->Projection(1);
+ // Int_t firstBinTime=his->FindFirstBinAbove(0);
+ //Int_t lastBinTime=his->FindLastBinAbove(0);
+ //hisTimeInput->GetAxis(1)->SetRange(firstBinTime, lastBinTime);
+ delete his;
+ //
+ //
+ his=hisTimeInput->Projection(2);
+ //Int_t firstBinZ=his->FindFirstBinAbove(0);
+ //Int_t lastBinZ=his->FindLastBinAbove(0);
+ //hisTimeInput->GetAxis(2)->SetRange(firstBinZ, lastBinZ);
+ delete his;
+ //
+ his=hisPadInput->Projection(2);
+ // Int_t firstBinY=his->FindFirstBinAbove(0);
+ //Int_t lastBinY=his->FindLastBinAbove(0);
+ //hisPadInput->GetAxis(2)->SetRange(firstBinY, lastBinY);
+ delete his;
+ //
+ //
+ //
+ THnSparse *hisPad0 = hisPadInput->Projection(5,indexes);
+ THnSparse *hisTime0 = hisTimeInput->Projection(5,indexes);
+ //
+ //
+ for (Int_t irow=0; irow<159; irow++){
+ entries[irow]=0;
+ if ((*(ltrp->GetVecSec()))[irow] <0) continue;
+ if ((*(ltrp->GetVecLX()))[irow] <80) continue;
+
+ hisPad0->GetAxis(4)->SetRange(irow+1,irow+1);
+ hisTime0->GetAxis(4)->SetRange(irow+1,irow+1);
+ //THnSparse *hisPad = hisPad0->Projection(4,indexes);
+ //THnSparse *hisTime = hisTime0->Projection(4,indexes);
+ THnSparse *hisPad = hisPad0;
+ THnSparse *hisTime = hisTime0;
+ //
+ // Get mean value of QA variables
+ //
+ // dEdx
+ his=hisTime->Projection(3);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meandEdx[irow] =his->GetMean();
+ sigmadEdx[irow]=his->GetRMS();
+ // Int_t bindedx0= his->FindBin(meandEdx[irow]-kSigma*sigmadEdx[irow]);
+ //Int_t bindedx1= his->FindBin(meandEdx[irow]+kSigma*sigmadEdx[irow]);
+ // hisPad->GetAxis(3)->SetRange(bindedx0,bindedx1);
+ //hisTime->GetAxis(3)->SetRange(bindedx0,bindedx1 );
+ delete his;
+ //
+ // sigma Time
+ //
+ his=hisTime->Projection(2);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()-kSigma*his->GetRMS());
+ meanSTime[irow] =his->GetMean();
+ sigmaSTime[irow]=his->GetRMS();
+ //Int_t binSTime0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
+ //Int_t binSTime1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
+ // hisTime->GetAxis(2)->SetRange(binSTime0, binSTime1);
+ delete his;
+ //
+ // sigma Pad
+ his=hisPad->Projection(2);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanSPad[irow] =his->GetMean();
+ sigmaSPad[irow]=his->GetRMS();
+ // Int_t binSPad0= his->FindBin(his->GetMean()-kSigma*his->GetRMS());
+ //Int_t binSPad1= his->FindBin(his->GetMean()+kSigma*his->GetRMS());
+ // hisPad->GetAxis(2)->SetRange(binSPad0, binSPad1);
+ delete his;
+ //
+ // apply selection on QA variables
+ //
+ //
+ //
+ // Y
+ his=hisPad->Projection(0);
+ entries[irow]=his->GetEntries();
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanY[irow] =his->GetMean();
+ sigmaY[irow]=his->GetRMS();
+ delete his;
+ // Z
+ his=hisTime->Projection(0);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanZ[irow] =his->GetMean();
+ sigmaZ[irow]=his->GetRMS();
+ delete his;
+ // Pad
+ his=hisPad->Projection(1);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanPad[irow] =his->GetMean();
+ meanDPad[irow] =his->GetMean()-Int_t(his->GetMean());
+ sigmaPad[irow]=his->GetRMS();
+ delete his;
+ // Time
+ his=hisTime->Projection(1);
+ his->GetXaxis()->SetRangeUser(his->GetMean()-kSigma*his->GetRMS(), his->GetMean()+kSigma*his->GetRMS());
+ meanTime[irow] = his->GetMean();
+ meanDTime[irow] = his->GetMean()-Int_t(his->GetMean());
+ sigmaTime[irow]=his->GetRMS();
+ delete his;
+ //
+ //delete hisTime;
+ //delete hisPad;
+ }
+ //
+ //
+ //
+ (*pcstream)<<"laserClusters"<<
+ "id="<<id<<
+ "run="<<run<<
+ "LTr.="<<ltrp<<
+ //
+ "entries.="<<&entries<<
+ "my.="<<&meanY<< //mean delta y
+ "rmsy.="<<&sigmaY<< //rms deltay
+ "mz.="<<&meanZ<< //mean deltaz
+ "rmsz.="<<&sigmaZ<< //rms z
+ //
+ "mPad.="<<&meanPad<< // mean pad
+ "mDPad.="<<&meanDPad<< // mead dpad
+ "rmsPad.="<<&sigmaPad<< // rms pad
+ "mTime.="<<&meanTime<<
+ "mDTime.="<<&meanTime<<
+ "rmsTime.="<<&sigmaTime<<
+ //
+ "mdEdx.="<<&meandEdx<< //mean dedx
+ "rmsdEdx.="<<&sigmadEdx<< //rms dedx
+ "mSPad.="<<&meanSPad<< //mean sigma pad
+ "rmsSPad.="<<&sigmaSPad<< //rms sigma pad
+ "mSTime.="<<&meanSTime<< //mean sigma time
+ "rmsSTime.="<<&sigmaSTime<<
+ "\n";
+ //
+ delete hisPad0;
+ delete hisTime0;
+ }
+ delete pcstream;
+
+ /*
+
+ */
+}
+
+void AliTPCcalibLaser::FitLaserClusters(Int_t run){
+ //
+ //
+ //input="TPCLaserObjects.root"
+ //Algorithm:
+ // 1. Select cluster candidates, remove outlyers
+ // edge clusters
+ // clusters with atypical spread (e.g due track overlaps)
+ // small amount of entries clusters (absolute minimal cut + raltive -to mean cut)
+ // 2. Fit the tracklets -per sector - in pad and time coordinate frame
+ // Remove outlyers
+ // Store info distance of track to pad, time center
+ // Fit the correction for distance to the center (sin,cos)
+ // 3. Do local fit
+ const Double_t kEpsilon=0.000001;
+ const Int_t kMinClusters=20;
+ const Double_t kEdgeCut=3;
+ const Double_t kDistCut=1.5; // cut distance to the ideal track
+ const Double_t kDistCutFit=0.5;
+ const Double_t kDistCutFitPad=0.25;
+ const Double_t kDistCutFitTime=0.25;
+ const Int_t kSmoothRow=5;
+ TFile f("hisLasers.root"); // Input file
+ TTree * treeInput=(TTree*)f.Get("laserClusters");
+ TTreeSRedirector *pcstream=new TTreeSRedirector("fitLasers.root");
+ TVectorD *vecN=0;
+ TVectorD *vecMY=0;
+ TVectorD *vecMZ=0;
+ TVectorD *vecPad=0;
+ TVectorD *vecTime=0;
+ TVectorD *vecSY=0;
+ TVectorD *vecSZ=0;
+ TVectorD *meandEdx=0;
+ TVectorD isOK(159);
+ TVectorD fitPad(159);
+ TVectorD fitTime(159);
+ TVectorD fitPadLocal(159);
+ TVectorD fitTimeLocal(159);
+ TVectorD fitDPad(159);
+ TVectorD fitDTime(159);
+ TVectorD fitIPad(159);
+ TVectorD fitITime(159);
+ Double_t chi2PadIROC=0;
+ Double_t chi2PadOROC=0;
+ //
+ treeInput->SetBranchAddress("my.",&vecMY);
+ treeInput->SetBranchAddress("mz.",&vecMZ);
+ treeInput->SetBranchAddress("mPad.",&vecPad);
+ treeInput->SetBranchAddress("mTime.",&vecTime);
+ treeInput->SetBranchAddress("rmsy.",&vecSY);
+ treeInput->SetBranchAddress("rmsz.",&vecSZ);
+ treeInput->SetBranchAddress("entries.",&vecN);
+ treeInput->SetBranchAddress("mdEdx.",&meandEdx);
+
+ AliTPCLaserTrack::LoadTracks();
+ //
+ //
+ TVectorD fitPadIROC(3), fitPadOROC(3);
+ TVectorD fitPadIROCSin(3), fitPadOROCSin(3);
+ TVectorD fitTimeIROC(3), fitTimeOROC(3);
+ TVectorD fitTimeIROCSin(3), fitTimeOROCSin(3);
+ //
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t refX=roc->GetPadRowRadii(0,roc->GetNRows(0)-1);
+
+ //
+ for (Int_t id=0; id<336; id++){
+ //
+ treeInput->GetEntry(id);
+ AliTPCLaserTrack *ltrp =(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ Int_t medianEntries = TMath::Nint(TMath::Median(159,vecN->GetMatrixArray()));
+ Double_t medianRMSY = TMath::Median(159,vecSY->GetMatrixArray());
+ Double_t rmsRMSY = TMath::RMS(159,vecSY->GetMatrixArray());
+ Double_t medianRMSZ = TMath::Median(159,vecSZ->GetMatrixArray());
+ Double_t rmsRMSZ = TMath::RMS(159,vecSZ->GetMatrixArray());
+ Double_t mdEdx = TMath::Median(159,meandEdx->GetMatrixArray());
+ Int_t sectorInner= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[63/2]);
+ Int_t sectorOuter= TMath::Nint(ltrp->GetVecSec()->GetMatrixArray()[64+96/2]);
+ TLinearFitter fitterY(2,"pol1");
+ TLinearFitter fitterZ(2,"pol1");
+ TLinearFitter fitterPad(2,"pol1");
+ TLinearFitter fitterTime(2,"pol1");
+ TLinearFitter fitterPadSin(2,"hyp1");
+ TLinearFitter fitterTimeSin(3,"hyp2");
+ //
+ //
+ for (UInt_t irow=0; irow<159; irow++){
+ fitPad[irow]=0; fitIPad[irow]=0; fitDPad[irow]=0;
+ fitTime[irow]=0; fitITime[irow]=0; fitDTime[irow]=0;
+ Double_t sign=(ltrp->GetZ()>0) ? 1.:-1.;
+ isOK[irow]=kFALSE;
+ fitPad[irow]=0;
+ fitTime[irow]=0;
+ Int_t sector=(irow<roc->GetNRows(0))? sectorInner:sectorOuter;
+ Int_t npads=(irow<roc->GetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0));
+ (*vecPad)[irow]-=npads*0.5;
+ //
+ if ((irow<roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if ((irow>=roc->GetNRows(0)) &&TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ //
+ if (TMath::Abs((*vecMY)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecMZ)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecPad)[irow])<kEpsilon) continue; //not determined position
+ if (TMath::Abs((*vecTime)[irow])<kEpsilon) continue; //not determined position
+ if ((*vecN)[irow]<0.5*medianEntries) continue; //small amount of clusters
+ if ((*vecSY)[irow]>medianRMSY+3*rmsRMSY) continue; //big sigma
+ if ((*vecSZ)[irow]>medianRMSZ+3*rmsRMSZ) continue; //big sigma
+ Double_t dEdge= TMath::Abs((*(ltrp->GetVecLY()))[irow])-(*(ltrp->GetVecLX()))[irow]*TMath::Tan(TMath::Pi()/18.); //edge cut
+ if (TMath::Abs(dEdge)<kEdgeCut) continue;
+ if (irow<roc->GetNRows(0)){
+ if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.4)>kDistCut) continue;
+ }
+ if (irow>roc->GetNRows(0)){
+ if (TMath::Abs(((*ltrp->GetVecLY())[irow])-sign*(*vecPad)[irow]*0.6)>kDistCut) continue;
+ }
+
+ isOK[irow]=kTRUE;
+ }
+ //
+ //fit OROC - get delta pad and delta time
+ //
+ fitterPad.ClearPoints();
+ fitterTime.ClearPoints();
+ fitterPadSin.ClearPoints();
+ fitterTimeSin.ClearPoints();
+ {for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ fitterPad.AddPoint(&x,y);
+ fitterTime.AddPoint(&x,z);
+ }}
+ chi2PadOROC=0;
+ if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
+ fitterPad.Eval();
+ fitterTime.Eval();
+ chi2PadOROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorOuter)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
+ fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
+ fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
+ fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
+ if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (isOK[irow]>0){
+ Double_t xxxPad[2]={TMath::Sin(2*TMath::Pi()*fitIPad[irow])};
+ Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
+ TMath::Cos(2*TMath::Pi()*fitITime[irow])};
+ fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
+ fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
+ }
+ }
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPadSin.FixParameter(0,0);
+ fitterTimeSin.FixParameter(0,0);
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ //
+ fitterPad.GetParameters(fitPadOROC);
+ fitterTime.GetParameters(fitTimeOROC);
+ fitterPadSin.GetParameters(fitPadOROCSin);
+ fitterTimeSin.GetParameters(fitTimeOROCSin);
+ }
+ //
+ //
+ //fit IROC
+ //
+ fitterPad.ClearPoints();
+ fitterTime.ClearPoints();
+ fitterPadSin.ClearPoints();
+ fitterTimeSin.ClearPoints();
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ fitterPad.AddPoint(&x,y);
+ fitterTime.AddPoint(&x,z);
+ }
+ chi2PadIROC=0;
+ if (fitterPad.GetNpoints()>kMinClusters&&fitterTime.GetNpoints()>kMinClusters){
+ fitterPad.Eval();
+ fitterTime.Eval();
+ chi2PadIROC=TMath::Sqrt(fitterPad.GetChisquare()/fitterPad.GetNpoints());
+ for (Int_t irow=2; irow<157; irow++){
+ if (isOK[irow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sectorInner)>0.1) continue;
+ if (TMath::Abs(ltrp->GetVecLX()->GetMatrixArray()[irow])<80) continue;
+ Double_t y=(*vecPad)[irow];
+ Double_t z=(*vecTime)[irow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[irow]-refX;
+ Double_t fitP=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ Double_t fitT=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitPad[irow]=fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x;
+ fitTime[irow]=fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x;
+ fitDPad[irow]=y-(fitterPad.GetParameter(0)+fitterPad.GetParameter(1)*x);
+ fitDTime[irow]=z-(fitterTime.GetParameter(0)+fitterTime.GetParameter(1)*x);
+ fitIPad[irow]=fitP-TMath::Nint(fitP-0.5);
+ fitITime[irow]=fitT-TMath::Nint(fitT-0.5);
+ if (fitDPad[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (fitDTime[irow]>kDistCutFit) isOK[irow]=kFALSE;
+ if (isOK[irow]>0.5){
+ Double_t xxxPad[3]={TMath::Sin(2*TMath::Pi()*fitIPad[irow]),
+ TMath::Cos(2*TMath::Pi()*fitIPad[irow])};
+ Double_t xxxTime[3]={TMath::Sin(2*TMath::Pi()*fitITime[irow]),
+ TMath::Cos(2*TMath::Pi()*fitITime[irow])};
+ fitterPadSin.AddPoint(xxxPad,fitDPad[irow]);
+ fitterTimeSin.AddPoint(xxxTime,fitDTime[irow]);
+ }
+ }
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPadSin.FixParameter(0,0);
+ fitterTimeSin.FixParameter(0,0);
+ fitterPadSin.Eval();
+ fitterTimeSin.Eval();
+ fitterPad.GetParameters(fitPadIROC);
+ fitterTime.GetParameters(fitTimeIROC);
+ fitterPadSin.GetParameters(fitPadIROCSin);
+ fitterTimeSin.GetParameters(fitTimeIROCSin);
+ }
+ for (Int_t irow=0; irow<159; irow++){
+ if (TMath::Abs(fitDPad[irow])<kEpsilon) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDTime[irow])<kEpsilon) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDPad[irow])>kDistCutFitPad) isOK[irow]=kFALSE;
+ if (TMath::Abs(fitDTime[irow])>kDistCutFitTime) isOK[irow]=kFALSE;
+ }
+ for (Int_t irow=kSmoothRow/2; irow<159-kSmoothRow/2; irow++){
+ fitPadLocal[irow]=0;
+ fitTimeLocal[irow]=0;
+ if (isOK[irow]<0.5) continue;
+ Int_t sector=(irow<Int_t(roc->GetNRows(0)))? sectorInner:sectorOuter;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[irow]-sector)>0.1) continue;
+ //
+ TLinearFitter fitterPadLocal(2,"pol1");
+ TLinearFitter fitterTimeLocal(2,"pol1");
+ Double_t xref=ltrp->GetVecLX()->GetMatrixArray()[irow];
+ for (Int_t delta=-kSmoothRow; delta<=kSmoothRow; delta++){
+ Int_t jrow=irow+delta;
+ if (jrow<0) jrow=0;
+ if (jrow>159) jrow=159;
+ if (isOK[jrow]<0.5) continue;
+ if (TMath::Abs(ltrp->GetVecSec()->GetMatrixArray()[jrow]-sector)>0.1) continue;
+ Double_t y=(*vecPad)[jrow];
+ Double_t z=(*vecTime)[jrow];
+ Double_t x=ltrp->GetVecLX()->GetMatrixArray()[jrow]-xref;
+ fitterPadLocal.AddPoint(&x,y);
+ fitterTimeLocal.AddPoint(&x,z);
+ }
+ if (fitterPadLocal.GetNpoints()<kSmoothRow) continue;
+ fitterPadLocal.Eval();
+ fitterTimeLocal.Eval();
+ fitPadLocal[irow]=fitterPadLocal.GetParameter(0);
+ fitTimeLocal[irow]=fitterTimeLocal.GetParameter(0);
+ }
+ //
+ //
+ (*pcstream)<<"fit"<<
+ "run="<<run<<
+ "id="<<id<<
+ "chi2PadIROC="<<chi2PadIROC<<
+ "chi2PadOROC="<<chi2PadOROC<<
+ "mdEdx="<<mdEdx<<
+ "LTr.="<<ltrp<<
+ "isOK.="<<&isOK<<
+ // mean measured-ideal values
+ "mY.="<<vecMY<<
+ "mZ.="<<vecMZ<<
+ // local coordinate fit
+ "mPad.="<<vecPad<<
+ "mTime.="<<vecTime<<
+ "fitPad.="<<&fitPad<<
+ "fitTime.="<<&fitTime<<
+ "fitPadLocal.="<<&fitPadLocal<<
+ "fitTimeLocal.="<<&fitTimeLocal<<
+ "fitDPad.="<<&fitDPad<<
+ "fitDTime.="<<&fitDTime<<
+ "fitIPad.="<<&fitIPad<<
+ "fitITime.="<<&fitITime<<
+ //
+ "fitPadIROC.="<<&fitPadIROC<< // pad fit linear IROC
+ "fitPadIROCSin.="<<&fitPadIROCSin<< // pad fit linear+ pad correction
+ "fitPadOROC.="<<&fitPadOROC<<
+ "fitPadOROCSin.="<<&fitPadOROCSin<<
+ //
+ "fitTimeIROC.="<<&fitTimeIROC<<
+ "fitTimeIROCSin.="<<&fitTimeIROCSin<<
+ "fitTimeOROC.="<<&fitTimeOROC<<
+ "fitTimeOROCSin.="<<&fitTimeOROCSin<<
+ "\n";
+ }
+ delete pcstream;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff