PndScrutAnaTask.cxx

Go to the documentation of this file.

// ************************************************************************
//
// Template for an Analysis Task, 
// 
// for the December 2013 CM Tutorial 
//
// For further info see also
//
// http://panda-wiki.gsi.de/cgi-bin/viewauth/Computing/PandaRootRhoTutorial
// http://panda-wiki.gsi.de/cgi-bin/view/Computing/PandaRootAnalysisJuly13
//
// K.Goetzen 11/2013
//
// ************************************************************************


// The header file
#include "PndScrutAnaTask.h"

// C++ headers
#include <string>
#include <iostream>

// FAIR headers
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"

// ROOT headers
#include "TClonesArray.h"
#include "TVector3.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TDatabasePDG.h"

// RHO headers
#include "RhoCandidate.h"
#include "RhoHistogram/RhoTuple.h"
#include "RhoFactory.h"
#include "RhoMassParticleSelector.h"
#include "RhoTuple.h"

// Analysis headers
#include "PndAnalysis.h"
#include "Rho4CFitter.h"
#include "RhoKinVtxFitter.h"
#include "RhoKinFitter.h"
#include "RhoVtxPoca.h"
#include "PndRhoTupleQA.h"
#include "PndEventShape.h"
                
                
using std::cout;
using std::endl;


// -----   Default constructor   -------------------------------------------
PndScrutAnaTask::PndScrutAnaTask(double pbarmom, TString outname) :
  FairTask("Panda Scrutiny Analysis Task") 
{ 
        double mp=0.938272;
        fIni.SetXYZT(0,0,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);

        fOutName = outname;
}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
PndScrutAnaTask::~PndScrutAnaTask() 
{ 
        delete fFile;
}
// -------------------------------------------------------------------------


// -----   Public method Init   --------------------------------------------
InitStatus PndScrutAnaTask::Init() 
{               
        // *** initialize PndAnalysis object
        fAnalysis = new PndAnalysis();
        
        // *** reset the event counter
        fEvtCount = 0;

        // *******
        // ******* PREPARE/CREATE THE STUFF YOU NEED
        // *******
        
        fPdg = TDatabasePDG::Instance();

        // *** Mass selector for the jpsi cands
        double m0_jpsi = fPdg->GetParticle("J/psi")->Mass();   // Get nominal PDG mass of the J/psi
        jpsiMassSel    = new RhoMassParticleSelector("jpsi"   ,m0_jpsi,0.5);
        jpsiPreMassSel = new RhoMassParticleSelector("jpsipre",m0_jpsi,2.0);
        
        // ***
        // *** Prepare RhoTuple output  
        // ***
        TDirectory *dir = gDirectory;   
    fFile=new TFile(fOutName,"RECREATE");
        fFile->cd();
        
        // *** create some ntuples
        ntp1 = new RhoTuple("ntp1", "jpsi analysis");
        ntp2 = new RhoTuple("ntp2", "psi(2S) analysis");
        nmc  = new RhoTuple("nmc",  "mctruth info");

        // assign RhoTuples to outfile
        if (ntp1) ntp1->GetInternalTree()->SetDirectory(gDirectory);
        if (ntp2) ntp2->GetInternalTree()->SetDirectory(gDirectory);
        if (nmc)  nmc->GetInternalTree()->SetDirectory(gDirectory);

        // *** restore original gDirectory
        dir->cd();
        
        return kSUCCESS;
}

// -------------------------------------------------------------------------
        
void PndScrutAnaTask::SetParContainers() 
{
  // *** Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");
}

// -------------------------------------------------------------------------

// -----   Public method Exec   --------------------------------------------
void PndScrutAnaTask::Exec(Option_t*)
{
        // *** some variables
        int i=0,j=0;
        
        // *** necessary to read the next event
        fAnalysis->GetEvent();
        
        // *** print event counter
        if (!(++fEvtCount%100)) cout << "evt "<<fEvtCount<<endl;
        
        // *******
        // ******* PUT ANALYSIS CODE HERE                       
        // *******
        
        TString pidalg = "PidChargedProbability";
        
        // *** QA tool for simple dumping of analysis results in RhoRuple
        // *** WIKI: https://panda-wiki.gsi.de/foswiki/bin/view/Computing/PandaRootAnalysisJuly13#PndRhoTupleQA
        PndRhoTupleQA qa(fAnalysis,fIni.P());

        // *** RhoCandLists for the analysis
        RhoCandList muplus, muminus, piplus, piminus, jpsi, psi2s, all, mclist;
        
        
        // *** Select with PID info pidalg and ('All'); type and mass are set           
        fAnalysis->FillList(muplus,  "MuonAllPlus",  pidalg);
        fAnalysis->FillList(muminus, "MuonAllMinus", pidalg);
        fAnalysis->FillList(piplus,  "PionAllPlus",  pidalg);
        fAnalysis->FillList(piminus, "PionAllMinus", pidalg);

        
        // *** Setup event shape object
        fAnalysis->FillList(all,   "All", pidalg);
        PndEventShape evsh(all, fIni, 0.05, 0.1);       
        
        
        // *** store MC info in ntuple
        fAnalysis->FillList(mclist, "McTruth");
        
        nmc->Column("ev", (Int_t) fEvtCount);
        qa.qaMcList("",   mclist, nmc);
        nmc->DumpData();
        
        
        // *** combinatorics for J/psi -> mu+ mu-
        jpsi.Combine(muplus, muminus);
        jpsi.Select(jpsiPreMassSel);
        jpsi.SetType(443);
        int njmct = fAnalysis->McTruthMatch(jpsi); // match the whole list to count #matches (should be only 1)
                        
        // *** write ntuple for the jpsi reconstruction
        for (j=0;j<jpsi.GetLength();++j) 
        {
                // some general info about event (actually written for each candidate!)
                ntp1->Column("ev",              (Float_t) i);
                ntp1->Column("cand",    (Float_t) j);
                ntp1->Column("ncand",   (Float_t) jpsi.GetLength());
                ntp1->Column("nmct",    (Float_t) njmct);
        
                // store info about initial 4-vector
                qa.qaP4("beam", fIni, ntp1);
                        
                // dump information about composite candidate tree recursively (see analysis/AnalysisTools/PndRhoTupleQA)
                qa.qaComp("j", jpsi[j], ntp1);
                
                // dump info about event shapes
                qa.qaEventShapeShort("es",&evsh, ntp1);
                
                // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
                RhoCandidate *truth = jpsi[j]->GetMcTruth();
                TLorentzVector lv;
                if (truth) lv = truth->P4();
                qa.qaP4("trj", lv, ntp1);
                
                ntp1->DumpData();
        }
        
        
        // *** some rough mass selection on J/psi
        jpsi.Select(jpsiMassSel);
        
        // *** combinatorics for psi(2S) -> J/psi pi+ pi-
        psi2s.Combine(jpsi, piplus, piminus);
        psi2s.SetType(100443);
        int npsimct = fAnalysis->McTruthMatch(psi2s); // match the whole list to count #matches (should be only 1)

        // *** write ntuple for the psi(2S) reconstruction
        for (j=0;j<psi2s.GetLength();++j) 
        {
                // some general info about event (actually written for each candidate!)
                ntp2->Column("ev",              (Float_t) i);
                ntp2->Column("cand",    (Float_t) j);
                ntp2->Column("ncand",   (Float_t) psi2s.GetLength());
                ntp2->Column("nmct",    (Float_t) npsimct);
                
                RhoKinFitter kinfit(psi2s[j]);
                kinfit.Add4MomConstraint(fIni);
                kinfit.Fit();
                
                RhoCandidate *psifit=psi2s[j]->GetFit();
                // store info about initial 4-vector
                qa.qaP4("beam", fIni, ntp2);
                
                // dump information about composite candidate tree recursively (see analysis/AnalysisTools/PndRhoTupleQA)
                qa.qaComp("psi", psi2s[j], ntp2);
                qa.qaComp("fpsi",psifit, ntp2);
                ntp2->Column("fchi2", (Float_t) kinfit.GetChi2());

                // dump info about event shapes
                qa.qaEventShapeShort("es",&evsh, ntp2);
                
                // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
                RhoCandidate *truth = psi2s[j]->GetMcTruth();
                TLorentzVector lv;
                if (truth) lv = truth->P4();
                qa.qaP4("trpsi", lv, ntp2);
                
                ntp2->DumpData();
        }                       
        
}


void PndScrutAnaTask::Finish()
{
        
        // *******
        // ******* STORE YOUR HISTOS AND TUPLES
        // *******
        
        fFile->Write();
        fFile->Close(); 

}

ClassImp(PndScrutAnaTask)