X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCcalibLaser.cxx;h=0eee70d827162ddaad62287df17da2b9b95dda02;hp=c9a51a2d7e7485889fe2e4802799e554b9dba70c;hb=d4e3890b93e7d6a22d4d1715c32a22dbe1bb2734;hpb=c03e32506516b3215faeb7b4aaec81c605f0b521 diff --git a/TPC/AliTPCcalibLaser.cxx b/TPC/AliTPCcalibLaser.cxx index c9a51a2d7e7..0eee70d8271 100644 --- a/TPC/AliTPCcalibLaser.cxx +++ b/TPC/AliTPCcalibLaser.cxx @@ -22,49 +22,66 @@ // // 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) + // (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("CalibObjects.root"); - TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); - AliTPCcalibLaser * laser = ( AliTPCcalibLaser *)array->FindObject("laserTPC"); - laser->DumpMeanInfo(-0.4) + 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>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); - TFile f("laserTPCDebug.root"); - TTree * treeT = (TTree*)f.Get("Track"); - - - // LASER scan + // + // + // Analyze LASER scan + // + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros"); gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+") - AliXRDPROOFtoolkit tool; + 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(); - treeT->Draw("(atan2(x1,x0)-atan2(lx1,lx0))*250.:fBundle","fSide==1&&fRod==0"+cutA,"prof") - - gSystem->Load("libSTAT.so") - TStatToolkit toolkit; - Double_t chi2; - TVectorD fitParam; - TMatrixD covMatrix; - Int_t npoints; - TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix); + 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(); */ @@ -78,6 +95,8 @@ #include "AliESDtrack.h" #include "AliTPCTracklet.h" #include "TH1D.h" +#include "TH1F.h" +#include "TProfile.h" #include "TVectorD.h" #include "TTreeStream.h" #include "TFile.h" @@ -86,13 +105,30 @@ #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 #include #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; @@ -101,62 +137,452 @@ ClassImp(AliTPCcalibLaser) AliTPCcalibLaser::AliTPCcalibLaser(): AliTPCcalibBase(), fESD(0), - fESDfriend(), + fESDfriend(0), + fNtracks(0), fTracksMirror(336), fTracksEsd(336), fTracksEsdParam(336), fTracksTPC(336), - fDeltaZ(336), // array of histograms of delta z for each track - fDeltaPhi(336), // array of histograms of delta z for each track - fDeltaPhiP(336), // array of histograms of delta z for each track - fSignals(336), // array of dedx signals - fFitAside(new TVectorD(3)), // drift fit - A side - fFitCside(new TVectorD(3)), // drift fit - C- side - fRun(0) + 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): +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 - fDeltaPhi(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 - fFitAside(new TVectorD(3)), // drift fit - A side - fFitCside(new TVectorD(3)), // drift fit - C- side - fRun(0) + // + // + 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); + 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; @@ -165,70 +591,101 @@ void AliTPCcalibLaser::Process(AliESDEvent * event) { if (!fESDfriend) { return; } + if (fESDfriend->TestSkipBit()) return; + if (fESD->GetNumberOfTracks()GetEventNumberInFile())); + // + // find CE background if present + // + if (AliTPCLaserTrack::GetTracks()==0) AliTPCLaserTrack::LoadTracks(); + TH1D hisCE("hhisCE","hhisCE",100,-100,100); + for (Int_t i=0;iGetNumberOfTracks();++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 run = fESD->GetRunNumber(); - fRun = run; + Int_t counter=0; for (Int_t i=0;iGetTrack(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(calibObject))) break; - if (track&&seed) FindMirror(track,seed); + 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 (counter0.3) continue; // tracks in saturation DumpLaser(id); - RefitLaser(id); - + if ( AcceptLaser(id) && fFitZ[id]<0.5){ + RefitLaserJW(id); + MakeDistHisto(id); + } } + } -void AliTPCcalibLaser::MakeDistHisto(){ +void AliTPCcalibLaser::MakeDistHisto(Int_t id){ // // // - for (Int_t id=0; id<336; id++){ + // for (Int_t id=0; id<336; id++){ // // - if (!fTracksEsdParam.At(id)) continue; - if (!AcceptLaser(id)) continue; + 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) UpdateFitHistos(); + 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); - - if (!hisdz){ - hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10); - hisdz->SetDirectory(0); - fDeltaZ.AddAt(hisdz,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),1000,0,1000); - hisSignal->SetDirectory(0); - fSignals.AddAt(hisSignal,id); - } + // AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id); AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id); @@ -249,141 +706,739 @@ void AliTPCcalibLaser::MakeDistHisto(){ 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(track->GetTPCsignal()); - } + 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::FitDriftV(){ // - // Fit drift velocity - linear approximation in the z and global y - // + // 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"); - fdriftA.ClearPoints(); - fdriftC.ClearPoints(); + 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++){ - if (!fTracksEsdParam.At(id)) continue; - if (!AcceptLaser(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 xxx[2] = {lxyz[2],lxyz[1]}; - if (ltrp->GetSide()==0){ - fdriftA.AddPoint(xxx,xyz[2],1); - }else{ - fdriftC.AddPoint(xxx,xyz[2],1); - } - } - Float_t chi2A = 0; - Float_t chi2C = 0; + 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; - // - if (fdriftA.GetNpoints()>10){ - fdriftA.Eval(); - fdriftA.EvalRobust(0.8); - fdriftA.GetParameters(*fFitAside); - npointsA= fdriftA.GetNpoints(); - chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints(); - } - if (fdriftC.GetNpoints()>10){ - fdriftC.Eval(); - fdriftC.EvalRobust(0.8); - fdriftC.GetParameters(*fFitCside); - npointsC= fdriftC.GetNpoints(); - chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints(); - } - - if (fStreamLevel>0){ - TTreeSRedirector *cstream = GetDebugStreamer(); - Int_t time = fESD->GetTimeStamp(); - if (cstream){ - (*cstream)<<"driftv"<< - "driftA.="<kSaturCut) continue; + if ( fClusterCounter[id]GetTPCsignal()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 (chi2AGetNoElements()<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 (chi2CGetNoElements()<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 (chi2ACGetXYZ(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; isensorNumSensors();isensor++){ + AliDCSSensor *gsensor = goofieArray->GetSensor(isensor); + if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp); + } + + if (cstream){ + (*cstream)<<"driftv"<< + "run="< + 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 -Bool_t AliTPCcalibLaser::AcceptLaser(Int_t id){ + 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 // - /* - TCut cutT("cutT","abs(Tr.fP[3])<0.06"); - TCut cutPt("cutPt","abs(Tr.fP[4])<0.1"); - TCut cutN("cutN","fTPCncls>70"); - TCut cutP("cutP","abs(atan2(x1,x0)-atan2(lx1,lx0))<0.03") - TCut cutA = cutT+cutPt+cutP; - */ - AliExternalTrackParam *param =(AliExternalTrackParam*)fTracksEsdParam.At(id); - AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id); - AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id); - - if (TMath::Abs(param->GetParameter()[4])>0.03) return kFALSE; - if (TMath::Abs(param->GetParameter()[3])>0.06) return kFALSE; - if (TMath::Abs(param->GetParameter()[2]-ltrp->GetParameter()[2])>0.06) return kFALSE; - if (TMath::Abs(param->GetParameter()[1]-ltrp->GetParameter()[1])>10) return kFALSE; + 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; + } // - // dedx cut // - if (TMath::Abs(track->GetTPCsignal())<20) return kFALSE; - if (TMath::Abs(track->GetTPCsignal())>800) return kFALSE; + 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 // - return kTRUE; -} - -Int_t AliTPCcalibLaser::FindMirror(AliESDtrack *track, AliTPCseed *seed){ // - // Find corresponding mirror - // add the corresponding tracks + 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]0) continue; + if (ltrp->GetSide()==0){ + npointsA++; + nbA[ltrp->GetBundle()]++; + } + if (ltrp->GetSide()==1){ + npointsC++; + nbC[ltrp->GetBundle()]++; + } + } // - Float_t kRadius0 = 252; - Float_t kRadius = 253.4; - if (!track->GetOuterParam()) return -1; - AliExternalTrackParam param(*(track->GetOuterParam())); - AliTracker::PropagateTrackTo(¶m,kRadius0,0.10566,3,kTRUE); - AliTracker::PropagateTrackTo(¶m,kRadius,0.10566,0.1,kTRUE); - AliTPCLaserTrack ltr; - AliTPCLaserTrack *ltrp=0x0; + // reject "bad events" // - Int_t id = AliTPCLaserTrack::IdentifyTrack(¶m); - if (id!=-1 && (AliTPCLaserTrack::GetTracks()->UncheckedAt(id))) - ltrp=(AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id); - else - ltrp=<r; + Bool_t isOK=kTRUE; + Int_t naA=0; + Int_t naC=0; + if (npointsAminFraction*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; - if (id>=0){ - // - // - Float_t radius=TMath::Abs(ltrp->GetX()); - AliTracker::PropagateTrackTo(¶m,radius,0.10566,0.01,kTRUE); + // + // + // + for (Int_t iter=0; iter<2; iter++){ + fdriftA.ClearPoints(); + fdriftC.ClearPoints(); + fdriftAC.ClearPoints(); + npointsA=0; npointsC=0; npointsAC=0; // - if (!fTracksMirror.At(id)) fTracksMirror.AddAt(ltrp,id); - fTracksEsdParam.AddAt(param.Clone(),id); - fTracksEsd.AddAt(track,id); - fTracksTPC.AddAt(seed,id); + 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]GetTPCsignal()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 (chi2AGetNoElements()<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 (chi2CGetNoElements()<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 (chi2ACGetXYZ(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; isensorNumSensors();isensor++){ + AliDCSSensor *gsensor = goofieArray->GetSensor(isensor); + if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp); + } + + if (cstream){ + (*cstream)<<"driftvN"<< + "run="<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){ + // + // 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,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; + // 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 (dist0GetPxPyPz(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(); Int_t time = fESD->GetTimeStamp(); Bool_t accept = AcceptLaser(id); if (cstream){ (*cstream)<<"Track"<< - "run="<GetClusterPointer(i); - if (c) nclusters[c->GetDetector()]++; - } - - for (Int_t isec=0; isec<72;isec++){ - if (nclusters[isec]GetClusterPointer(irow); - //if (c && RejectCluster(c)) continue; - if (c&&c->GetDetector()==isec) { - Double_t xd = c->GetX()-120;; - Double_t x[2]={xd,xd*xd}; - fy2.AddPoint(x,c->GetY()); - fz2.AddPoint(x,c->GetZ()); + + 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; icutGetInner() || !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(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(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()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()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="<5){ + TTreeSRedirector *cstream = GetDebugStreamer(); + if (cstream){ + Float_t dedx = track->GetdEdx(); + (*cstream)<<"Residuals"<< + "run="<-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 // - fy1.AddPoint(x,c->GetY()); - fz1.AddPoint(x,c->GetZ()); + // + 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 } - } - fy2.Eval(); - fz2.Eval(); - fy1.Eval(); - fz1.Eval(); - fy1.GetParameters(vecy1); - fy2.GetParameters(vecy2); - fz1.GetParameters(vecz1); - fz2.GetParameters(vecz2); - - if (fStreamLevel>0){ - TTreeSRedirector *cstream = GetDebugStreamer(); - if (cstream){ - Float_t dedx = track->GetdEdx(); - (*cstream)<<"Tracklet"<< - "LTr.="<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<160;++irow) { - // AliTPCclusterMI *c=track->GetClusterPointer(irow); - // if (c && RejectCluster(c)) continue; - // if (c&&c->GetDetector()==isec) { - // Double_t x[2]={c->GetX(),c->GetX()*c->GetX()}; - // fy2.AddPoint(&x,c->GetY()); - // fz2.AddPoint(&x,c->GetZ()); - // // - // fy1.AddPoint(&x,c->GetY()); - // fz1.AddPoint(&x,c->GetZ()); - // } - // } - + 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(Float_t bfield){ - // + +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 (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 * hisS = laser->fDeltaZ->At(id); - if (!hisphi) continue;; - Double_t entries = hisphi->GetEntries(); - if (entries<30) continue; + 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 (entriesUncheckedAt(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(); - 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); - 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); - hisZ->Fit(&fg,"","",hisZ->GetMean()-4*hisZ->GetRMS(),hisZ->GetMean()+4*hisZ->GetRMS()); - Double_t gz1 = fg.GetParameter(1); - Double_t gz2 = fg.GetParameter(2); - + 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); - // - (*pcstream)<<"Mean"<< - "entries="<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; -void AliTPCcalibLaser::Analyze(){ - // - // - // -} + 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]; -Long64_t AliTPCcalibLaser::Merge(TCollection *li) { + 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="<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; ientryDraw("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; ientryMakeIterator(); 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()); + // Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName()); return -1; } + AliDebug(5,Form("Marging number %d\n", counter0)); + counter0++; // - // fHistNTracks->Add(cal->fHistNTracks); -// fClusters->Add(cal->fClusters); -// fModules->Add(cal->fModules); -// fHistPt->Add(cal->fHistPt); -// fPtResolution->Add(cal->fPtResolution); -// fDeDx->Add(cal->fDeDx); - - + 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); + hm = (TH1F*)cal->fDeltaZ.At(id); h = (TH1F*)fDeltaZ.At(id); - if (!h) { - h=new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10); - fDeltaZ.AddAt(h,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 (hm) h->Add(hm); + if (h&&hm) h->Add(hm); + + // // merge fDeltaPhi histograms - hm = (TH1F*)cal->fDeltaPhi.At(id); + hm = (TH1F*)cal->fDeltaPhi.At(id); h = (TH1F*)fDeltaPhi.At(id); - if (!h) { - h= new TH1F(Form("hisdphi%d",id),Form("hisdphi%d",id),1000,-1,1); + if (!h &&hm &&hm->GetEntries()>0) { + h= (TH1F*)hm->Clone(); + h->SetDirectory(0); fDeltaPhi.AddAt(h,id); } - if (hm) h->Add(hm); + if (h&&hm) h->Add(hm); // merge fDeltaPhiP histograms - hm = (TH1F*)cal->fDeltaPhiP.At(id); + hm = (TH1F*)cal->fDeltaPhiP.At(id); h = (TH1F*)fDeltaPhiP.At(id); - if (!h) { - h=new TH1F(Form("hisdphiP%d",id),Form("hisdphiP%d",id),1000,-0.01,0.01); + if (!h&&hm &&hm->GetEntries()>0) { + h=(TH1F*)hm->Clone(); + h->SetDirectory(0); fDeltaPhiP.AddAt(h,id); } - if (hm) h->Add(hm); + if (h&&hm) h->Add(hm); // merge fSignals histograms - hm = (TH1F*)cal->fSignals.At(id); + hm = (TH1F*)cal->fSignals.At(id); h = (TH1F*)fSignals.At(id); - if (!h) { - h=new TH1F(Form("hisSignal%d",id),Form("hisSignal%d",id),1000,0,1000); + if (!h&&hm &&hm->GetEntries()>0) { + h=(TH1F*)hm->Clone(); + h->SetDirectory(0); fSignals.AddAt(h,id); } - if (hm) h->Add(hm); + 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 -/* - gSystem->Load("libSTAT.so") - TStatToolkit toolkit; - Double_t chi2; - TVectorD fitParam; - TMatrixD covMatrix; - Int_t npoints; - TCut cutA("gphi2<0.2&&abs(mphi-gphi1)<0.1") + 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 -// TString *strq0 = toolkit.FitPlane(treeT,"Tr.fP[1]-LTr.fP[1]","lx1++lx2", "fSide==1"+cutA, chi2,npoints,fitParam,covMatrix); + + // + // + // + 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); -TString fstring=""; -fstring+="LTr.fP[2]++" // 1 -fstring+="LTr.fP[2]*cos(atan2(lx1,lx0))++" // 2 -fstring+="LTr.fP[2]*sin(atan2(lx1,lx0))++" // 3 -fstring+="cos(atan2(lx1,lx0))++" // 4 -fstring+="sin(atan2(lx1,lx0))++" // 5 -fstring+="gphiP1++" // 6 + if (!hisdz){ + hisdz = new TH1F(Form("hisdz%d",id),Form("hisdz%d",id),1000,-10,10); + hisdz->SetDirectory(0); + fDeltaZ.AddAt(hisdz,id); -fstring+="bz++"; // 7 -fstring+="bz*cos(atan2(lx1,lx0))++" // 8 -fstring+="bz*sin(atan2(lx1,lx0))++" // 9 -fstring+="bz*sin(atan2(lx1,lx0))*LTr.fP[2]++" // 10 -fstring+="bz*cos(atan2(lx1,lx0))*LTr.fP[2]++" // 11 + 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); + } + } - TString *strq0 = toolkit.FitPlane(chain,"gphi1",fstring->Data(), "fSide==1&&fBundle==3"+cutA, chi2,npoints,fitParam,covMatrix); + // + // 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::UpdateFitHistos(){ + //create the fit histos and set the beam parameters(needs OCDB access) + MakeFitHistos(); + SetBeamParameters(fBeamOffsetZOuter, fBeamSlopeZOuter, fBeamSectorOuter,2); + SetBeamParameters(fBeamOffsetZInner, fBeamSlopeZInner, fBeamSectorInner,3); + SetBeamParameters(fBeamOffsetYOuter, fBeamSlopeYOuter, fBeamSectorOuter,0); + SetBeamParameters(fBeamOffsetYInner, fBeamSlopeYInner, fBeamSectorInner,1); +} + +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="<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 - chain->SetAlias("fit",strq0->Data()); + 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="<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=(irowGetNRows(0))? sectorInner:sectorOuter; + Int_t npads=(irowGetNRows(0))? roc->GetNPads(sector,irow):roc->GetNPads(sector,irow-roc->GetNRows(0)); + (*vecPad)[irow]-=npads*0.5; + // + if ((irowGetNRows(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])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)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])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=(irowGetNRows(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()