tutorials/rho/quickana.C

Go to the documentation of this file.

// ********************************************
// Quick analysis using PndSimpleCombinerTask/PndSimpleCombiner
// ********************************************

// The parameters are
// -------------------
// USAGE
// quickana.C( <input>, <mom>, <decay>, [nevt], [parms], [fastsim], [runST], [runnum], [mode] )
//    <input>   : input file name with PndPidCandidates
//    <mom>     : pbar momentum; negative values are interpreted as -E_cm
//    <decay>   : the decay pattern to be reconstructed, e.g. 'phi -> K+ K-; D_s+ -> phi pi-'
//    [nevt]    : number of events; default: 0 = all
//    [parms]   : parameters for the analysis, e.g. 'mwin=0.4:mwin(phi)=0.1:emin=0.1:pmin=0.1:qamc'
//    [fastsim] : set true, if running fast sim (sets the PID algos properly); default: false'
//    [runST]   : if 'true' runs Software Trigger (default: false)
//    [runnum]  : integer run number (default: 0)
//    [mode]    : arbitrary mode number (default: 0)
// -------------------

void quickana(TString Fname="", double Mom=0, TString anadecay="", int nevts=0, TString anaparms="", bool fastsim=false, bool runST=false, int run=0, int runmode=0)
{
        if (Fname=="" || Mom==0) 
        {
                cout << "USAGE:\n";
                cout << "quickana.C( <input>, <mom>, <decay>, [nevt], [parms], [fastsim], [runST], [runnum], [mode] )\n\n";
                cout << "   <input>   : input file name with PndPidCandidates\n";
                cout << "   <mom>     : pbar momentum; negative values are interpreted as -E_cm\n";
                cout << "   <decay>   : the decay pattern to be reconstructed, e.g. 'phi -> K+ K-; D_s+ -> phi pi-'\n";
                cout << "   [nevt]    : number of events; default: 0 = all\n";
                cout << "   [parms]   : parameters for the analysis, e.g. 'mwin=0.4:mwin(phi)=0.1:emin=0.1:pmin=0.1:qamc'; 'qapart' runs particle QA task\n";
                cout << "   [fastsim] : set true, if running fast sim (sets the PID algos properly); default: false'\n";
                cout << "   [runST]   : if 'true' runs Software Trigger (default: false)\n";
                cout << "   [runnum]  : integer run number (default: 0)\n";
                cout << "   [mode]    : arbitrary mode number (default: 0)\n\n";
                return;
        }
        
        // do some reconstruction ?
        bool doreco  = (anadecay != "" || anaparms.Contains("nevt"));

        // do particle QA?
        bool partQA  = (anaparms.Contains("qapart"));
        bool mc      = !(anaparms.Contains("!mc")) && !(anaparms.Contains("qamc"));
        bool neut    = !(anaparms.Contains("!neut"));
        bool chrg    = !(anaparms.Contains("!chrg"));
        
        // if Mom<0, interprete as -E_cm
        double mp = 0.938272;
        
        // if mom<0, it's -E_cm -> compute mom
        if (Mom<0)
        {
                double X = (Mom*Mom-2*mp*mp)/(2*mp);
                Mom = sqrt(X*X-mp*mp);
        }
        
        // PID algorithm for the PndSimpleCombinerTask (for Eventshape variables)
        TString pidalgo = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts;PidAlgoRich;PidAlgoSciT";
        if (fastsim) pidalgo = "PidChargedProbability";
        
        // allow shortcuts
        anadecay.ReplaceAll("pbp","pbarpSystem");
        anadecay.ReplaceAll("pbp0","pbarpSystem0");
        
        // Prevent generator from throwing a lot of warnings
        //TLorentzVector fIni(0,0,Mom,0.938272+sqrt(Mom*Mom+0.938272*0.938272));
        TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",3,kFALSE,0.1,0, "",88888);
        TDatabasePDG::Instance()->AddParticle("pbarpSystem0","pbarpSystem0",3,kFALSE,0.1,0, "",88880);
        
        // *** set this to your output path
        TString OutFile = Fname;//(Fname.Last('/')+1,Fname.Length()); // cut away input path
        OutFile.ReplaceAll(".root","_ana.root");
        
        // *** the output file for FairRunAna
        TString InFile  = Fname;
        if (!InFile.EndsWith(".root")) InFile+="_fast.root";
                                        
        // *** initialization
        FairLogger::GetLogger()->SetLogToFile(kFALSE);

        FairRunAna* fRun = new FairRunAna();
        fRun->SetGenerateRunInfo(kFALSE);
        fRun->SetSource(new FairFileSource(InFile));
        fRun->SetOutputFile(OutFile);

        // *** take constant field; needed for PocaVtx
        RhoCalculationTools::ForceConstantBz(20.0);
        
        // ***********************
        // *** SoftTriggerTask ***
        // ***********************
        
        if (runST)
        {       
                // this file contains the trigger line definitions
                TString      triggercfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/triggerlines.cfg";       // fullsim trigger definitions              
                if (fastsim) triggercfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/triggerlines_fsim.cfg";  // fastsim trigger definitions      
                
                PndSoftTriggerTask *stTask = new PndSoftTriggerTask(Mom, 0, run, triggercfg);
                
                if (fastsim) stTask->SetFastSimDefaults();
                else         stTask->SetFullSimDefaults();
                                
                fRun->AddTask(stTask);
        }
        
        // --------------------------------
        // *** Analysis Task ***
        // --------------------------------

        // *****************************
        // *** PndSimpleCombinerTask ***
        // *****************************
        
        if (doreco)
        {
                if (fastsim) anaparms+=":algo="+pidalgo;
                PndSimpleCombinerTask *scTask = new PndSimpleCombinerTask(anadecay, anaparms, Mom, run, runmode);
                scTask->SetPidAlgo(pidalgo);
                fRun->AddTask(scTask);
        }
        
        // *****************************
        // *** PndParticleQATask ***
        // *****************************
        
        if (partQA)
        {
                PndParticleQATask *partQaTask = new PndParticleQATask(fastsim,chrg,neut,mc); // particle QA task
                fRun->AddTask(partQaTask);
        }
                

        // *** and run analysis
        fRun->Init(); 
        fRun->Run(0,nevts);     
}