sim_filter_inv_mass.C

Go to the documentation of this file.

// Example for usage of event filter on invariant masses in combination with dpm direct
// Martin Galuska


sim_filter_inv_mass(Int_t nEvents = 5, TString  SimEngine ="TGeant3", Float_t mom = 6.231552)
{
        //-----User Settings:-----------------------------------------------
        TString  OutputFile     ="sim_filter_inv_mass.root";
        TString  ParOutputfile  ="simparams.root";
        TString  MediaFile      ="media_pnd.geo";
        gDebug                  = 0;
        TString digiFile        = "all.par"; //The emc run the hit producer directly
        // choose your event generator
        Bool_t UseEvtGen              =kFALSE;
        Bool_t UseEvtGenDirect      =kFALSE;
        Bool_t UseDpm         =kTRUE;
        Bool_t UseBoxGenerator      =kFALSE;

        Double_t BeamMomentum = 0.; // beam momentum ONLY for the scaling of the dipole field.
        if (UseBoxGenerator)
        {
                BeamMomentum   =15.0; // ** change HERE if you run Box generator
        }
        else
        {
                BeamMomentum = mom;  // for DPM/EvtGen BeamMomentum is always = mom
        }

        //------------------------------------------------------------------
        TStopwatch timer;
        timer.Start();
        gRandom->SetSeed();

        // Create the Simulation run manager--------------------------------
        FairRunSim *fRun = new FairRunSim();
        fRun->SetName(SimEngine.Data() );
        fRun->SetOutputFile(OutputFile.Data());
        fRun->SetGenerateRunInfo(kFALSE);
        fRun->SetBeamMom(BeamMomentum);
        fRun->SetMaterials(MediaFile.Data());
        FairRuntimeDb *rtdb=fRun->GetRuntimeDb();

        // Set the parameters
        //-------------------------------
        TString allDigiFile = gSystem->Getenv("VMCWORKDIR");
        allDigiFile += "/macro/params/";
        allDigiFile += digiFile;


        //-------Set the parameter output --------------------
        FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
        parIo1->open(allDigiFile.Data(),"in");
        rtdb->setFirstInput(parIo1);

        //---------------------Set Parameter output      ----------
        Bool_t kParameterMerged=kTRUE;
        FairParRootFileIo* output=new FairParRootFileIo(kParameterMerged);
        output->open(ParOutputfile.Data());
        rtdb->setOutput(output);

        // Create and add detectors

        //-------------------------  CAVE      -----------------

        FairModule *Cave= new PndCave("CAVE");
        Cave->SetGeometryFileName("pndcave.geo");
        fRun->AddModule(Cave);
        //-------------------------  Magnet   -----------------
        //      FairModule *Magnet= new PndMagnet("MAGNET");
        //      //Magnet->SetGeometryFileName("FullSolenoid_V842.root");
        //      Magnet->SetGeometryFileName("FullSuperconductingSolenoid_v831.root");
        //      fRun->AddModule(Magnet);
        //      FairModule *Dipole= new PndMagnet("MAGNET");
        //      Dipole->SetGeometryFileName("dipole.geo");
        //      fRun->AddModule(Dipole);
        //      //-------------------------  Pipe     -----------------
        //      FairModule *Pipe= new PndPipe("PIPE");
        //      Pipe->SetGeometryFileName("beampipe_201309.root");
        //      fRun->AddModule(Pipe);
        //      //-------------------------  STT       -----------------
        //      FairDetector *Stt= new PndStt("STT", kTRUE);
        //      Stt->SetGeometryFileName("straws_skewed_blocks_35cm_pipe.geo");
        //      fRun->AddModule(Stt);
        //      //-------------------------  MVD       -----------------
        //      FairDetector *Mvd = new PndMvdDetector("MVD", kTRUE);
        //      Mvd->SetGeometryFileName("Mvd-2.1_FullVersion.root");
        //      fRun->AddModule(Mvd);
        //      //-------------------------  GEM       -----------------
        //      FairDetector *Gem = new PndGemDetector("GEM", kTRUE);
        //      Gem->SetGeometryFileName("gem_3Stations.root");
        //      fRun->AddModule(Gem);
        //      //-------------------------  EMC       -----------------
        //      PndEmc *Emc = new PndEmc("EMC",kTRUE);
        //      Emc->SetGeometryVersion(1);
        //      Emc->SetStorageOfData(kFALSE);
        //      fRun->AddModule(Emc);
        //      //-------------------------  SCITIL    -----------------
        //      FairDetector *SciT = new PndSciT("SCIT",kTRUE);
        //      SciT->SetGeometryFileName("barrel-SciTil_07022013.root");
        //      fRun->AddModule(SciT);
        //      //-------------------------  DRC       -----------------
        //      PndDrc *Drc = new PndDrc("DIRC", kTRUE);
        //      Drc->SetGeometryFileName("dirc_l0_p0_updated.root");
        //      Drc->SetRunCherenkov(kFALSE);
        //      fRun->AddModule(Drc);
        //      //-------------------------  DISC      -----------------
        //      PndDsk* Dsk = new PndDsk("DSK", kTRUE);
        //      Dsk->SetStoreCerenkovs(kFALSE);
        //      Dsk->SetStoreTrackPoints(kFALSE);
        //      fRun->AddModule(Dsk);
        //      //-------------------------  MDT       -----------------
        //      PndMdt *Muo = new PndMdt("MDT",kTRUE);
        //      Muo->SetBarrel("fast");
        //      Muo->SetEndcap("fast");
        //      Muo->SetMuonFilter("fast");
        //      Muo->SetForward("fast");
        //      Muo->SetMdtMagnet(kTRUE);
        //      Muo->SetMdtMFIron(kTRUE);
        //      fRun->AddModule(Muo);
        //      //-------------------------  FTS       -----------------
        //      FairDetector *Fts= new PndFts("FTS", kTRUE);
        //      Fts->SetGeometryFileName("fts.geo");
        //      fRun->AddModule(Fts);
        //      //-------------------------  FTOF      -----------------
        //      FairDetector *FTof = new PndFtof("FTOF",kTRUE);
        //      FTof->SetGeometryFileName("ftofwall.root");
        //      fRun->AddModule(FTof);
        //      //-------------------------  RICH       ----------------
        //      FairDetector *Rich= new PndRich("RICH",kFALSE);
        //      Rich->SetGeometryFileName("rich_v2_shift.geo");
        //      fRun->AddModule(Rich);

        // Create and Set Event Generator
        //-------------------------------
        FairFilteredPrimaryGenerator* primGen = new FairFilteredPrimaryGenerator();
        fRun->SetGenerator(primGen);

        if(UseBoxGenerator){    // Box Generator
                FairBoxGenerator* boxGen = new FairBoxGenerator(22, 5); // 13 = muon; 1 = multipl.
                boxGen->SetPRange(mom,mom); // GeV/c
                boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
                boxGen->SetThetaRange(0., 90.); // Polar angle in lab system range [degree]
                boxGen->SetXYZ(0., 0., 0.); // cm
                primGen->AddGenerator(boxGen);
        }
        if(UseDpm){
                PndDpmDirect *Dpm= new PndDpmDirect(mom,1);
                primGen->AddGenerator(Dpm);
        }
        if(UseEvtGen){
                TString  EvtInput =gSystem->Getenv("VMCWORKDIR");
                EvtInput+="/input/psi2s_jpsi2pi_1k.evt";
                FairEvtGenGenerator* evtGen = new FairEvtGenGenerator(EvtInput.Data());
                primGen->AddGenerator(evtGen);
        }
        if(UseEvtGenDirect){
                TString  EvtInput =gSystem->Getenv("VMCWORKDIR");
                EvtInput+="/macro/run/psi2s_Jpsi2pi_Jpsi_mumu.dec";
                PndEvtGenDirect *EvtGen = new PndEvtGenDirect("pbarpSystem", EvtInput.Data(), mom);
                EvtGen->SetStoreTree(kTRUE);
                primGen->AddGenerator(EvtGen);
        }

        // now add the event filters

        //              primGen->SetVerbose();//highest commenting level of the FairPrimaryGenerator
        //              primGen->SetFilterMaxTries(10000);

        // filter on e+ e- invariant mass
        PndEvtFilterOnInvMassCounts* eeInv= new PndEvtFilterOnInvMassCounts("eeInvMFilter");
        eeInv->SetVerbose();//highest commenting level of the FairEvtFilterOnCounts
        eeInv->SetPdgCodesToCombine( 11, -11);
        eeInv->SetMinMaxInvMass( 2.0, 4.0 );
        eeInv->SetMinMaxCounts(1,10000);
        primGen->AndFilter(eeInv);//add filter to fFilterList

        //---------------------Create and Set the Field(s)----------
        PndMultiField *fField= new PndMultiField("AUTO");
        fRun->SetField(fField);

        // EMC Hit producer
        //-------------------------------
        PndEmcHitProducer* emcHitProd = new PndEmcHitProducer();
        fRun->AddTask(emcHitProd);

        //-------------------------  Initialize the RUN  -----------------
        fRun->Init();
        //-------------------------  Run the Simulation  -----------------
        fRun->Run(nEvents);
        primGen->WriteEvtFilterStatsToRootFile();
        //-------------------------  Save the parameters -----------------
        rtdb->saveOutput();
        //------------------------Print some info and exit----------------
        timer.Stop();
        Double_t rtime = timer.RealTime();
        Double_t ctime = timer.CpuTime();
        printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

        cout << " Test passed" << endl;
        cout << " All ok " << endl;

        //exit(0);

}