From 51cfc50dcd5101f7906a9e32c961e27b11e7a1dd Mon Sep 17 00:00:00 2001 From: dainese Date: Mon, 10 May 2010 22:19:21 +0000 Subject: [PATCH] Update of ITS-TPC alignment code (Mikolaj) --- PWG1/AliAnalysisTaskITSTPCalignment.cxx | 680 +++++++++++++++--- PWG1/AliAnalysisTaskITSTPCalignment.h | 64 +- PWG1/AliRelAlignerKalmanArray.cxx | 462 ++++++++---- PWG1/AliRelAlignerKalmanArray.h | 67 +- PWG1/macros/AddTaskITSTPCalignment.C | 80 +++ PWG1/macros/plotAnalysisTaskITSTPCalignment.C | 674 +++++++++++++++++ PWG1/macros/runAnalysisTaskITSTPCalignment.C | 294 ++++++++ 7 files changed, 2022 insertions(+), 299 deletions(-) create mode 100644 PWG1/macros/AddTaskITSTPCalignment.C create mode 100644 PWG1/macros/plotAnalysisTaskITSTPCalignment.C create mode 100644 PWG1/macros/runAnalysisTaskITSTPCalignment.C diff --git a/PWG1/AliAnalysisTaskITSTPCalignment.cxx b/PWG1/AliAnalysisTaskITSTPCalignment.cxx index 7119a06914c..a77b9e3e298 100644 --- a/PWG1/AliAnalysisTaskITSTPCalignment.cxx +++ b/PWG1/AliAnalysisTaskITSTPCalignment.cxx @@ -4,16 +4,30 @@ // Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch //////////////////////////////////////////////////////////////////////////////// -#include "TChain.h" -#include "TTree.h" -#include "TH2.h" -#include "AliAnalysisTask.h" -#include "AliAnalysisManager.h" -#include "AliESDEvent.h" -#include "AliESDInputHandler.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include "AliRelAlignerKalman.h" #include "AliRelAlignerKalmanArray.h" #include "AliAnalysisTaskITSTPCalignment.h" +#include "AliLog.h" ClassImp(AliAnalysisTaskITSTPCalignment) @@ -21,36 +35,74 @@ ClassImp(AliAnalysisTaskITSTPCalignment) AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment(): AliAnalysisTask(), fESD(0), - fArray(0), - fYZResidualsHist(0), - fPLResidualsHist(0), - fListOfHistos(0), - fSaveInterval(600), - fTimeMatchingTolerance(20), - fDoQA(kFALSE) + fESDfriend(0), + fMC(0), + fArrayITSglobal(0), + fArrayITSsa(0), + fDebugTree(0), + fAligner(0), + fList(0), + fFillDebugTree(kFALSE), + fDoQA(kTRUE), + fT0(0), + fTend(0), + fSlotWidth(0), + fMinPt(0.4), + fMinPointsVol1(3), + fMinPointsVol2(80), + fRejectOutliers(kTRUE), + fOutRejSigma(1.), + fRejectOutliersSigma2Median(kFALSE), + fOutRejSigma2Median(5.), + fOutRejSigmaOnMerge(10.), + fUseITSoutGlobalTrack(kFALSE), + fUseITSoutITSSAtrack(kTRUE) { //dummy ctor + // Define input and output slots here + // Input slot #0 works with a TChain + //DefineInput(0, TChain::Class()); + //DefineOutput(0, TTree::Class()); + //DefineOutput(1, TList::Class()); + //DefineOutput(2, AliRelAlignerKalmanArray::Class()); } //________________________________________________________________________ AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment(const char *name): AliAnalysisTask(name,name), fESD(0), - fArray(0), - fYZResidualsHist(0), - fPLResidualsHist(0), - fListOfHistos(0), - fSaveInterval(600), - fTimeMatchingTolerance(20), - fDoQA(kFALSE) + fESDfriend(0), + fMC(0), + fArrayITSglobal(0), + fArrayITSsa(0), + fDebugTree(0), + fAligner(0), + fList(0), + fFillDebugTree(kFALSE), + fDoQA(kTRUE), + fT0(0), + fTend(0), + fSlotWidth(0), + fMinPt(0.4), + fMinPointsVol1(3), + fMinPointsVol2(80), + fRejectOutliers(kTRUE), + fOutRejSigma(1.), + fRejectOutliersSigma2Median(kFALSE), + fOutRejSigma2Median(5.), + fOutRejSigmaOnMerge(10.), + fUseITSoutGlobalTrack(kFALSE), + fUseITSoutITSSAtrack(kTRUE) { // Constructor // Define input and output slots here // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); - DefineOutput(0, AliRelAlignerKalmanArray::Class()); + DefineOutput(0, TTree::Class()); DefineOutput(1, TList::Class()); + DefineOutput(2, AliRelAlignerKalmanArray::Class()); + DefineOutput(3, AliRelAlignerKalmanArray::Class()); } //________________________________________________________________________ @@ -60,18 +112,53 @@ void AliAnalysisTaskITSTPCalignment::ConnectInputData(Option_t *) TTree* tree = dynamic_cast (GetInputData(0)); if (!tree) { - Printf("ERROR: Could not read chain from input slot 0"); + AliError("Could not read chain from input slot 0"); + return; + } + + AliESDInputHandler *esdH = dynamic_cast + (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + if (!esdH) + { + AliError("Could not get ESDInputHandler"); + return; } else + fESD = esdH->GetEvent(); + + AliMCEventHandler* mcH = dynamic_cast + (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); + if (!mcH) { + AliWarning("Could not retrieve MC event handler"); + } + else + fMC = mcH->MCEvent(); + + //esdH->SetFriendFileName("AliESDfriends.root"); + //fESDfriend = esdH->GetESDfriend(); + //attach the ESD friend + tree->SetBranchStatus("ESDfriend*", 1); + fESDfriend = (AliESDfriend*)fESD->FindListObject("AliESDfriend"); + if (!fESDfriend) { - AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); - if (!esdH) - { - Printf("ERROR: Could not get ESDInputHandler"); - } - else - fESD = esdH->GetEvent(); + // works for both, we just want to avoid setting the branch adress twice + // in case of the new ESD + tree->SetBranchAddress("ESDfriend.",&fESDfriend); + } +} + +//________________________________________________________________________ +Bool_t AliAnalysisTaskITSTPCalignment::Notify() +{ + //ON INPUT FILE/TREE CHANGE + + //fArray->Print(); + if (fFillDebugTree) + { + if (fAligner->GetNUpdates()>0) fDebugTree->Fill(); } + + return kTRUE; } //________________________________________________________________________ @@ -79,17 +166,241 @@ void AliAnalysisTaskITSTPCalignment::CreateOutputObjects() { // Create output objects // Called once + fArrayITSglobal = new AliRelAlignerKalmanArray(); + fArrayITSglobal->SetName("outputArrayITSglobal"); + fArrayITSglobal->SetupArray(fT0,fTend,fSlotWidth); + fArrayITSglobal->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge); + //set up the template + AliRelAlignerKalman* templ = fArrayITSglobal->GetAlignerTemplate(); + templ->SetRejectOutliers(fRejectOutliers); + templ->SetOutRejSigma(fOutRejSigma); + templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median); + templ->SetOutRejSigma2Median(fOutRejSigma2Median); + fArrayITSsa = new AliRelAlignerKalmanArray(); + fArrayITSsa->SetName("outputArrayITSsa"); + fArrayITSsa->SetupArray(fT0,fTend,fSlotWidth); + fArrayITSsa->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge); + //set up the template + templ = fArrayITSsa->GetAlignerTemplate(); + templ->SetRejectOutliers(fRejectOutliers); + templ->SetOutRejSigma(fOutRejSigma); + templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median); + templ->SetOutRejSigma2Median(fOutRejSigma2Median); - fArray = new AliRelAlignerKalmanArray("ITSTPCalignmentArray"); - fArray->SetSaveInterval(fSaveInterval); - fArray->SetTimeMatchingTolerance(fTimeMatchingTolerance); + fList = new TList(); - fYZResidualsHist = new TH2F("fYZResidualsHist", "YZ residuals", 50, -0.5, 0.5, 50, -2., 2. ); - fPLResidualsHist = new TH2F("fPLResidualsHist", "sin(phi) tan(lambda) residuals", 50, -.05, 0.05, 50, -0.05, 0.05 ); + TH2F* pZYAResidualsHistBpos = new TH2F("fZYAResidualsHistBpos","z-r\\phi residuals side A (z>0), Bpos", 100, -4., 4., 100, -2.0, 2.0 ); + pZYAResidualsHistBpos->SetXTitle("\\deltaz [cm]"); + pZYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pZYCResidualsHistBpos = new TH2F("fZYCResidualsHistBpos","z-r\\phi residuals side C (z<0), Bpos", 100, -4., 4., 100, -2.0, 2.0 ); + pZYCResidualsHistBpos->SetXTitle("\\deltaz [cm]"); + pZYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLPAResidualsHistBpos = new TH2F("fLPAResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bpos", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPAResidualsHistBpos->SetXTitle("\\deltasin(\\phi)"); + pLPAResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pLPCResidualsHistBpos = new TH2F("fLPCResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bpos", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPCResidualsHistBpos->SetXTitle("\\deltasin(\\phi)"); + pLPCResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pPhiYAResidualsHistBpos = new TH2F("fPhiYAResidualsHistBpos","\\phi-\\deltar\\phi side A (z>0), Bpos",36,0.0,6.3,100,-2.0,2.0); + pPhiYAResidualsHistBpos->SetXTitle("\\phi [rad]"); + pPhiYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiYCResidualsHistBpos = new TH2F("fPhiYCResidualsHistBpos","\\phi-\\deltar\\phi side C (z<0), Bpos",36,0.0,6.3,100,-2.0,2.0); + pPhiYCResidualsHistBpos->SetXTitle("\\phi [rad]"); + pPhiYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiZAResidualsHistBpos = new TH2F("fPhiZAResidualsHistBpos","\\phi-\\deltaz side A (z>0), Bpos",36,0.0,6.3,100,-4.0,4.0); + pPhiZAResidualsHistBpos->SetXTitle("\\phi [rad]"); + pPhiZAResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TH2F* pPhiZCResidualsHistBpos = new TH2F("fPhiZCResidualsHistBpos","\\phi-\\deltaz side C (z<0), Bpos",36,0.0,6.3,100,-4.0,4.0); + pPhiZCResidualsHistBpos->SetXTitle("\\phi [rad]"); + pPhiZCResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TH2F* pPtYAResidualsHistBpos = new TH2F("fPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",20,.3,8.,80,-2.0,2.0); + pPtYAResidualsHistBpos->SetXTitle("Pt [rad]"); + pPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtYCResidualsHistBpos = new TH2F("fPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",20,.3,8.,80,-2.0,2.0); + pPtYCResidualsHistBpos->SetXTitle("Pt [rad]"); + pPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtZAResidualsHistBpos = new TH2F("fPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",20,.3,8.,80,-4.0,4.0); + pPtZAResidualsHistBpos->SetXTitle("Pt"); + pPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TH2F* pPtZCResidualsHistBpos = new TH2F("fPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",20,.3,8.,80,-4.0,4.0); + pPtZCResidualsHistBpos->SetXTitle("Pt"); + pPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtYAResidualsHistBpos = new TH2F("fLowPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",100,0.0,1.5,80,-2.0,2.0); + pLowPtYAResidualsHistBpos->SetXTitle("Pt [rad]"); + pLowPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtYCResidualsHistBpos = new TH2F("fLowPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",100,0.0,1.5,80,-2.0,2.0); + pLowPtYCResidualsHistBpos->SetXTitle("Pt [rad]"); + pLowPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtZAResidualsHistBpos = new TH2F("fLowPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",100,0.0,1.5,80,-4.0,4.0); + pLowPtZAResidualsHistBpos->SetXTitle("Pt"); + pLowPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtZCResidualsHistBpos = new TH2F("fLowPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",100,0.0,1.5,80,-4.0,4.0); + pLowPtZCResidualsHistBpos->SetXTitle("Pt"); + pLowPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]"); + TList* listBpos = new TList(); + fList->Add(listBpos); + listBpos->Add(pZYAResidualsHistBpos); + listBpos->Add(pZYCResidualsHistBpos); + listBpos->Add(pLPAResidualsHistBpos); + listBpos->Add(pLPCResidualsHistBpos); + listBpos->Add(pPhiYAResidualsHistBpos); + listBpos->Add(pPhiYCResidualsHistBpos); + listBpos->Add(pPhiZAResidualsHistBpos); + listBpos->Add(pPhiZCResidualsHistBpos); + listBpos->Add(pPtYAResidualsHistBpos); + listBpos->Add(pPtYCResidualsHistBpos); + listBpos->Add(pPtZAResidualsHistBpos); + listBpos->Add(pPtZCResidualsHistBpos); + listBpos->Add(pLowPtYAResidualsHistBpos); + listBpos->Add(pLowPtYCResidualsHistBpos); + listBpos->Add(pLowPtZAResidualsHistBpos); + listBpos->Add(pLowPtZCResidualsHistBpos); - fListOfHistos = new TList(); - fListOfHistos->Add(fYZResidualsHist); - fListOfHistos->Add(fPLResidualsHist); + TH2F* pZYAResidualsHistBneg = new TH2F("fZYAResidualsHistBneg","z-r\\phi residuals side A (z>0), Bneg", 100, -4., 4., 100, -2.0, 2.0 ); + pZYAResidualsHistBneg->SetXTitle("\\deltaz [cm]"); + pZYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pZYCResidualsHistBneg = new TH2F("fZYCResidualsHistBneg","z-r\\phi residuals side C (z<0), Bneg", 100, -4., 4., 100, -2.0, 2.0 ); + pZYCResidualsHistBneg->SetXTitle("\\deltaz [cm]"); + pZYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLPAResidualsHistBneg = new TH2F("fLPAResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bneg", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPAResidualsHistBneg->SetXTitle("\\deltasin(\\phi)"); + pLPAResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pLPCResidualsHistBneg = new TH2F("fLPCResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bneg", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPCResidualsHistBneg->SetXTitle("\\deltasin(\\phi)"); + pLPCResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pPhiYAResidualsHistBneg = new TH2F("fPhiYAResidualsHistBneg","\\phi-\\deltar\\phi side A (z>0), Bneg",36,0.0,6.3,100,-2.0,2.0); + pPhiYAResidualsHistBneg->SetXTitle("\\phi [rad]"); + pPhiYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiYCResidualsHistBneg = new TH2F("fPhiYCResidualsHistBneg","\\phi-\\deltar\\phi side C (z<0), Bneg",36,0.0,6.3,100,-2.0,2.0); + pPhiYCResidualsHistBneg->SetXTitle("\\phi [rad]"); + pPhiYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiZAResidualsHistBneg = new TH2F("fPhiZAResidualsHistBneg","\\phi-\\deltaz side A (z>0), Bneg",36,0.0,6.3,100,-4.0,4.0); + pPhiZAResidualsHistBneg->SetXTitle("\\phi [rad]"); + pPhiZAResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TH2F* pPhiZCResidualsHistBneg = new TH2F("fPhiZCResidualsHistBneg","\\Pt-\\deltaz side C (z<0), Bneg",36,0.0,6.3,100,-4.0,4.0); + pPhiZCResidualsHistBneg->SetXTitle("\\phi [rad]"); + pPhiZCResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TH2F* pPtYAResidualsHistBneg = new TH2F("fPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",20,.3,8.,80,-2.0,2.0); + pPtYAResidualsHistBneg->SetXTitle("Pt [rad]"); + pPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtYCResidualsHistBneg = new TH2F("fPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",20,.3,8.,80,-2.0,2.0); + pPtYCResidualsHistBneg->SetXTitle("Pt [rad]"); + pPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtZAResidualsHistBneg = new TH2F("fPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",20,.3,8.,80,-4.0,4.0); + pPtZAResidualsHistBneg->SetXTitle("Pt"); + pPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TH2F* pPtZCResidualsHistBneg = new TH2F("fPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",20,.3,8.,80,-4.0,4.0); + pPtZCResidualsHistBneg->SetXTitle("Pt"); + pPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtYAResidualsHistBneg = new TH2F("fLowPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",100,0.0,1.5,80,-2.0,2.0); + pLowPtYAResidualsHistBneg->SetXTitle("Pt [rad]"); + pLowPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtYCResidualsHistBneg = new TH2F("fLowPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",100,0.0,1.5,80,-2.0,2.0); + pLowPtYCResidualsHistBneg->SetXTitle("Pt [rad]"); + pLowPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtZAResidualsHistBneg = new TH2F("fLowPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",100,0.0,1.5,80,-4.0,4.0); + pLowPtZAResidualsHistBneg->SetXTitle("Pt"); + pLowPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtZCResidualsHistBneg = new TH2F("fLowPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",100,0.0,1.5,80,-4.0,4.0); + pLowPtZCResidualsHistBneg->SetXTitle("Pt"); + pLowPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]"); + TList* listBneg = new TList(); + fList->Add(listBneg); + listBneg->Add(pZYAResidualsHistBneg); + listBneg->Add(pZYCResidualsHistBneg); + listBneg->Add(pLPAResidualsHistBneg); + listBneg->Add(pLPCResidualsHistBneg); + listBneg->Add(pPhiYAResidualsHistBneg); + listBneg->Add(pPhiYCResidualsHistBneg); + listBneg->Add(pPhiZAResidualsHistBneg); + listBneg->Add(pPhiZCResidualsHistBneg); + listBneg->Add(pPtYAResidualsHistBneg); + listBneg->Add(pPtYCResidualsHistBneg); + listBneg->Add(pPtZAResidualsHistBneg); + listBneg->Add(pPtZCResidualsHistBneg); + listBneg->Add(pLowPtYAResidualsHistBneg); + listBneg->Add(pLowPtYCResidualsHistBneg); + listBneg->Add(pLowPtZAResidualsHistBneg); + listBneg->Add(pLowPtZCResidualsHistBneg); + + TH2F* pZYAResidualsHistBnil = new TH2F("fZYAResidualsHistBnil","z-r\\phi residuals side A (z>0), Bnil", 100, -4., 4., 100, -2.0, 2.0 ); + pZYAResidualsHistBnil->SetXTitle("\\deltaz [cm]"); + pZYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pZYCResidualsHistBnil = new TH2F("fZYCResidualsHistBnil","z-r\\phi residuals side C (z<0), Bnil", 100, -4., 4., 100, -2.0, 2.0 ); + pZYCResidualsHistBnil->SetXTitle("\\deltaz [cm]"); + pZYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLPAResidualsHistBnil = new TH2F("fLPAResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bnil", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPAResidualsHistBnil->SetXTitle("\\deltasin(\\phi)"); + pLPAResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pLPCResidualsHistBnil = new TH2F("fLPCResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bnil", 100, -.07, 0.07, 100, -0.07, 0.07 ); + pLPCResidualsHistBnil->SetXTitle("\\deltasin(\\phi)"); + pLPCResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)"); + TH2F* pPhiYAResidualsHistBnil = new TH2F("fPhiYAResidualsHistBnil","\\phi-\\deltar\\phi side A (z>0), Bnil",36,0.0,6.3,100,-2.0,2.0); + pPhiYAResidualsHistBnil->SetXTitle("\\phi [rad]"); + pPhiYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiYCResidualsHistBnil = new TH2F("fPhiYCResidualsHistBnil","\\phi-\\deltar\\phi side C (z<0), Bnil",36,0.0,6.3,100,-2.0,2.0); + pPhiYCResidualsHistBnil->SetXTitle("\\phi [rad]"); + pPhiYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPhiZAResidualsHistBnil = new TH2F("fPhiZAResidualsHistBnil","\\phi-\\deltaz side A (z>0), Bnil",36,0.0,6.3,100,-4.0,4.0); + pPhiZAResidualsHistBnil->SetXTitle("\\phi [rad]"); + pPhiZAResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TH2F* pPhiZCResidualsHistBnil = new TH2F("fPhiZCResidualsHistBnil","\\Pt-\\deltaz side C (z<0), Bnil",36,0.0,6.3,100,-4.0,4.0); + pPhiZCResidualsHistBnil->SetXTitle("\\phi [rad]"); + pPhiZCResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TH2F* pPtYAResidualsHistBnil = new TH2F("fPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",20,.3,8.,80,-2.0,2.0); + pPtYAResidualsHistBnil->SetXTitle("Pt [rad]"); + pPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtYCResidualsHistBnil = new TH2F("fPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",20,.3,8.,80,-2.0,2.0); + pPtYCResidualsHistBnil->SetXTitle("Pt [rad]"); + pPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pPtZAResidualsHistBnil = new TH2F("fPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",20,.3,8.,80,-4.0,4.0); + pPtZAResidualsHistBnil->SetXTitle("Pt"); + pPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TH2F* pPtZCResidualsHistBnil = new TH2F("fPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",20,.3,8.,80,-4.0,4.0); + pPtZCResidualsHistBnil->SetXTitle("Pt"); + pPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtYAResidualsHistBnil = new TH2F("fLowPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",100,0.0,1.5,80,-2.0,2.0); + pLowPtYAResidualsHistBnil->SetXTitle("Pt [rad]"); + pLowPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtYCResidualsHistBnil = new TH2F("fLowPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",100,0.0,1.5,80,-2.0,2.0); + pLowPtYCResidualsHistBnil->SetXTitle("Pt [rad]"); + pLowPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]"); + TH2F* pLowPtZAResidualsHistBnil = new TH2F("fLowPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",100,0.0,1.5,80,-4.0,4.0); + pLowPtZAResidualsHistBnil->SetXTitle("Pt"); + pLowPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TH2F* pLowPtZCResidualsHistBnil = new TH2F("fLowPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",100,0.0,1.5,80,-4.0,4.0); + pLowPtZCResidualsHistBnil->SetXTitle("Pt"); + pLowPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]"); + TList* listBnil = new TList(); + fList->Add(listBnil); + listBnil->Add(pZYAResidualsHistBnil); + listBnil->Add(pZYCResidualsHistBnil); + listBnil->Add(pLPAResidualsHistBnil); + listBnil->Add(pLPCResidualsHistBnil); + listBnil->Add(pPhiYAResidualsHistBnil); + listBnil->Add(pPhiYCResidualsHistBnil); + listBnil->Add(pPhiZAResidualsHistBnil); + listBnil->Add(pPhiZCResidualsHistBnil); + listBnil->Add(pPtYAResidualsHistBnil); + listBnil->Add(pPtYCResidualsHistBnil); + listBnil->Add(pPtZAResidualsHistBnil); + listBnil->Add(pPtZCResidualsHistBnil); + listBnil->Add(pLowPtYAResidualsHistBnil); + listBnil->Add(pLowPtYCResidualsHistBnil); + listBnil->Add(pLowPtZAResidualsHistBnil); + listBnil->Add(pLowPtZCResidualsHistBnil); + TH1F* pNmatchingEff=new TH1F("pNmatchingEff","matching efficiency",50,0.,1.); + fList->Add(pNmatchingEff); + + fAligner = new AliRelAlignerKalman(); + fAligner->SetRejectOutliers(fRejectOutliers); + fAligner->SetOutRejSigma(fOutRejSigma); + fAligner->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median); + fAligner->SetOutRejSigma2Median(fOutRejSigma2Median); + fList->Add(fAligner); + + fDebugTree = new TTree("debugTree","tree with debug info"); + fDebugTree->Branch("aligner","AliRelAlignerKalman",&fAligner); } //________________________________________________________________________ @@ -97,84 +408,259 @@ void AliAnalysisTaskITSTPCalignment::Exec(Option_t *) { // Main loop // Called for each event + + //check for ESD and friends if (!fESD) { - Printf("ERROR: fESD not available"); + AliError("no ESD"); return; } - AliRelAlignerKalman* aligner = fArray->GetAligner(); - - Int_t lastrunnumber = aligner->GetRunNumber(); - UInt_t currentTimeStamp = fESD->GetTimeStamp(); - - //for a new run reset TPC errors - if (lastrunnumber != fESD->GetRunNumber()) + if (!fESDfriend) { - aligner->ResetTPCparamsCovariance(); - fArray->SetCurrentTimeBin(currentTimeStamp); + AliError("no ESD friend"); + return; } - //if time jumps back reset all - if (currentTimeStamp < aligner->GetTimeStamp()) + fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD + + //Update the parmeters + AnalyzeESDevent(); + + // Post output data. + PostData(0, fDebugTree); + PostData(1, fList); + PostData(2, fArrayITSglobal); + PostData(3, fArrayITSsa); +} + +//________________________________________________________________________ +void AliAnalysisTaskITSTPCalignment::AnalyzeESDevent() +{ + //analyze an ESD event with track matching + Int_t ntracks = fESD->GetNumberOfTracks(); + if (ntracks==0) return; + + if (fUseITSoutGlobalTrack) { - aligner->Reset(); - fArray->SetCurrentTimeBin(currentTimeStamp); + AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(fESD); + if (alignerGlobal) + alignerGlobal->AddESDevent(fESD); } - //Update the parmeters - if (fArray->AddCosmicEvent(fESD)) - if (fDoQA) + if (fUseITSoutITSSAtrack) + { + //get the aligner for the event + //do it with matching + TObjArray arrayParamITS(ntracks); + TObjArray arrayParamTPC(ntracks); + + Int_t n = FindMatchingTracks(arrayParamITS, arrayParamTPC, fESD); + AliInfo(Form("n matched: %i\n",n)); + + TH1F* nMatchingEff=dynamic_cast(fList->At(3)); + Float_t ratio = (float)n/(float)ntracks; + nMatchingEff->Fill(ratio); + + AliRelAlignerKalman* alignerSA = NULL; + if (n>0) alignerSA = fArrayITSsa->GetAligner(fESD); + + for (Int_t i=0;iFindCosmicTrackletNumbersInEvent( - trackTArrITS, trackTArrTPC, fESD )) + AliExternalTrackParam* paramsITS=dynamic_cast(arrayParamITS[i]); + AliExternalTrackParam* paramsTPC=dynamic_cast(arrayParamTPC[i]); + + //QA + if (fDoQA) DoQA(paramsITS,paramsTPC); + + //debug + if (fAligner->AddTrackParams(paramsITS, paramsTPC)) { - AliESDtrack* ptrack; - const AliExternalTrackParam* pconstparams1; - const AliExternalTrackParam* pconstparams2; - AliExternalTrackParam params1; - AliExternalTrackParam params2; - - //////////////////////////////// - for (Int_t i=0;iGetTrack(trackTArrITS[i]); - pconstparams1 = ptrack->GetOuterParam(); - if (!pconstparams1) continue; - params1 = *pconstparams1; //make copy to be safe - - //TPC track - ptrack = fESD->GetTrack(trackTArrTPC[i]); - pconstparams2 = ptrack->GetInnerParam(); - if (!pconstparams2) continue; - params2 = *pconstparams2; //make copy - params2.Rotate(params1.GetAlpha()); - params2.PropagateTo( params1.GetX(), aligner->GetMagField() ); - - Float_t resy = params2.GetY() - params1.GetY(); - Float_t resz = params2.GetZ() - params1.GetZ(); - Float_t ressnp = params2.GetSnp() - params1.GetSnp(); - Float_t restgl = params2.GetTgl() - params1.GetTgl(); - fYZResidualsHist->Fill(resy,resz); - fPLResidualsHist->Fill(ressnp,restgl); - } - }//if DoQA - - }//if AddEvent + fAligner->SetRunNumber(fESD->GetRunNumber()); + fAligner->SetMagField(fESD->GetMagneticField()); + fAligner->SetTimeStamp(fESD->GetTimeStamp()); + } - // Post output data. - PostData(0, fArray); - PostData(1, fListOfHistos); + //alignment check + if (alignerSA) alignerSA->AddTrackParams(paramsITS, paramsTPC); + } + arrayParamITS.Delete(); + arrayParamTPC.Delete(); + } } //________________________________________________________________________ void AliAnalysisTaskITSTPCalignment::Terminate(Option_t *) { // Called once at the end of the query - fArray->Merge(new TList()); //final cleanup +} + +//________________________________________________________________________ +Int_t AliAnalysisTaskITSTPCalignment::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD) +{ + //find matching tracks and return tobjarrays with the params + //fUniqueID of param object contains tracknumber in event. + + Int_t ntracks = pESD->GetNumberOfTracks(); + Double_t magfield = pESD->GetMagneticField(); + + Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found + for (Int_t i=0;iGetTrack(i); + if (!track1) continue; + + if (track1->GetNcls(0) < fMinPointsVol1) continue; + + if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) && + (!track1->IsOn(AliESDtrack::kTPCin)) )) continue; + + const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack(); + if (!constfriendtrack1) {AliInfo(Form("no friendTrack1\n"));continue;} + AliESDfriendTrack friendtrack1(*constfriendtrack1); + + if (!friendtrack1.GetITSOut()) continue; + AliExternalTrackParam params1(*(friendtrack1.GetITSOut())); + + Double_t bestd = 1000.; //best distance + Bool_t newi = kTRUE; //whether we start with a new i + for (Int_t j=0; j0 && matchedArr[j]!=i) continue; //already matched, everything tried + //get track2 and friend + AliESDtrack* track2 = pESD->GetTrack(j); + if (!track2) continue; + if (track1==track2) continue; + if (!(track2->IsOn(AliESDtrack::kTPCin))) continue; //must be TPC + if (!(track2->IsOn(AliESDtrack::kITSout))) continue; //must have ITS + + //if (track2->GetNcls(0) != track1->GetNcls(0)) continue; + //if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue; + if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC + if (track2->GetTgl() > 1.) continue; //acceptance + //cut crossing tracks + if (!track2->GetInnerParam()) continue; + if (!track2->GetOuterParam()) continue; + if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue; + if (track2->GetInnerParam()->GetX()>90) continue; + if (TMath::Abs(track2->GetInnerParam()->GetZ())<5.) continue; //too close to membrane? + + AliExternalTrackParam params2(*(track2->GetInnerParam())); + + //bring to same reference plane + if (!params2.Rotate(params1.GetAlpha())) continue; + if (!params2.PropagateTo(params1.GetX(), magfield)) continue; + + //pt cut, only for data with magfield + if ((params2.Pt()1.)) continue; + //else + //if (fMC) + // { + // AliStack* stack = fMC->Stack(); + // Int_t label = track2->GetLabel(); + // if (label<0) continue; + // TParticle* particle = stack->Particle(label); + // if (particle->Pt() 2.0) continue; + if (dz > 10.0) continue; + if (dphi > 0.1 ) continue; + if (dlam > 0.1 ) continue; + + //best match only + Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam); + if ( d >= bestd) continue; + bestd = d; + matchedArr[j]=i; //j-th track matches i-th (ITS) track + if (newi) iMatched++; newi=kFALSE; //increment at most once per i + params1.SetUniqueID(i); //store tracknummer + params2.SetUniqueID(j); + if (arrITS[iMatched] && arrTPC[iMatched]) + { + *(arrITS[iMatched]) = params1; + *(arrTPC[iMatched]) = params2; + } + else + { + arrITS[iMatched] = new AliExternalTrackParam(params1); + arrTPC[iMatched] = new AliExternalTrackParam(params2); + }//else + }//for j + }//for i + return iMatched+1; +} + +//________________________________________________________________________ +void AliAnalysisTaskITSTPCalignment::DoQA(AliExternalTrackParam* paramsITS, + AliExternalTrackParam* paramsTPC) +{ + //fill qa histograms in a given list, per field direction (5,-5,0) + Float_t resy = paramsTPC->GetY() - paramsITS->GetY(); + Float_t resz = paramsTPC->GetZ() - paramsITS->GetZ(); + Float_t ressnp = paramsTPC->GetSnp() - paramsITS->GetSnp(); + Float_t restgl = paramsTPC->GetTgl() - paramsITS->GetTgl(); + + TList* pList=NULL; + Double_t field = fESD->GetMagneticField(); + if (field >= 1.) pList = dynamic_cast(fList->At(0)); + if (field <= -1.) pList = dynamic_cast(fList->At(1)); + if (field < 1. && field > -1.) pList = dynamic_cast(fList->At(2)); + if (!pList) return; + + TH2F* pZYAResidualsHist = dynamic_cast(pList->At(0)); + TH2F* pZYCResidualsHist = dynamic_cast(pList->At(1)); + TH2F* pLPAResidualsHist = dynamic_cast(pList->At(2)); + TH2F* pLPCResidualsHist = dynamic_cast(pList->At(3)); + TH2F* pPhiYAResidualsHist = dynamic_cast(pList->At(4)); + TH2F* pPhiYCResidualsHist = dynamic_cast(pList->At(5)); + TH2F* pPhiZAResidualsHist = dynamic_cast(pList->At(6)); + TH2F* pPhiZCResidualsHist = dynamic_cast(pList->At(7)); + TH2F* pPtYAResidualsHist = dynamic_cast(pList->At(8)); + TH2F* pPtYCResidualsHist = dynamic_cast(pList->At(9)); + TH2F* pPtZAResidualsHist = dynamic_cast(pList->At(10)); + TH2F* pPtZCResidualsHist = dynamic_cast(pList->At(11)); + TH2F* pLowPtYAResidualsHist = dynamic_cast(pList->At(12)); + TH2F* pLowPtYCResidualsHist = dynamic_cast(pList->At(13)); + TH2F* pLowPtZAResidualsHist = dynamic_cast(pList->At(14)); + TH2F* pLowPtZCResidualsHist = dynamic_cast(pList->At(15)); + + if (paramsITS->GetZ() > 0.0) + { + pZYAResidualsHist->Fill(resz,resy); + pLPAResidualsHist->Fill(restgl,ressnp); + pPhiYAResidualsHist->Fill(paramsTPC->Phi(),resy); + pPhiZAResidualsHist->Fill(paramsTPC->Phi(),resz); + pPtYAResidualsHist->Fill(paramsTPC->Pt(),resy); + pPtZAResidualsHist->Fill(paramsTPC->Pt(),resz); + pLowPtYAResidualsHist->Fill(paramsTPC->Pt(),resy); + pLowPtZAResidualsHist->Fill(paramsTPC->Pt(),resz); + } + else + { + pZYCResidualsHist->Fill(resz,resy); + pLPCResidualsHist->Fill(restgl,ressnp); + pPhiYCResidualsHist->Fill(paramsTPC->Phi(),resy); + pPhiZCResidualsHist->Fill(paramsTPC->Phi(),resz); + pPtYCResidualsHist->Fill(paramsTPC->Pt(),resy); + pPtZCResidualsHist->Fill(paramsTPC->Pt(),resz); + pLowPtYCResidualsHist->Fill(paramsTPC->Pt(),resy); + pLowPtZCResidualsHist->Fill(paramsTPC->Pt(),resz); + } } diff --git a/PWG1/AliAnalysisTaskITSTPCalignment.h b/PWG1/AliAnalysisTaskITSTPCalignment.h index 8a63eccbe8e..006908e8848 100644 --- a/PWG1/AliAnalysisTaskITSTPCalignment.h +++ b/PWG1/AliAnalysisTaskITSTPCalignment.h @@ -8,9 +8,13 @@ // Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch /////////////////////////////////////////////////////////////////////////// -#include +#include +class TList; class TTree; +class AliExternalTrackParam; class AliESDEvent; +class AliESDfriend; +class AliMCEvent; class AliRelAlignerKalman; class AliRelAlignerKalmanArray; class TH2F; @@ -22,30 +26,62 @@ public: AliAnalysisTaskITSTPCalignment(const char *name); virtual ~AliAnalysisTaskITSTPCalignment() {} - void SetSaveInterval( const UInt_t si ) {fSaveInterval = si;} - void SetTimeMatchingTolerance( const UInt_t tol ) {fTimeMatchingTolerance = tol; } - void SetDoQA( const Bool_t qa ) {fDoQA=qa;} + void SetupAlignerArray( Int_t t0, Int_t tend, Int_t slotwidth ) + { fT0=t0; fTend=tend; fSlotWidth=slotwidth; } + void SetFillDebugTree(Bool_t m=kTRUE) {fFillDebugTree=m;} + void SetDoQA(Bool_t d=kTRUE) {fDoQA=d;} + void SetMinPt(Double_t m) {fMinPt=m;} + void SetMinNclsITS(Int_t m) {fMinPointsVol1=m;} + void SetMinNclsTPC(Int_t m) {fMinPointsVol2=m;} + void DoQA(AliExternalTrackParam* paramsITS,AliExternalTrackParam* paramsTPC); + void SetRejectOutliers(Bool_t set=kTRUE){fRejectOutliers=set;} + void SetRejectOutliersSigma2Median(Bool_t set=kTRUE){fRejectOutliersSigma2Median=set;} + void SetOutRejSigma(Double_t d){fOutRejSigma=d;} + void SetOutRejSigma2Median(Double_t d){fOutRejSigma2Median=d;} + void SetOutRejSigmaOnMerge(Double_t d){fOutRejSigmaOnMerge=d;} + void SetUseITSoutGlobalTrack(Bool_t b){fUseITSoutGlobalTrack=b;} + void SetUseITSoutITSSAtrack(Bool_t b){fUseITSoutITSSAtrack=b;} + + Int_t FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD); + void AnalyzeESDevent(); virtual void ConnectInputData(Option_t *); virtual void CreateOutputObjects(); virtual void Exec(Option_t *option); virtual void Terminate(Option_t *); + virtual Bool_t Notify(); private: - AliESDEvent* fESD; //ESD object - AliRelAlignerKalmanArray* fArray; //array of aligners - TH2F* fYZResidualsHist; //2D histogram with the yz residuals - TH2F* fPLResidualsHist; //2D histogram with the phi lambda residuals - TList* fListOfHistos; //list with QA histograms - UInt_t fSaveInterval; //save interveal - UInt_t fTimeMatchingTolerance; //time matching tolerance - - Bool_t fDoQA; //whether to fill QA histograms + AliESDEvent* fESD; //! ESD object + AliESDfriend* fESDfriend; //! ESD friend + AliMCEvent* fMC; //! mc event + AliRelAlignerKalmanArray* fArrayITSglobal; //! array of aligners with ITS global + AliRelAlignerKalmanArray* fArrayITSsa; //! array of aligners ITS standalone + TTree* fDebugTree; //! tree + AliRelAlignerKalman* fAligner; //! aligner + TList* fList; //! list with QA histograms + + Bool_t fFillDebugTree; // do we write the debug tree? + Bool_t fDoQA; // do QA? + Int_t fT0; // t0 + Int_t fTend; // tend + Int_t fSlotWidth; // slotwidth + Double_t fMinPt; // min pt for tracks + Int_t fMinPointsVol1; // min clusters its + Int_t fMinPointsVol2; // min clusters tpc + Bool_t fRejectOutliers; // reject outliers? + Double_t fOutRejSigma; // how many outliers + Bool_t fRejectOutliersSigma2Median; // + Double_t fOutRejSigma2Median; // + Double_t fOutRejSigmaOnMerge; // + Bool_t fUseITSoutGlobalTrack; // + Bool_t fUseITSoutITSSAtrack; // AliAnalysisTaskITSTPCalignment(const AliAnalysisTaskITSTPCalignment&); // not implemented AliAnalysisTaskITSTPCalignment& operator=(const AliAnalysisTaskITSTPCalignment&); // not implemented - ClassDef(AliAnalysisTaskITSTPCalignment, 1); + ClassDef(AliAnalysisTaskITSTPCalignment, 2); }; #endif + diff --git a/PWG1/AliRelAlignerKalmanArray.cxx b/PWG1/AliRelAlignerKalmanArray.cxx index 1f9b88a8621..85918d2a540 100644 --- a/PWG1/AliRelAlignerKalmanArray.cxx +++ b/PWG1/AliRelAlignerKalmanArray.cxx @@ -24,10 +24,13 @@ ////////////////////////////////////////////////////////////////////////////// #include +#include +#include +#include #include #include #include -#include "AliESDEvent.h" +#include #include "AliRelAlignerKalman.h" #include "AliRelAlignerKalmanArray.h" @@ -35,240 +38,391 @@ ClassImp(AliRelAlignerKalmanArray) //______________________________________________________________________________ AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(): - TNamed(), - fArray(new TObjArray()), - fSaveInterval(600), //default every 10minutes - fTimeMatchingTolerance(20), - fCurrentTimeBin(0), - fAligner(new AliRelAlignerKalman()) + TNamed("alignerArray","array of aligners"), + fT0(0), + fTimebinWidth(0), + fSize(0), + fOutRejSigmaOnMerge(10.), + fOutRejSigmaOnSmooth(1.), + fPArray(NULL) { //ctor } //______________________________________________________________________________ -AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(const char* name): - TNamed(name, name), - fArray(new TObjArray()), - fSaveInterval(600), //default every 10 minutes - fTimeMatchingTolerance(60), - fCurrentTimeBin(0), - fAligner(new AliRelAlignerKalman()) +AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(Int_t t0, Int_t tend, Int_t slotwidth): + TNamed("alignerArray","array of aligners"), + fT0(t0), + fTimebinWidth(slotwidth), + fSize(0), + fOutRejSigmaOnMerge(10.), + fOutRejSigmaOnSmooth(1.), + fPArray(NULL) { //ctor + if (slotwidth==0) fSize = 1; + else fSize = (tend-t0)/slotwidth; + fPArray = new AliRelAlignerKalman* [fSize]; + if (fPArray) for (Int_t i=0;iSetOwner(); - delete fArray; - delete fAligner; + ClearContents(); + delete [] fPArray; +} + +//______________________________________________________________________________ +void AliRelAlignerKalmanArray::SetupArray(Int_t t0, Int_t tend, Int_t slotwidth) +{ + //setup array - clears old content + ClearContents(); + fT0=t0; + fTimebinWidth=slotwidth; + Int_t tmpsize; + if (slotwidth==0) tmpsize = 1; + else tmpsize = (tend-t0)/slotwidth; + if (tmpsize!=fSize) + { + delete [] fPArray; + fPArray=new AliRelAlignerKalman* [tmpsize]; + if (fPArray) fSize=tmpsize; + else fSize=0; + } + for (Int_t i=0;i(next())) ) { - if (arrayFromList==this) continue; - - fArray->AddAll(arrayFromList->fArray); //put all objects in one array - - //do the merge - fArray->Sort(); - TObjArray* tmpArray = SortedMerge(fArray); - tmpArray->SetOwner(kTRUE); - - TObjArray* newArray = dynamic_cast(tmpArray->Clone()); - delete fArray; //takes care of all loaded objects - fArray = newArray; + if (fT0 != arrayFromList->fT0) continue; + if (fTimebinWidth != arrayFromList->fTimebinWidth) continue; + if (fSize != arrayFromList->fSize) continue; - fArray->AddLast(arrayFromList->fAligner); //add the endofrun aligner + for (Int_t i=0; ifPArray[i]; + if (a1 && a2) + { + a1->SetRejectOutliers(kFALSE); + a1->SetOutRejSigma(fOutRejSigmaOnMerge); a1->Merge(a2); + } + else + if (!a1 && a2) fPArray[i] = new AliRelAlignerKalman(*a2); + } } + return 0; +} - //TODO: this can be done better! - //Add own endofrun aligner and clean up - fArray->AddLast(fAligner); - fArray->Sort(); - //TObjArray* tmpArray = SortedMerge(fArray); - //tmpArray->SetOwner(kTRUE); - //TObjArray* newArray = dynamic_cast(tmpArray->Clone()); - //delete fArray; //takes care of all loaded objects - //fArray = newArray; - - return fArray->GetEntriesFast(); +//______________________________________________________________________________ +Int_t AliRelAlignerKalmanArray::Timebin( UInt_t timestamp ) const +{ + //calculate binnumber for given timestamp + if (fTimebinWidth==0) return 0; + Int_t slot = (timestamp-fT0)/fTimebinWidth; //it's all integers! + return slot; } //______________________________________________________________________________ -TObjArray* AliRelAlignerKalmanArray::SortedMerge( TObjArray* input ) +AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts) { - //Merges the adjacent aligners if close enough - //input needs to be already sorted - - UInt_t timeStampIn; - AliRelAlignerKalman* alignerIn; - AliRelAlignerKalman* alignerOut = dynamic_cast(input->At(0)); - TObjArray* output = new TObjArray(); //empty array - output->AddLast(alignerOut); //first object in: copy of first input element - - timeStampIn = alignerOut->GetTimeStamp(); - SetCurrentTimeBin( timeStampIn ); - TIter next(input); - while ( (alignerIn = dynamic_cast(next())) ) + //get the aligner for specified timestamp + Int_t tb = Timebin(ts); + if (tb<0) return NULL; + if (tb>=fSize) return NULL; + if (!fPArray) return NULL; + if (!fPArray[tb]) { - timeStampIn = alignerIn->GetTimeStamp(); - if ( IsInCurrentTimeBin(timeStampIn) ) - { - alignerOut->Merge(alignerIn); - } - else - { - alignerOut = alignerIn; - output->AddLast(alignerOut); - }//if - SetCurrentTimeBin( timeStampIn ); + fPArray[tb] = new AliRelAlignerKalman( *fAlignerTemplate ); + fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth); } - return output; + return fPArray[tb]; } //______________________________________________________________________________ -void AliRelAlignerKalmanArray::SetCurrentTimeBin( UInt_t timestamp ) +AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event) { - //set the current timebin - fCurrentTimeBin = TimeBin(timestamp); + //get the aligner for this event and set relevant info + AliRelAlignerKalman* a = GetAligner(event->GetTimeStamp()); + if (a) a->SetRunNumber(event->GetRunNumber()); + if (a) a->SetMagField(event->GetMagneticField()); + return a; } //______________________________________________________________________________ -UInt_t AliRelAlignerKalmanArray::TimeBin( UInt_t timestamp ) const +AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in) { - return (timestamp+(fSaveInterval/2))/fSaveInterval*fSaveInterval; //it's all integers! + //assignment operator + if (fSize!=in.fSize) + { + //if sizes different, delete curent and make a new one with new size + ClearContents(); + delete [] fPArray; + fPArray = new AliRelAlignerKalman* [in.fSize]; + } + + fOutRejSigmaOnMerge=in.fOutRejSigmaOnMerge; + fOutRejSigmaOnSmooth=in.fOutRejSigmaOnSmooth; + + if (!fPArray) fSize=0; + else fSize = in.fSize; + fTimebinWidth = in.fTimebinWidth; + fT0 = in.fT0; + for (Int_t i=0;i=fCurrentTimeBin)?timestamp-fCurrentTimeBin: - fCurrentTimeBin-timestamp; - return (timeDiff < fTimeMatchingTolerance); + //clear contents of array + for (Int_t i=0;iGetRunNumber() != fAligner->GetRunNumber()) -// { -// //what to do when a new run starts -// } -// -//} -// -////______________________________________________________________________________ -Bool_t AliRelAlignerKalmanArray::AddCosmicEvent( AliESDEvent* event ) +//______________________________________________________________________________ +AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const { - if (!fAligner->AddCosmicEvent(event)) return kFALSE; - - UInt_t currentTimeStamp = event->GetTimeStamp(); - UInt_t timeFromLastBinCentre = currentTimeStamp - fCurrentTimeBin; - UInt_t binCentre = TimeBin(currentTimeStamp); - UInt_t timeFromBinCentre = currentTimeStamp-binCentre; - UInt_t nIntervals = timeFromLastBinCentre/fSaveInterval; - - ////////////////////////////////////////////////////////////////////////////// - //only ONE SAVE PER TIMEBIN!!!!, as close as possible to the bin center - ////////////////////////////////////////////////////////////////////////////// - if ( (nIntervals == 1) && //is in next time bin passed centre - (timeFromBinCentre < fTimeMatchingTolerance) ) //and close to it - { - AddLast(new AliRelAlignerKalman(*fAligner)); - } - else if ( (nIntervals > 2) ) //TODO: don't hardwire stuff! - { - //if missed a few windows save anyway at current bin centre - fAligner->SetTimeStamp(binCentre); - AddLast(new AliRelAlignerKalman(*fAligner)); - } - ////////////////////////////////////////////////////////////////////////////// + //mimic TObjArray::At( Int_t i ) + if (i>=fSize) {printf("AliRelAlignerKalmanArray::At index %i out of bounds, size=%i\n",i,fSize); return NULL;} + if (i<0) return NULL; + return fPArray[i]; +} - return kTRUE; +//______________________________________________________________________________ +AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const +{ + //mimic TObjArray::operator[](Int_t) + if (i>=fSize) {printf("AliRelAlignerKalmanArray::operator[] index %i out of bounds, size=%i\n",i,fSize); return NULL;} + if (i<0) return NULL; + return fPArray[i]; } //______________________________________________________________________________ -AliRelAlignerKalmanArray::AliRelAlignerKalmanArray( const AliRelAlignerKalmanArray& in): - TNamed(in.GetName(), in.GetTitle()), - fArray(NULL), - fSaveInterval(in.fSaveInterval), - fTimeMatchingTolerance(in.fTimeMatchingTolerance), - fCurrentTimeBin(in.fCurrentTimeBin), - fAligner(new AliRelAlignerKalman(*in.fAligner)) +AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i) { - //copy ctor - fArray = static_cast(in.Clone()); + //mimic TObjArray::operator[](Int_t) can be used as lvalue + if (i>=fSize||i<0) {Error("operator[]", "index %i out of bounds, size=%i\n",i,fSize); return fPArray[0];} + return fPArray[i]; } //______________________________________________________________________________ -AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in) +AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const { - //assignment operator - fArray = static_cast(in.Clone()); - fSaveInterval = in.fSaveInterval; - fTimeMatchingTolerance = in.fTimeMatchingTolerance; - return *this; + //return aligner in last filled slot + if (fSize==0) return NULL; + return fPArray[fSize-1]; } //______________________________________________________________________________ -Bool_t AliRelAlignerKalmanArray::SetSaveInterval( const UInt_t s ) +void AliRelAlignerKalmanArray::Print(Option_t* option) const { - //only set if array empty - if (fArray->GetEntriesFast()) return kFALSE; - fSaveInterval = s; - return kTRUE; + // print some info + TString optionStr(option); + printf( "%s: t0: %i, tend: %i, width: %i, size: %i, entries: %i\n", + GetName(), + fT0, (fT0+(fSize)*fTimebinWidth), fTimebinWidth, + fSize, GetEntries() ); + if (optionStr.Contains("a")) + for (Int_t i=0; iPrint(); + } } //______________________________________________________________________________ -Bool_t AliRelAlignerKalmanArray::SetTimeMatchingTolerance( const UInt_t m ) +Int_t AliRelAlignerKalmanArray::GetEntries() const { - //only set if array empty - if (fArray->GetEntriesFast()) return kFALSE; - fTimeMatchingTolerance = m; - return kTRUE; + //get number of filled slots + if (!fPArray) return 0; + Int_t entries=0; + for (Int_t i=0; i(fArray->At(i)); + AliRelAlignerKalman* al = NULL; + tree->Branch("aligner","AliRelAlignerKalman",&al); + //fill a tree with filled slots + for (Int_t i=0; iFill(); + } } //______________________________________________________________________________ -void AliRelAlignerKalmanArray::AddLast( AliRelAlignerKalman* al ) +TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const { - //mimic TObjArray::AddLast( TObject* obj ) - fArray->AddLast( al ); - SetCurrentTimeBin(al->GetTimeStamp()); + //return a graph for the iparam-th parameter + if (iparam>8 || iparam<0) + { + printf("wrong parameter number. must be from 0-8"); + return NULL; + } + + TVectorD vx(GetEntries()); + TVectorD vy(GetEntries()); + TVectorD vex(GetEntries()); + TVectorD vey(GetEntries()); + Int_t entry=0; + for (Int_t i=0; iGetTimeStamp(); + vy(entry) = al->GetStateArr()[iparam]; + TMatrixDSym* cm = al->GetStateCov(); + vey(entry) = TMath::Sqrt((*cm)(iparam,iparam)); + entry++; + } + + TGraphErrors* graph = new TGraphErrors(vx,vy,vex,vey); + TString graphtitle; + TString graphtitley; + switch (iparam) + { + case 0: + graphtitle="rotation \\psi"; + graphtitley="\\psi [rad]"; + break; + case 1: + graphtitle="rotation \\theta"; + graphtitley="\\theta [rad]"; + break; + case 2: + graphtitle="rotation \\phi"; + graphtitley="\\phi [rad]"; + break; + case 3: + graphtitle="shift x"; + graphtitley="x [cm]"; + break; + case 4: + graphtitle="shift y"; + graphtitley="y [cm]"; + break; + case 5: + graphtitle="shift z"; + graphtitley="z [cm]"; + break; + case 6: + graphtitle="TPC vd correction"; + graphtitley="correction factor []"; + break; + case 7: + graphtitle="TPC t0 correction"; + graphtitley="t0 correction [\\micros]"; + break; + case 8: + graphtitle="TPC dv/dy"; + graphtitley="dv/dy [cm/\\micros/m]"; + break; + } + graph->SetTitle(graphtitle); + TAxis* xas = graph->GetXaxis(); + TAxis* yas = graph->GetYaxis(); + xas->SetTitle("time"); + xas->SetTimeDisplay(1); + yas->SetTitle(graphtitley); + return graph; } //______________________________________________________________________________ -AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const +AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const { - //mimic TObjArray::Last() - return dynamic_cast(fArray->Last()); + //return a smoothed version of the data + AliRelAlignerKalmanArray* outputarr = new AliRelAlignerKalmanArray(fT0, + fT0+fSize*fTimebinWidth,fTimebinWidth); + + AliRelAlignerKalman* tmpaligner = outputarr->GetAlignerTemplate(); + tmpaligner->SetOutRejSigma(fOutRejSigmaOnSmooth); + //copy the first filled slot + Int_t n=0; + while (!fPArray[n]) {n++;} + if (n==fSize) return NULL; + *tmpaligner = *fPArray[n]; + if (fPArray[n]->GetNUpdates()>10) + (*outputarr)[n] = new AliRelAlignerKalman(*(fPArray[n])); + //for the rest of slots use merge + for (Int_t i=n+1;iGetTimeStamp()); + if (tmpaligner->Merge(fPArray[i])) + (*outputarr)[i] = new AliRelAlignerKalman(*tmpaligner); + else + (*outputarr)[i] = new AliRelAlignerKalman(*(fPArray[i])); + } + + tmpaligner->Reset(); //clean up + + return outputarr; } //______________________________________________________________________________ -AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const +void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const { - //mimic TObjArray::operator[](Int_t) - return dynamic_cast(fArray->At(i)); + //propagate an aligner in time + TMatrixDSym* cov = al->GetStateCov(); + TMatrixDSym corr(TMatrixDSym::kZero,*cov); + UInt_t dt = (timestamp>al->GetTimeStamp())? + timestamp-al->GetTimeStamp():al->GetTimeStamp()-timestamp; + //the propagation matrix + //to be added to the covariance matrix of kalman filter + corr(6,6) = dt*1e-8; //vdrift + (*cov) += corr; } diff --git a/PWG1/AliRelAlignerKalmanArray.h b/PWG1/AliRelAlignerKalmanArray.h index e3e19abf4c1..a0f8951c07c 100644 --- a/PWG1/AliRelAlignerKalmanArray.h +++ b/PWG1/AliRelAlignerKalmanArray.h @@ -10,58 +10,57 @@ // ////////////////////////////////////////////////////////////////////////////// -#ifndef AliRelAlignerKalmanArray_h -#define AliRelAlignerKalmanArray_h - - -class TString; -class TCollection; -class AliESDEvent; class TObjArray; +class TGraphErrors; class AliRelAlignerKalman; class TNamed; +class TTree; +class TCollection; +class AliESDEvent; class AliRelAlignerKalmanArray:public TNamed { public: AliRelAlignerKalmanArray(); - AliRelAlignerKalmanArray(const char* name); + AliRelAlignerKalmanArray(Int_t t0, Int_t tend, Int_t slotwidth); virtual ~AliRelAlignerKalmanArray(); AliRelAlignerKalmanArray& operator=(const AliRelAlignerKalmanArray& a ); AliRelAlignerKalmanArray(const AliRelAlignerKalmanArray& a); + void SetupArray(Int_t t0, Int_t tend, Int_t slotwidth); + AliRelAlignerKalman* GetAligner(UInt_t timestamp); + AliRelAlignerKalman* GetAligner(AliESDEvent* event); + AliRelAlignerKalman* GetAlignerTemplate(); Long64_t Merge( TCollection* list ); - //Bool_t AddESDEvent( AliESDEvent* event ); - Bool_t AddCosmicEvent( AliESDEvent* event ); - void AddLast( AliRelAlignerKalman* al ); AliRelAlignerKalman* At( Int_t i ) const; - AliRelAlignerKalman* Last() const; - Int_t GetEntries() const {return fArray->GetEntriesFast();} AliRelAlignerKalman* operator[](Int_t i) const; - Bool_t SetTimeMatchingTolerance( const UInt_t m ); - Bool_t SetSaveInterval( const UInt_t s ); - UInt_t GetTimeMatchingTolerance() const {return fTimeMatchingTolerance;} - UInt_t GetSaveInterval() const {return fSaveInterval;} - UInt_t TimeBin( UInt_t timebin ) const; - void SetCurrentTimeBin( UInt_t timestamp ); - UInt_t GetCurrentTimeBin() const {return fCurrentTimeBin;} - Bool_t IsInCurrentTimeBin( UInt_t timestamp ) const; - AliRelAlignerKalman* GetAligner() const {return fAligner;} - TObjArray* SortedMerge ( TObjArray* input ); - //void SetResetAllAtNewRun( Bool_t s ) {fResetAllAtNewRun = s;} - //void SetResetTPCAtNewRun( Bool_t s ) {fResetTPCAtNewRun = s;} + AliRelAlignerKalman*& operator[](Int_t i); + Int_t GetEntries() const; + Int_t GetSize() const {return fSize;} + AliRelAlignerKalman* Last() const; + UInt_t GetT0() const {return fT0;} + UInt_t GetTimebinWidth() const {return fTimebinWidth;} + Int_t Timebin( UInt_t timestamp ) const; + virtual void Print(Option_t* option="") const; + void FillTree( TTree* tree )const ; + TGraphErrors* MakeGraph(Int_t iparam) const; + AliRelAlignerKalmanArray* MakeSmoothArray() const; + void SetOutRejSigmaOnMerge(Double_t s) {fOutRejSigmaOnMerge=s;} + void SetOutRejSigmaOnSmooth(Double_t s) {fOutRejSigmaOnSmooth=s;} private: - TObjArray* fArray; //an array of aligners - UInt_t fSaveInterval; //how often to save (in seconds) - UInt_t fTimeMatchingTolerance; //tolerance for matching timestamps - UInt_t fCurrentTimeBin; //current timebin - AliRelAlignerKalman* fAligner; //aligner object - //Bool_t fResetAllAtNewRun; - //Bool_t fResetTPCAtNewRun; + void ClearContents(); + void PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const; + + UInt_t fT0; //time of first time slot + Int_t fTimebinWidth; //width of the time bin in seconds + Int_t fSize; //size + Double_t fOutRejSigmaOnMerge; //how much outlier rejection on merge + Double_t fOutRejSigmaOnSmooth; //how much outlier rejection on Smooth + AliRelAlignerKalman* fAlignerTemplate; //template + AliRelAlignerKalman** fPArray; //[fSize] an array of aligners - ClassDef(AliRelAlignerKalmanArray,1) //AliRelAlignerKalman class + ClassDef(AliRelAlignerKalmanArray,4); //AliRelAlignerKalman class }; -#endif diff --git a/PWG1/macros/AddTaskITSTPCalignment.C b/PWG1/macros/AddTaskITSTPCalignment.C new file mode 100644 index 00000000000..9b104a2507b --- /dev/null +++ b/PWG1/macros/AddTaskITSTPCalignment.C @@ -0,0 +1,80 @@ +AliAnalysisTaskITSTPCalignment *AddTaskITSTPCalignment() +{ + //add the ITS TPC alignemtn task to the manager + //Mikolaj Krzewicki, mikolaj@nikhef.nl + + //______________________________________________________________________________ + // Get the pointer to the existing analysis manager via the static access method. + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) + { + ::Error("AddTaskITSTPCalignment", "No analysis manager to connect to."); + return NULL; + } + + //______________________________________________________________________________ + // Check the analysis type using the event handlers connected to the analysis manager. + if (!mgr->GetInputEventHandler()) + { + ::Error("AddTaskITSTPCalignment", "This task requires an input event handler"); + return NULL; + } + TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + if (type!="ESD") + { + ::Error("AddTaskITSTPCalignment", "This task only works on ESDs"); + return NULL; + } + + //___________________________________________________________________________ + // Create the task, add it to manager and configure it. + AliAnalysisTaskITSTPCalignment *task = new AliAnalysisTaskITSTPCalignment("taskITSTPCalignment"); + TTimeStamp t0(2009,11,1,0,0,0); + TTimeStamp tend(2010,12,31,0,0,0); + Int_t slotwidth = 3600; + task->SetupAlignerArray(t0.GetSec(),tend.GetSec(),slotwidth); + task->SetFillDebugTree(kFALSE); + task->SetDoQA(kTRUE); + task->SetMinPt(0.4); + task->SetMinNclsITS(4); + task->SetMinNclsTPC(80); + task->SetRejectOutliers(kTRUE); //internal KF outlier rejection (kalman update-based) + task->SetRejectOutliersSigma2Median(kTRUE); //input data outlier removal + task->SetOutRejSigma(1.); //max distance the kf state is allowed to jump + task->SetOutRejSigmaOnMerge(50.); //outlier rejection when merging vertically + task->SetOutRejSigma2Median(5.); //max distance from median for input data + task->SetUseITSoutGlobalTrack(kFALSE); + task->SetUseITSoutITSSAtrack(kTRUE); + + mgr->AddTask(task); + + //______________________________________________________________________________ + //connect output + TString outputFilename = "outputITSTPCalignment.root"; + + AliAnalysisDataContainer* coutput0 = mgr->CreateContainer("outputTree", + TTree::Class(), + AliAnalysisManager::kOutputContainer, + outputFilename.Data()); + AliAnalysisDataContainer* coutput1 = mgr->CreateContainer("outputList", + TList::Class(), + AliAnalysisManager::kOutputContainer, + outputFilename.Data()); + AliAnalysisDataContainer* coutput2 = mgr->CreateContainer("outputArrayITSglobal", + AliRelAlignerKalmanArray::Class(), + AliAnalysisManager::kOutputContainer, + outputFilename.Data()); + AliAnalysisDataContainer* coutput3 = mgr->CreateContainer("outputArrayITSSA", + AliRelAlignerKalmanArray::Class(), + AliAnalysisManager::kOutputContainer, + outputFilename.Data()); + + mgr->ConnectInput(task,0,mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task,0,coutput0); + mgr->ConnectOutput(task,1,coutput1); + mgr->ConnectOutput(task,2,coutput2); + mgr->ConnectOutput(task,3,coutput3); + + return task; +} + diff --git a/PWG1/macros/plotAnalysisTaskITSTPCalignment.C b/PWG1/macros/plotAnalysisTaskITSTPCalignment.C new file mode 100644 index 00000000000..9f49deaf640 --- /dev/null +++ b/PWG1/macros/plotAnalysisTaskITSTPCalignment.C @@ -0,0 +1,674 @@ +void plotAnalysisTaskITSTPCalignment(const char* option = "b") +{ + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetCanvasColor(0); + gStyle->SetHistFillColor(17); + + TString optionstr(option); + gROOT->LoadMacro("AliRelAlignerKalmanArray.cxx++"); + TFile f("outputITSTPCalignment.root"); + + ////////////////////////////////////////////////////////////////////////////////////////////////// + AliRelAlignerKalmanArray* fArray = dynamic_cast(f.Get("outputArrayITSsa")); + if (!fArray) + { + printf("fArray cannot be read!\n"); + return; + } + else + { + TCanvas* c1 = new TCanvas("c1","psi in time"); + fArray->MakeGraph(0)->Draw("A*"); + c1->SaveAs("graphpsi.eps"); + + TCanvas* c2 = new TCanvas("c2","theta in time"); + fArray->MakeGraph(1)->Draw("A*"); + c2->SaveAs("graphtheta.eps"); + + TCanvas* c3 = new TCanvas("c3","phi in time"); + fArray->MakeGraph(2)->Draw("A*"); + c3->SaveAs("graphphi.eps"); + + TCanvas* c4 = new TCanvas("c4","x in time"); + fArray->MakeGraph(3)->Draw("A*"); + c4->SaveAs("graphx.eps"); + + TCanvas* c5 = new TCanvas("c5","y in time"); + fArray->MakeGraph(4)->Draw("A*"); + c5->SaveAs("graphy.eps"); + + TCanvas* c6 = new TCanvas("c6","z in time"); + fArray->MakeGraph(5)->Draw("A*"); + c6->SaveAs("graphz.eps"); + + TCanvas* c7 = new TCanvas("c7","TPC vd correction in time"); + fArray->MakeGraph(6)->Draw("A*"); + c7->SaveAs("graphvd.eps"); + + TCanvas* c8 = new TCanvas("c8","TPC t0 correction in time"); + fArray->MakeGraph(7)->Draw("A*"); + c8->SaveAs("grapht0.eps"); + + TCanvas* c9 = new TCanvas("c9","TPC dv/dy in time"); + fArray->MakeGraph(8)->Draw("A*"); + c9->SaveAs("graphdvdy.eps"); + } + + ////////////////////////////////////////////////////////////////////////////////////////////////// + TList* fList = dynamic_cast(f.Get("outputList")); + TList* pList = NULL; + + TH1F* pMatchingEfficiency = dynamic_cast(fList->At(3)); + TCanvas* canvasMEff = new TCanvas(); + pMatchingEfficiency->DrawCopy(); + + + //------------------------------------------------------------------------------------ + pList = dynamic_cast(fList->At(0)); + TH2F* pZYAResidualsHistBpos = dynamic_cast(pList->At(0)); + TH2F* pZYCResidualsHistBpos = dynamic_cast(pList->At(1)); + TH2F* pLPAResidualsHistBpos = dynamic_cast(pList->At(2)); + TH2F* pLPCResidualsHistBpos = dynamic_cast(pList->At(3)); + TH2F* pPhiYAResidualsHistBpos = dynamic_cast(pList->At(4)); + TH2F* pPhiYCResidualsHistBpos = dynamic_cast(pList->At(5)); + TH2F* pPhiZAResidualsHistBpos = dynamic_cast(pList->At(6)); + TH2F* pPhiZCResidualsHistBpos = dynamic_cast(pList->At(7)); + TH2F* pPtYAResidualsHistBpos = dynamic_cast(pList->At(8)); + TH2F* pPtYCResidualsHistBpos = dynamic_cast(pList->At(9)); + TH2F* pPtZAResidualsHistBpos = dynamic_cast(pList->At(10)); + TH2F* pPtZCResidualsHistBpos = dynamic_cast(pList->At(11)); + TH2F* pLowPtYAResidualsHistBpos = dynamic_cast(pList->At(12)); + TH2F* pLowPtYCResidualsHistBpos = dynamic_cast(pList->At(13)); + TH2F* pLowPtZAResidualsHistBpos = dynamic_cast(pList->At(14)); + TH2F* pLowPtZCResidualsHistBpos = dynamic_cast(pList->At(15)); + + TH1D* resycBpos = pZYCResidualsHistBpos->ProjectionY(); + resycBpos->SetTitle("r\\phi residual distribution side C (z<0), Bpos"); + TH1D* reszcBpos = pZYCResidualsHistBpos->ProjectionX(); + reszcBpos->SetTitle("z residual distribution side C (z<0), Bpos"); + TH1D* respcBpos = pLPCResidualsHistBpos->ProjectionY(); + respcBpos->SetTitle("sin(\\phi) residual distribution side C (z<0), Bpos"); + TH1D* reslcBpos = pLPCResidualsHistBpos->ProjectionX(); + reslcBpos->SetTitle("tan(\\lambda) residual distribution side C (z<0), Bpos"); + TProfile* profphiycBpos = pPhiYCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profphiycBpos->SetTitle("\\phi profile of r\\phi residuals side C (z<0), Bpos"); + profphiycBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizcBpos = pPhiZCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profphizcBpos->SetTitle("\\phi profile of z residuals side C (z<0), Bpos"); + profphizcBpos->SetYTitle("\\deltaz [cm]"); + TProfile* profptycBpos = pPtYCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profptycBpos->SetTitle("pt profile of r\\phi residuals side C (z<0), Bpos"); + profptycBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzcBpos = pPtZCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profptzcBpos->SetTitle("pt profile of z residuals side C (z<0), Bpos"); + profptzcBpos->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptycBpos = pLowPtYCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + proflowptycBpos->SetTitle("pt profile of r\\phi residuals side C (z<0), Bpos"); + proflowptycBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzcBpos = pLowPtZCResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + proflowptzcBpos->SetTitle("pt profile of z residuals side C (z<0), Bpos"); + proflowptzcBpos->SetYTitle("\\deltaz [cm]"); + + TH1D* resyaBpos = pZYAResidualsHistBpos->ProjectionY(); + resyaBpos->SetTitle("r\\phi residual distribution side A (z>0), Bpos"); + TH1D* reszaBpos = pZYAResidualsHistBpos->ProjectionX(); + reszaBpos->SetTitle("z residual distribution side A (z>0), Bpos"); + TH1D* respaBpos = pLPAResidualsHistBpos->ProjectionY(); + respaBpos->SetTitle("sin(\\phi) residual distribution side A (z>0), Bpos"); + TH1D* reslaBpos = pLPAResidualsHistBpos->ProjectionX(); + reslaBpos->SetTitle("tan(\\lambda) residual distribution side A (z>0), Bpos"); + TProfile* profphiyaBpos = pPhiYAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profphiyaBpos->SetTitle("\\phi profile of r\\phi residuals side A (z>0), Bpos"); + profphiyaBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizaBpos = pPhiZAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profphizaBpos->SetTitle("\\phi profile of z residuals side A (z>0), Bpos"); + profphizaBpos->SetYTitle("\\deltaz [cm]"); + TProfile* profptyaBpos = pPtYAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profptyaBpos->SetTitle("pt profile of r\\phi residuals side A (z>0), Bpos"); + profptyaBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzaBpos = pPtZAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + profptzaBpos->SetTitle("pt profile of z residuals side A (z>0), Bpos"); + profptzaBpos->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptyaBpos = pLowPtYAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + proflowptyaBpos->SetTitle("pt profile of r\\phi residuals side A (z>0), Bpos"); + proflowptyaBpos->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzaBpos = pLowPtZAResidualsHistBpos->ProfileX("_pfx",1,-1,"s"); + proflowptzaBpos->SetTitle("pt profile of z residuals side A (z>0), Bpos"); + proflowptzaBpos->SetYTitle("\\deltaz [cm]"); + + if (pZYAResidualsHistBpos->GetEntries()>0) + { + TCanvas* c1Bpos = new TCanvas("c1Bpos","Residuals ZY 2D, Bpos",800,300); + c1Bpos->Divide(2,1); + c1Bpos->cd(1); pZYAResidualsHistBpos->DrawCopy("col"); + c1Bpos->cd(2); pZYCResidualsHistBpos->DrawCopy("col"); + c1Bpos->cd(0); + c1Bpos->SaveAs("residualsZY2D-Bpos.eps"); + + TCanvas* c2Bpos = new TCanvas("c2Bpos","Residuals LP 2D, Bpos",800,300); + c2Bpos->Divide(2,1); + c2Bpos->cd(1); pLPAResidualsHistBpos->DrawCopy("col"); + c2Bpos->cd(2); pLPCResidualsHistBpos->DrawCopy("col"); + c2Bpos->cd(0); + c2Bpos->SaveAs("residualsLP2D-Bpos.eps"); + + TCanvas* c3Bpos = new TCanvas("c3Bpos","Residuals z 1D, Bpos",800,300); + c3Bpos->Divide(2,1); + c3Bpos->cd(1); reszaBpos->DrawCopy(); + c3Bpos->cd(2); reszcBpos->DrawCopy(); + c3Bpos->cd(0); + c3Bpos->SaveAs("residualsZ-Bpos.eps"); + + TCanvas* c4Bpos = new TCanvas("c4Bpos","Residuals y 1D, Bpos",800,300); + c4Bpos->Divide(2,1); + c4Bpos->cd(1); resyaBpos->DrawCopy(); + c4Bpos->cd(2); resycBpos->DrawCopy(); + c4Bpos->cd(0); + c4Bpos->SaveAs("residualsY-Bpos.eps"); + + TCanvas* c5Bpos = new TCanvas("c5Bpos","Residuals phi 1D, Bpos",800,300); + c5Bpos->Divide(2,1); + c5Bpos->cd(1); respaBpos->DrawCopy(); + c5Bpos->cd(2); respcBpos->DrawCopy(); + c5Bpos->cd(0); + c5Bpos->SaveAs("residualsPhi-Bpos.eps"); + + TCanvas* c6Bpos = new TCanvas("c6Bpos","Residuals lambda 1D, Bpos",800,300); + c6Bpos->Divide(2,1); + c6Bpos->cd(1); reslaBpos->DrawCopy(); + c6Bpos->cd(2); reslcBpos->DrawCopy(); + c6Bpos->cd(0); + c6Bpos->SaveAs("residualsLambda-Bpos.eps"); + + TCanvas* c7Bpos = new TCanvas("c7Bpos","Profiles z in phi, Bpos",800,300); + c7Bpos->Divide(2,1); + c7Bpos->cd(1); profphizaBpos->DrawCopy(); + c7Bpos->cd(2); profphizcBpos->DrawCopy(); + c7Bpos->cd(0); + c7Bpos->SaveAs("residualsProfilePhiZ-Bpos.eps"); + + TCanvas* c8Bpos = new TCanvas("c8Bpos","Profiles y in phi, Bpos",800,300); + c8Bpos->Divide(2,1); + c8Bpos->cd(1); profphiyaBpos->DrawCopy(); + c8Bpos->cd(2); profphiycBpos->DrawCopy(); + c8Bpos->cd(0); + c8Bpos->SaveAs("residualsProfilePhiY-Bpos.eps"); + + TCanvas* c9Bpos = new TCanvas("c9Bpos", "Residuals Phi-Z 2D, Bpos",800,300); + c9Bpos->Divide(2,1); + c9Bpos->cd(1); pPhiZAResidualsHistBpos->DrawCopy("col"); + c9Bpos->cd(2); pPhiZCResidualsHistBpos->DrawCopy("col"); + c9Bpos->cd(0); + c9Bpos->SaveAs("residualsPhiZ2D-Bpos.eps"); + + TCanvas* c1Bpos0 = new TCanvas("c1Bpos0", "Residuals Phi-Y 2D, Bpos",800,300); + c1Bpos0->Divide(2,1); + c1Bpos0->cd(1); pPhiYAResidualsHistBpos->DrawCopy("col"); + c1Bpos0->cd(2); pPhiYCResidualsHistBpos->DrawCopy("col"); + c1Bpos0->cd(0); + c1Bpos0->SaveAs("residualsPhiY2D-Bpos.eps"); + + TCanvas* c11Bpos = new TCanvas("c11Bpos", "Residuals Pt-Z 2D, Bpos",800,300); + c11Bpos->Divide(2,1); + c11Bpos->cd(1); pPtZAResidualsHistBpos->DrawCopy("col"); + c11Bpos->cd(2); pPtZCResidualsHistBpos->DrawCopy("col"); + c11Bpos->cd(0); + c11Bpos->SaveAs("residualsPtZ2D-Bpos.eps"); + + TCanvas* c12Bpos = new TCanvas("c12Bpos", "Residuals Pt-Y 2D, Bpos",800,300); + c12Bpos->Divide(2,1); + c12Bpos->cd(1); pPtYAResidualsHistBpos->DrawCopy("col"); + c12Bpos->cd(2); pPtYCResidualsHistBpos->DrawCopy("col"); + c12Bpos->cd(0); + c12Bpos->SaveAs("residualsPtY2D-Bpos.eps"); + + TCanvas* c13Bpos = new TCanvas("c13Bpos","Profiles Pt-Z, Bpos",800,300); + c13Bpos->Divide(2,1); + c13Bpos->cd(1); profptzaBpos->DrawCopy(); + c13Bpos->cd(2); profptzcBpos->DrawCopy(); + c13Bpos->cd(0); + c13Bpos->SaveAs("residualsProfilePtZ-Bpos.eps"); + + TCanvas* c14Bpos = new TCanvas("c14Bpos","Profiles Pt-Y, Bpos",800,300); + c14Bpos->Divide(2,1); + c14Bpos->cd(1); profptyaBpos->DrawCopy(); + c14Bpos->cd(2); profptycBpos->DrawCopy(); + c14Bpos->cd(0); + c14Bpos->SaveAs("residualsProfilePtY-Bpos.eps"); + + TCanvas* c15Bpos = new TCanvas("c15Bpos", "Residuals low Pt-Z 2D, Bpos",800,300); + c15Bpos->Divide(2,1); + c15Bpos->cd(1); pLowPtZAResidualsHistBpos->DrawCopy("col"); + c15Bpos->cd(2); pLowPtZCResidualsHistBpos->DrawCopy("col"); + c15Bpos->cd(0); + c15Bpos->SaveAs("residualsLowPtZ2D-Bpos.eps"); + + TCanvas* c16Bpos = new TCanvas("c16Bpos", "Residuals in low PtY 2D, Bpos",800,300); + c16Bpos->Divide(2,1); + c16Bpos->cd(1); pLowPtYAResidualsHistBpos->DrawCopy("col"); + c16Bpos->cd(2); pLowPtYCResidualsHistBpos->DrawCopy("col"); + c16Bpos->cd(0); + c16Bpos->SaveAs("residualsLowPtY2D-Bpos.eps"); + + TCanvas* c17Bpos = new TCanvas("c17Bpos","Profiles low Pt-Z, Bpos",800,300); + c17Bpos->Divide(2,1); + c17Bpos->cd(1); proflowptzaBpos->DrawCopy(); + c17Bpos->cd(2); proflowptzcBpos->DrawCopy(); + c17Bpos->cd(0); + c17Bpos->SaveAs("residualsProfileLowPtZ-Bpos.eps"); + + TCanvas* c18Bpos = new TCanvas("c18Bpos","Profiles low Pt-Y, Bpos",800,300); + c18Bpos->Divide(2,1); + c18Bpos->cd(1); proflowptyaBpos->DrawCopy(); + c18Bpos->cd(2); proflowptycBpos->DrawCopy(); + c18Bpos->cd(0); + c18Bpos->SaveAs("residualsProfileLowPtY-Bpos.eps"); + } + + //------------------------------------------------------------------------------------ + pList = dynamic_cast(fList->At(1)); + TH2F* pZYAResidualsHistBneg = dynamic_cast(pList->At(0)); + TH2F* pZYCResidualsHistBneg = dynamic_cast(pList->At(1)); + TH2F* pLPAResidualsHistBneg = dynamic_cast(pList->At(2)); + TH2F* pLPCResidualsHistBneg = dynamic_cast(pList->At(3)); + TH2F* pPhiYAResidualsHistBneg = dynamic_cast(pList->At(4)); + TH2F* pPhiYCResidualsHistBneg = dynamic_cast(pList->At(5)); + TH2F* pPhiZAResidualsHistBneg = dynamic_cast(pList->At(6)); + TH2F* pPhiZCResidualsHistBneg = dynamic_cast(pList->At(7)); + TH2F* pPtYAResidualsHistBneg = dynamic_cast(pList->At(8)); + TH2F* pPtYCResidualsHistBneg = dynamic_cast(pList->At(9)); + TH2F* pPtZAResidualsHistBneg = dynamic_cast(pList->At(10)); + TH2F* pPtZCResidualsHistBneg = dynamic_cast(pList->At(11)); + TH2F* pLowPtYAResidualsHistBneg = dynamic_cast(pList->At(12)); + TH2F* pLowPtYCResidualsHistBneg = dynamic_cast(pList->At(13)); + TH2F* pLowPtZAResidualsHistBneg = dynamic_cast(pList->At(14)); + TH2F* pLowPtZCResidualsHistBneg = dynamic_cast(pList->At(15)); + + TH1D* resycBneg = pZYCResidualsHistBneg->ProjectionY(); + resycBneg->SetTitle("r\\phi residual distribution side C (z<0), Bneg"); + TH1D* reszcBneg = pZYCResidualsHistBneg->ProjectionX(); + reszcBneg->SetTitle("z residual distribution side C (z<0), Bneg"); + TH1D* respcBneg = pLPCResidualsHistBneg->ProjectionY(); + respcBneg->SetTitle("sin(\\phi) residual distribution side C (z<0), Bneg"); + TH1D* reslcBneg = pLPCResidualsHistBneg->ProjectionX(); + reslcBneg->SetTitle("tan(\\lambda) residual distribution side C (z<0), Bneg"); + TProfile* profphiycBneg = pPhiYCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profphiycBneg->SetTitle("\\phi profile of r\\phi residuals side C (z<0), Bneg"); + profphiycBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizcBneg = pPhiZCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profphizcBneg->SetTitle("\\phi profile of z residuals side C (z<0), Bneg"); + profphizcBneg->SetYTitle("\\deltaz [cm]"); + TProfile* profptycBneg = pPtYCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profptycBneg->SetTitle("pt profile of r\\phi residuals side C (z<0), Bneg"); + profptycBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzcBneg = pPtZCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profptzcBneg->SetTitle("pt profile of z residuals side C (z<0), Bneg"); + profptzcBneg->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptycBneg = pLowPtYCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + proflowptycBneg->SetTitle("pt profile of r\\phi residuals side C (z<0), Bneg"); + proflowptycBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzcBneg = pLowPtZCResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + proflowptzcBneg->SetTitle("pt profile of z residuals side C (z<0), Bneg"); + proflowptzcBneg->SetYTitle("\\deltaz [cm]"); + + TH1D* resyaBneg = pZYAResidualsHistBneg->ProjectionY(); + resyaBneg->SetTitle("r\\phi residual distribution side A (z>0), Bneg"); + TH1D* reszaBneg = pZYAResidualsHistBneg->ProjectionX(); + reszaBneg->SetTitle("z residual distribution side A (z>0), Bneg"); + TH1D* respaBneg = pLPAResidualsHistBneg->ProjectionY(); + respaBneg->SetTitle("sin(\\phi) residual distribution side A (z>0), Bneg"); + TH1D* reslaBneg = pLPAResidualsHistBneg->ProjectionX(); + reslaBneg->SetTitle("tan(\\lambda) residual distribution side A (z>0), Bneg"); + TProfile* profphiyaBneg = pPhiYAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profphiyaBneg->SetTitle("\\phi profile of r\\phi residuals side A (z>0), Bneg"); + profphiyaBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizaBneg = pPhiZAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profphizaBneg->SetTitle("\\phi profile of z residuals side A (z>0), Bneg"); + profphizaBneg->SetYTitle("\\deltaz [cm]"); + TProfile* profptyaBneg = pPtYAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profptyaBneg->SetTitle("pt profile of r\\phi residuals side A (z>0), Bneg"); + profptyaBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzaBneg = pPtZAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + profptzaBneg->SetTitle("pt profile of z residuals side A (z>0), Bneg"); + profptzaBneg->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptyaBneg = pLowPtYAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + proflowptyaBneg->SetTitle("pt profile of r\\phi residuals side A (z>0), Bneg"); + proflowptyaBneg->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzaBneg = pLowPtZAResidualsHistBneg->ProfileX("_pfx",1,-1,"s"); + proflowptzaBneg->SetTitle("pt profile of z residuals side A (z>0), Bneg"); + proflowptzaBneg->SetYTitle("\\deltaz [cm]"); + + if (pZYAResidualsHistBneg->GetEntries()>0) + { + TCanvas* c1Bneg = new TCanvas("c1Bneg","Residuals ZY 2D, Bneg",800,300); + c1Bneg->Divide(2,1); + c1Bneg->cd(1); pZYAResidualsHistBneg->DrawCopy("col"); + c1Bneg->cd(2); pZYCResidualsHistBneg->DrawCopy("col"); + c1Bneg->cd(0); + c1Bneg->SaveAs("residualsZY2D-Bneg.eps"); + + TCanvas* c2Bneg = new TCanvas("c2Bneg","Residuals LP 2D, Bneg",800,300); + c2Bneg->Divide(2,1); + c2Bneg->cd(1); pLPAResidualsHistBneg->DrawCopy("col"); + c2Bneg->cd(2); pLPCResidualsHistBneg->DrawCopy("col"); + c2Bneg->cd(0); + c2Bneg->SaveAs("residualsLP2D-Bneg.eps"); + + TCanvas* c3Bneg = new TCanvas("c3Bneg","Residuals z 1D, Bneg",800,300); + c3Bneg->Divide(2,1); + c3Bneg->cd(1); reszaBneg->DrawCopy(); + c3Bneg->cd(2); reszcBneg->DrawCopy(); + c3Bneg->cd(0); + c3Bneg->SaveAs("residualsZ-Bneg.eps"); + + TCanvas* c4Bneg = new TCanvas("c4Bneg","Residuals y 1D, Bneg",800,300); + c4Bneg->Divide(2,1); + c4Bneg->cd(1); resyaBneg->DrawCopy(); + c4Bneg->cd(2); resycBneg->DrawCopy(); + c4Bneg->cd(0); + c4Bneg->SaveAs("residualsY-Bneg.eps"); + + TCanvas* c5Bneg = new TCanvas("c5Bneg","Residuals phi 1D, Bneg",800,300); + c5Bneg->Divide(2,1); + c5Bneg->cd(1); respaBneg->DrawCopy(); + c5Bneg->cd(2); respcBneg->DrawCopy(); + c5Bneg->cd(0); + c5Bneg->SaveAs("residualsPhi-Bneg.eps"); + + TCanvas* c6Bneg = new TCanvas("c6Bneg","Residuals lambda 1D, Bneg",800,300); + c6Bneg->Divide(2,1); + c6Bneg->cd(1); reslaBneg->DrawCopy(); + c6Bneg->cd(2); reslcBneg->DrawCopy(); + c6Bneg->cd(0); + c6Bneg->SaveAs("residualsLambda-Bneg.eps"); + + TCanvas* c7Bneg = new TCanvas("c7Bneg","Profiles z in phi, Bneg",800,300); + c7Bneg->Divide(2,1); + c7Bneg->cd(1); profphizaBneg->DrawCopy(); + c7Bneg->cd(2); profphizcBneg->DrawCopy(); + c7Bneg->cd(0); + c7Bneg->SaveAs("residualsProfilePhiZ-Bneg.eps"); + + TCanvas* c8Bneg = new TCanvas("c8Bneg","Profiles y in phi, Bneg",800,300); + c8Bneg->Divide(2,1); + c8Bneg->cd(1); profphiyaBneg->DrawCopy(); + c8Bneg->cd(2); profphiycBneg->DrawCopy(); + c8Bneg->cd(0); + c8Bneg->SaveAs("residualsProfilePhiY-Bneg.eps"); + + TCanvas* c9Bneg = new TCanvas("c9Bneg", "Residuals Phi-Z 2D, Bneg",800,300); + c9Bneg->Divide(2,1); + c9Bneg->cd(1); pPhiZAResidualsHistBneg->DrawCopy("col"); + c9Bneg->cd(2); pPhiZCResidualsHistBneg->DrawCopy("col"); + c9Bneg->cd(0); + c9Bneg->SaveAs("residualsPhiZ2D-Bneg.eps"); + + TCanvas* c1Bneg0 = new TCanvas("c1Bneg0", "Residuals Phi-Y 2D, Bneg",800,300); + c1Bneg0->Divide(2,1); + c1Bneg0->cd(1); pPhiYAResidualsHistBneg->DrawCopy("col"); + c1Bneg0->cd(2); pPhiYCResidualsHistBneg->DrawCopy("col"); + c1Bneg0->cd(0); + c1Bneg0->SaveAs("residualsPhiY2D-Bneg.eps"); + + TCanvas* c11Bneg = new TCanvas("c11Bneg", "Residuals Pt-Z 2D, Bneg",800,300); + c11Bneg->Divide(2,1); + c11Bneg->cd(1); pPtZAResidualsHistBneg->DrawCopy("col"); + c11Bneg->cd(2); pPtZCResidualsHistBneg->DrawCopy("col"); + c11Bneg->cd(0); + c11Bneg->SaveAs("residualsPtZ2D-Bneg.eps"); + + TCanvas* c12Bneg = new TCanvas("c12Bneg", "Residuals Pt-Y 2D, Bneg",800,300); + c12Bneg->Divide(2,1); + c12Bneg->cd(1); pPtYAResidualsHistBneg->DrawCopy("col"); + c12Bneg->cd(2); pPtYCResidualsHistBneg->DrawCopy("col"); + c12Bneg->cd(0); + c12Bneg->SaveAs("residualsPtY2D-Bneg.eps"); + + TCanvas* c13Bneg = new TCanvas("c13Bneg","Profiles Pt-Z, Bneg",800,300); + c13Bneg->Divide(2,1); + c13Bneg->cd(1); profptzaBneg->DrawCopy(); + c13Bneg->cd(2); profptzcBneg->DrawCopy(); + c13Bneg->cd(0); + c13Bneg->SaveAs("residualsProfilePtZ-Bneg.eps"); + + TCanvas* c14Bneg = new TCanvas("c14Bneg","Profiles Pt-Y, Bneg",800,300); + c14Bneg->Divide(2,1); + c14Bneg->cd(1); profptyaBneg->DrawCopy(); + c14Bneg->cd(2); profptycBneg->DrawCopy(); + c14Bneg->cd(0); + c14Bneg->SaveAs("residualsProfilePtY-Bneg.eps"); + + TCanvas* c15Bneg = new TCanvas("c15Bneg", "Residuals low Pt-Z 2D, Bneg",800,300); + c15Bneg->Divide(2,1); + c15Bneg->cd(1); pLowPtZAResidualsHistBneg->DrawCopy("col"); + c15Bneg->cd(2); pLowPtZCResidualsHistBneg->DrawCopy("col"); + c15Bneg->cd(0); + c15Bneg->SaveAs("residualsLowPtZ2D-Bneg.eps"); + + TCanvas* c16Bneg = new TCanvas("c16Bneg", "Residuals in low PtY 2D, Bneg",800,300); + c16Bneg->Divide(2,1); + c16Bneg->cd(1); pLowPtYAResidualsHistBneg->DrawCopy("col"); + c16Bneg->cd(2); pLowPtYCResidualsHistBneg->DrawCopy("col"); + c16Bneg->cd(0); + c16Bneg->SaveAs("residualsLowPtY2D-Bneg.eps"); + + TCanvas* c17Bneg = new TCanvas("c17Bneg","Profiles low Pt-Z, Bneg",800,300); + c17Bneg->Divide(2,1); + c17Bneg->cd(1); proflowptzaBneg->DrawCopy(); + c17Bneg->cd(2); proflowptzcBneg->DrawCopy(); + c17Bneg->cd(0); + c17Bneg->SaveAs("residualsProfileLowPtZ-Bneg.eps"); + + TCanvas* c18Bneg = new TCanvas("c18Bneg","Profiles low Pt-Y, Bneg",800,300); + c18Bneg->Divide(2,1); + c18Bneg->cd(1); proflowptyaBneg->DrawCopy(); + c18Bneg->cd(2); proflowptycBneg->DrawCopy(); + c18Bneg->cd(0); + c18Bneg->SaveAs("residualsProfileLowPtY-Bneg.eps"); + } + + //------------------------------------------------------------------------------------ + pList = dynamic_cast(fList->At(2)); + TH2F* pZYAResidualsHistBnil = dynamic_cast(pList->At(0)); + TH2F* pZYCResidualsHistBnil = dynamic_cast(pList->At(1)); + TH2F* pLPAResidualsHistBnil = dynamic_cast(pList->At(2)); + TH2F* pLPCResidualsHistBnil = dynamic_cast(pList->At(3)); + TH2F* pPhiYAResidualsHistBnil = dynamic_cast(pList->At(4)); + TH2F* pPhiYCResidualsHistBnil = dynamic_cast(pList->At(5)); + TH2F* pPhiZAResidualsHistBnil = dynamic_cast(pList->At(6)); + TH2F* pPhiZCResidualsHistBnil = dynamic_cast(pList->At(7)); + TH2F* pPtYAResidualsHistBnil = dynamic_cast(pList->At(8)); + TH2F* pPtYCResidualsHistBnil = dynamic_cast(pList->At(9)); + TH2F* pPtZAResidualsHistBnil = dynamic_cast(pList->At(10)); + TH2F* pPtZCResidualsHistBnil = dynamic_cast(pList->At(11)); + TH2F* pLowPtYAResidualsHistBnil = dynamic_cast(pList->At(12)); + TH2F* pLowPtYCResidualsHistBnil = dynamic_cast(pList->At(13)); + TH2F* pLowPtZAResidualsHistBnil = dynamic_cast(pList->At(14)); + TH2F* pLowPtZCResidualsHistBnil = dynamic_cast(pList->At(15)); + + TH1D* resycBnil = pZYCResidualsHistBnil->ProjectionY(); + resycBnil->SetTitle("r\\phi residual distribution side C (z<0), Bnil"); + TH1D* reszcBnil = pZYCResidualsHistBnil->ProjectionX(); + reszcBnil->SetTitle("z residual distribution side C (z<0), Bnil"); + TH1D* respcBnil = pLPCResidualsHistBnil->ProjectionY(); + respcBnil->SetTitle("sin(\\phi) residual distribution side C (z<0), Bnil"); + TH1D* reslcBnil = pLPCResidualsHistBnil->ProjectionX(); + reslcBnil->SetTitle("tan(\\lambda) residual distribution side C (z<0), Bnil"); + TProfile* profphiycBnil = pPhiYCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profphiycBnil->SetTitle("\\phi profile of r\\phi residuals side C (z<0), Bnil"); + profphiycBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizcBnil = pPhiZCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profphizcBnil->SetTitle("\\phi profile of z residuals side C (z<0), Bnil"); + profphizcBnil->SetYTitle("\\deltaz [cm]"); + TProfile* profptycBnil = pPtYCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profptycBnil->SetTitle("pt profile of r\\phi residuals side C (z<0), Bnil"); + profptycBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzcBnil = pPtZCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profptzcBnil->SetTitle("pt profile of z residuals side C (z<0), Bnil"); + profptzcBnil->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptycBnil = pLowPtYCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + proflowptycBnil->SetTitle("pt profile of r\\phi residuals side C (z<0), Bnil"); + proflowptycBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzcBnil = pLowPtZCResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + proflowptzcBnil->SetTitle("pt profile of z residuals side C (z<0), Bnil"); + proflowptzcBnil->SetYTitle("\\deltaz [cm]"); + + TH1D* resyaBnil = pZYAResidualsHistBnil->ProjectionY(); + resyaBnil->SetTitle("r\\phi residual distribution side A (z>0), Bnil"); + TH1D* reszaBnil = pZYAResidualsHistBnil->ProjectionX(); + reszaBnil->SetTitle("z residual distribution side A (z>0), Bnil"); + TH1D* respaBnil = pLPAResidualsHistBnil->ProjectionY(); + respaBnil->SetTitle("sin(\\phi) residual distribution side A (z>0), Bnil"); + TH1D* reslaBnil = pLPAResidualsHistBnil->ProjectionX(); + reslaBnil->SetTitle("tan(\\lambda) residual distribution side A (z>0), Bnil"); + TProfile* profphiyaBnil = pPhiYAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profphiyaBnil->SetTitle("\\phi profile of r\\phi residuals side A (z>0), Bnil"); + profphiyaBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profphizaBnil = pPhiZAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profphizaBnil->SetTitle("\\phi profile of z residuals side A (z>0), Bnil"); + profphizaBnil->SetYTitle("\\deltaz [cm]"); + TProfile* profptyaBnil = pPtYAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profptyaBnil->SetTitle("pt profile of r\\phi residuals side A (z>0), Bnil"); + profptyaBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* profptzaBnil = pPtZAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + profptzaBnil->SetTitle("pt profile of z residuals side A (z>0), Bnil"); + profptzaBnil->SetYTitle("\\deltaz [cm]"); + TProfile* proflowptyaBnil = pLowPtYAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + proflowptyaBnil->SetTitle("pt profile of r\\phi residuals side A (z>0), Bnil"); + proflowptyaBnil->SetYTitle("\\deltar\\phi [cm]"); + TProfile* proflowptzaBnil = pLowPtZAResidualsHistBnil->ProfileX("_pfx",1,-1,"s"); + proflowptzaBnil->SetTitle("pt profile of z residuals side A (z>0), Bnil"); + proflowptzaBnil->SetYTitle("\\deltaz [cm]"); + + if (pZYAResidualsHistBnil->GetEntries()>0) + { + TCanvas* c1Bnil = new TCanvas("c1Bnil","Residuals ZY 2D, Bnil",800,300); + c1Bnil->Divide(2,1); + c1Bnil->cd(1); pZYAResidualsHistBnil->DrawCopy("col"); + c1Bnil->cd(2); pZYCResidualsHistBnil->DrawCopy("col"); + c1Bnil->cd(0); + c1Bnil->SaveAs("residualsZY2D-Bnil.eps"); + + TCanvas* c2Bnil = new TCanvas("c2Bnil","Residuals LP 2D, Bnil",800,300); + c2Bnil->Divide(2,1); + c2Bnil->cd(1); pLPAResidualsHistBnil->DrawCopy("col"); + c2Bnil->cd(2); pLPCResidualsHistBnil->DrawCopy("col"); + c2Bnil->cd(0); + c2Bnil->SaveAs("residualsLP2D-Bnil.eps"); + + TCanvas* c3Bnil = new TCanvas("c3Bnil","Residuals z 1D, Bnil",800,300); + c3Bnil->Divide(2,1); + c3Bnil->cd(1); reszaBnil->DrawCopy(); + c3Bnil->cd(2); reszcBnil->DrawCopy(); + c3Bnil->cd(0); + c3Bnil->SaveAs("residualsZ-Bnil.eps"); + + TCanvas* c4Bnil = new TCanvas("c4Bnil","Residuals y 1D, Bnil",800,300); + c4Bnil->Divide(2,1); + c4Bnil->cd(1); resyaBnil->DrawCopy(); + c4Bnil->cd(2); resycBnil->DrawCopy(); + c4Bnil->cd(0); + c4Bnil->SaveAs("residualsY-Bnil.eps"); + + TCanvas* c5Bnil = new TCanvas("c5Bnil","Residuals phi 1D, Bnil",800,300); + c5Bnil->Divide(2,1); + c5Bnil->cd(1); respaBnil->DrawCopy(); + c5Bnil->cd(2); respcBnil->DrawCopy(); + c5Bnil->cd(0); + c5Bnil->SaveAs("residualsPhi-Bnil.eps"); + + TCanvas* c6Bnil = new TCanvas("c6Bnil","Residuals lambda 1D, Bnil",800,300); + c6Bnil->Divide(2,1); + c6Bnil->cd(1); reslaBnil->DrawCopy(); + c6Bnil->cd(2); reslcBnil->DrawCopy(); + c6Bnil->cd(0); + c6Bnil->SaveAs("residualsLambda-Bnil.eps"); + + TCanvas* c7Bnil = new TCanvas("c7Bnil","Profiles z in phi, Bnil",800,300); + c7Bnil->Divide(2,1); + c7Bnil->cd(1); profphizaBnil->DrawCopy(); + c7Bnil->cd(2); profphizcBnil->DrawCopy(); + c7Bnil->cd(0); + c7Bnil->SaveAs("residualsProfilePhiZ-Bnil.eps"); + + TCanvas* c8Bnil = new TCanvas("c8Bnil","Profiles y in phi, Bnil",800,300); + c8Bnil->Divide(2,1); + c8Bnil->cd(1); profphiyaBnil->DrawCopy(); + c8Bnil->cd(2); profphiycBnil->DrawCopy(); + c8Bnil->cd(0); + c8Bnil->SaveAs("residualsProfilePhiY-Bnil.eps"); + + TCanvas* c9Bnil = new TCanvas("c9Bnil", "Residuals Phi-Z 2D, Bnil",800,300); + c9Bnil->Divide(2,1); + c9Bnil->cd(1); pPhiZAResidualsHistBnil->DrawCopy("col"); + c9Bnil->cd(2); pPhiZCResidualsHistBnil->DrawCopy("col"); + c9Bnil->cd(0); + c9Bnil->SaveAs("residualsPhiZ2D-Bnil.eps"); + + TCanvas* c1Bnil0 = new TCanvas("c1Bnil0", "Residuals Phi-Y 2D, Bnil",800,300); + c1Bnil0->Divide(2,1); + c1Bnil0->cd(1); pPhiYAResidualsHistBnil->DrawCopy("col"); + c1Bnil0->cd(2); pPhiYCResidualsHistBnil->DrawCopy("col"); + c1Bnil0->cd(0); + c1Bnil0->SaveAs("residualsPhiY2D-Bnil.eps"); + + TCanvas* c11Bnil = new TCanvas("c11Bnil", "Residuals Pt-Z 2D, Bnil",800,300); + c11Bnil->Divide(2,1); + c11Bnil->cd(1); pPtZAResidualsHistBnil->DrawCopy("col"); + c11Bnil->cd(2); pPtZCResidualsHistBnil->DrawCopy("col"); + c11Bnil->cd(0); + c11Bnil->SaveAs("residualsPtZ2D-Bnil.eps"); + + TCanvas* c12Bnil = new TCanvas("c12Bnil", "Residuals Pt-Y 2D, Bnil",800,300); + c12Bnil->Divide(2,1); + c12Bnil->cd(1); pPtYAResidualsHistBnil->DrawCopy("col"); + c12Bnil->cd(2); pPtYCResidualsHistBnil->DrawCopy("col"); + c12Bnil->cd(0); + c12Bnil->SaveAs("residualsPtY2D-Bnil.eps"); + + TCanvas* c13Bnil = new TCanvas("c13Bnil","Profiles Pt-Z, Bnil",800,300); + c13Bnil->Divide(2,1); + c13Bnil->cd(1); profptzaBnil->DrawCopy(); + c13Bnil->cd(2); profptzcBnil->DrawCopy(); + c13Bnil->cd(0); + c13Bnil->SaveAs("residualsProfilePtZ-Bnil.eps"); + + TCanvas* c14Bnil = new TCanvas("c14Bnil","Profiles Pt-Y, Bnil",800,300); + c14Bnil->Divide(2,1); + c14Bnil->cd(1); profptyaBnil->DrawCopy(); + c14Bnil->cd(2); profptycBnil->DrawCopy(); + c14Bnil->cd(0); + c14Bnil->SaveAs("residualsProfilePtY-Bnil.eps"); + + TCanvas* c15Bnil = new TCanvas("c15Bnil", "Residuals low Pt-Z 2D, Bnil",800,300); + c15Bnil->Divide(2,1); + c15Bnil->cd(1); pLowPtZAResidualsHistBnil->DrawCopy("col"); + c15Bnil->cd(2); pLowPtZCResidualsHistBnil->DrawCopy("col"); + c15Bnil->cd(0); + c15Bnil->SaveAs("residualsLowPtZ2D-Bnil.eps"); + + TCanvas* c16Bnil = new TCanvas("c16Bnil", "Residuals in low PtY 2D, Bnil",800,300); + c16Bnil->Divide(2,1); + c16Bnil->cd(1); pLowPtYAResidualsHistBnil->DrawCopy("col"); + c16Bnil->cd(2); pLowPtYCResidualsHistBnil->DrawCopy("col"); + c16Bnil->cd(0); + c16Bnil->SaveAs("residualsLowPtY2D-Bnil.eps"); + + TCanvas* c17Bnil = new TCanvas("c17Bnil","Profiles low Pt-Z, Bnil",800,300); + c17Bnil->Divide(2,1); + c17Bnil->cd(1); proflowptzaBnil->DrawCopy(); + c17Bnil->cd(2); proflowptzcBnil->DrawCopy(); + c17Bnil->cd(0); + c17Bnil->SaveAs("residualsProfileLowPtZ-Bnil.eps"); + + TCanvas* c18Bnil = new TCanvas("c18Bnil","Profiles low Pt-Y, Bnil",800,300); + c18Bnil->Divide(2,1); + c18Bnil->cd(1); proflowptyaBnil->DrawCopy(); + c18Bnil->cd(2); proflowptycBnil->DrawCopy(); + c18Bnil->cd(0); + c18Bnil->SaveAs("residualsProfileLowPtY-Bnil.eps"); + } +} + diff --git a/PWG1/macros/runAnalysisTaskITSTPCalignment.C b/PWG1/macros/runAnalysisTaskITSTPCalignment.C new file mode 100644 index 00000000000..2a7fc9e3ea5 --- /dev/null +++ b/PWG1/macros/runAnalysisTaskITSTPCalignment.C @@ -0,0 +1,294 @@ +void runAnalysisTaskITSTPCalignment() +{ + TStopwatch timer; + timer.Start(); + + //runProof("/ITS/dainesea/run104070#esdTree"); + //runProof("/ITS/dainesea/run104070_newTPCalign#esdTree"); + //runProof("/ALIREC/aliprod/run104671#esdTree"); + //runProof("/ITS/dainesea/run104070_newTPCcalib#esdTree"); + //runProof("/COMMON/COMMON/LHC09d9a_0.9TeV_0.5T#esdTree"); + //runLocal("find /data/alice3/mikolaj/runs_pp/pass1/ -name AliESDs.root","tree"); + //runLocal("find /data/alice3/mikolaj/ITSmisal -path */0.001/* -name AliESDs.root","tree"); + //runLocal("find /data/alice3/mikolaj/TPCfullmisalignmentB0 -name AliESDs.root","tree"); + //runLocal("find /data/alice3/mikolaj/TPCfullmisalignmentB5/ -name AliESDs.root -path \"*187824*\" ","tree"); + //runLocal("find /data/alice3/mikolaj/TPCfullmisalignmentB5/ -name AliESDs.root","tree"); + //runLocal("find /data/alice3/mikolaj/run104892test -name AliESDs.root","tree"); + runAlienPlugin("terminate"); + + timer.Stop(); + timer.Print(); +} + +//______________________________________________________________________________ +void runLocal(TString findCommand = "", TString options="") +{ + + TString outputFilename = "outputITSTPCalignment.root"; + + //setupPar("STEERBase"); + //setupPar("ESD"); + //setupPar("AOD"); + //setupPar("ANALYSIS"); + //setupPar("ANALYSISalice"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); // The plugin is here + + gSystem->AddIncludePath("-I$ALICE_ROOT/include"); + gROOT->LoadMacro("AliRelAlignerKalman.cxx++"); + gROOT->LoadMacro("AliRelAlignerKalmanArray.cxx++g"); + gROOT->LoadMacro("AliAnalysisTaskITSTPCalignment.cxx++g"); + gROOT->LoadMacro("AddTaskITSTPCalignment.C"); + gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C"); + + //add input files from list saved in a file + findCommand += ">ESDfiles.txt"; + gSystem->Exec(findCommand); + TChain* chain = CreateESDChain("ESDfiles.txt",100000); + + // analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("ITSTPCalignmentAnalysisManager"); + + // input handlers + AliVEventHandler* esdH = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdH); + + //AliMCEventHandler* MCH = new AliMCEventHandler(); + //mgr->SetMCtruthEventHandler(MCH); + + //add the task + AliAnalysisTaskITSTPCalignment *task = AddTaskITSTPCalignment(); + + //start analysis + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + + mgr->StartAnalysis("local",chain); +} + +//______________________________________________________________________________ +void runProof(const char* dataset, TString options="" ) +{ + + TString outputFilename = "outputITSTPCalignment.root"; + + printf("****** Connect to PROOF *******\n"); + gEnv->SetValue("XSec.GSI.DelegProxy","2"); + TProof::Open("mkrzewic@alicecaf.cern.ch"); + //gProof->SetParallel(); + gProof->ClearPackages(); + + // Enable the Analysis Package + gProof->UploadPackage("STEERBase.par"); + gProof->EnablePackage("STEERBase"); + gProof->UploadPackage("ESD.par"); + gProof->EnablePackage("ESD"); + gProof->UploadPackage("AOD.par"); + gProof->EnablePackage("AOD"); + gProof->UploadPackage("ANALYSIS.par"); + gProof->EnablePackage("ANALYSIS"); + gProof->UploadPackage("ANALYSISalice.par"); + gProof->EnablePackage("ANALYSISalice"); + gProof->UploadPackage("CORRFW.par"); + gProof->EnablePackage("CORRFW"); + + //gProof->UploadPackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-18-12-AN/AF-v4-18-12-AN"); + //gProof->EnablePackage("/afs/cern.ch/alice/caf/sw/ALICE/PARs/v4-18-12-AN/AF-v4-18-12-AN"); + + gProof->Load("AliRelAlignerKalman.cxx++g"); + gProof->Load("AliRelAlignerKalmanArray.cxx++g"); + gProof->Load("AliAnalysisTaskITSTPCalignment.cxx++g"); + gProof->Load("AddTaskITSTPCalignment.C"); + + // analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("ITSTPCalignmentAnalysisManager"); + + // input handlers + AliVEventHandler* esdH = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdH); + + //add the task + AliAnalysisTaskITSTPCalignment *task = AddTaskITSTPCalignment(); + + //start analysis + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + + mgr->StartAnalysis("proof",dataset,10000,0); +} + +//______________________________________________________________________________ +void runAlienPlugin(const char* pluginmode="full") +{ + //must be configured in CreateAlienHandler.C file + + TString outputFilename = "outputITSTPCalignment.root"; + + // Load common libraries + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libSTEERBase"); + gSystem->Load("libESD"); + gSystem->Load("libAOD"); + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); // The plugin is here + + // Use AliRoot includes to compile our task + gSystem->AddIncludePath("-I$ALICE_ROOT/include"); + + ///////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////////////////// + // Create and configure the alien handler plugin + // Check if user has a valid token, otherwise make one. This has limitations. + // One can always follow the standard procedure of calling alien-token-init then + // source /tmp/gclient_env_$UID in the current shell. + if (!AliAnalysisGrid::CreateToken()) return NULL; + AliAnalysisAlien *plugin = new AliAnalysisAlien(); + // Set the run mode (can be "full", "test", "offline", "submit" or "terminate") + plugin->SetRunMode(pluginmode); + // Set versions of used packages + plugin->SetAPIVersion("V1.1x"); + plugin->SetROOTVersion("v5-26-00b-4"); + plugin->SetAliROOTVersion("v4-19-11-AN"); + // Declare input data to be processed. + // Method 1: Create automatically XML collections using alien 'find' command. + // Define production directory LFN + plugin->SetGridDataDir("/alice/data/2010/LHC10b"); + // Set data search pattern + plugin->SetDataPattern("*/pass1/*AliESDs.root"); + // ...then add run numbers to be considered + plugin->SetRunPrefix("000"); + plugin->AddRunNumber(117222); + plugin->AddRunNumber(117220); + plugin->AddRunNumber(117112); + plugin->AddRunNumber(117099); + plugin->AddRunNumber(117048); + plugin->AddRunNumber(117048); + plugin->AddRunNumber(116288); + plugin->AddRunNumber(115322); + plugin->AddRunNumber(114931); + //plugin->AddRunNumber(117051); + //plugin->SetRunRange(117029,117121); + // Method 2: Declare existing data files (raw collections, xml collections, root file) + // If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir()) + // XML collections added via this method can be combined with the first method if + // the content is compatible (using or not tags) + // plugin->AddDataFile("/alice/cern.ch/user/m/mkrzewic/pp2009pass2.xml"); + // plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root"); + // Define alien work directory where all files will be copied. Relative to alien $HOME. + plugin->SetGridWorkingDir("analysisGridLHC10bPass1_3runs"); + // Declare alien output directory. Relative to working directory. + plugin->SetGridOutputDir("output"); // In this case will be $HOME/work/output + // Declare the analysis source files names separated by blancs. To be compiled runtime + // using ACLiC on the worker nodes. + plugin->SetAnalysisSource("AliRelAlignerKalman.cxx AliRelAlignerKalmanArray.cxx AliAnalysisTaskITSTPCalignment.cxx"); + // Declare all libraries (other than the default ones for the framework. These will be + // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here. + plugin->SetAdditionalLibs("AliRelAlignerKalman.h AliRelAlignerKalman.cxx AliRelAlignerKalmanArray.h AliRelAlignerKalmanArray.cxx AliAnalysisTaskITSTPCalignment.h AliAnalysisTaskITSTPCalignment.cxx");// AddTaskITSTPCalignment.C"); + // Declare the output file names separated by blancs. + // (can be like: file.root or file.root@ALICE::Niham::File) + plugin->SetOutputFiles("outputITSTPCalignment.root"); + // Optionally define the files to be archived. + // plugin->SetOutputArchive("log_archive.zip:stdout,stderr@ALICE::NIHAM::File root_archive.zip:*.root@ALICE::NIHAM::File"); + plugin->SetOutputArchive("log_archive.zip:stdout,stderr"); + // Optionally set a name for the generated analysis macro (default MyAnalysis.C) + plugin->SetAnalysisMacro("AnalysisITSTPCalignmentGenerated.C"); + // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore) + plugin->SetSplitMaxInputFileNumber(20); + // Optionally set number of failed jobs that will trigger killing waiting sub-jobs. + plugin->SetMaxInitFailed(5); + // Optionally resubmit threshold. + plugin->SetMasterResubmitThreshold(90); + // Optionally set time to live (default 30000 sec) + plugin->SetTTL(100000); + // Optionally set input format (default xml-single) + plugin->SetInputFormat("xml-single"); + // Optionally modify the name of the generated JDL (default analysis.jdl) + plugin->SetJDLName("taskITSTPCalignment.jdl"); + // Optionally modify job price (default 1) + plugin->SetPrice(1); + // Optionally modify split mode (default 'se') + plugin->SetSplitMode("se"); + ///////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////////////////// + + // analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("ITSTPCalignmentAnalysisManager"); + + // input handlers + mgr->SetGridHandler(plugin); + + //input handler + AliESDInputHandler* esdH = new AliESDInputHandler(); + esdH->SetInactiveBranches("Calo FMD"); + //esdH->SetReadFriends(kTRUE); + mgr->SetInputEventHandler(esdH); + + gROOT->LoadMacro("AliRelAlignerKalman.cxx++g"); + gROOT->LoadMacro("AliRelAlignerKalmanArray.cxx++g"); + gROOT->LoadMacro("AliAnalysisTaskITSTPCalignment.cxx++g"); + gROOT->LoadMacro("AddTaskITSTPCalignment.C"); + + //add the task + AliAnalysisTaskITSTPCalignment *task = AddTaskITSTPCalignment(); + + //mgr->SetDebugLevel(2); + + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + mgr->StartAnalysis("grid"); +} + +//______________________________________________________________________________ +Int_t setupPar(const char* pararchivename) +{ + /////////////////// + // Setup PAR File// + /////////////////// + if (pararchivename) + { + char processline[1024]; + sprintf(processline,".! tar xvzf %s.par",pararchivename); + gROOT->ProcessLine(processline); + const char* ocwd = gSystem->WorkingDirectory(); + gSystem->ChangeDirectory(pararchivename); + + // check for BUILD.sh and execute + if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) + { + printf("*******************************\n"); + printf("*** Building PAR archive ***\n"); + printf("*******************************\n"); + + if (gSystem->Exec("PROOF-INF/BUILD.sh")) + { + Error("runAnalysis","Cannot Build the PAR Archive! - Abort!"); + return -1; + } + } + // check for SETUP.C and execute + if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) + { + printf("*******************************\n"); + printf("*** Setup PAR archive ***\n"); + printf("*******************************\n"); + gROOT->Macro("PROOF-INF/SETUP.C"); + } + + gSystem->ChangeDirectory("../"); + } + return 1; +} + + -- 2.43.0