// QA macro for testing event filter implementation with certain test conditions on output from dpm standalone // Martin Galuska, Katja Kleeberg sim_filter_qa(Int_t nEvents = 20, TString SimEngine ="TGeant3", Float_t mom = 15.) { //-----User Settings:----------------------------------------------- TString OutputFile ="sim_complete.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 =kFALSE; 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->SetWriteRunInfoFile(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 //------------------------------- FairPrimaryGenerator* primGen = new FairPrimaryGenerator(); //FilterMaxTries=99999, fSimulatedEvents=0, fFailedEvents=0, fVerbose=4 fRun->SetGenerator(primGen); gSystem->Load("libPGen"); PndDpmGenerator* dpmGen = new PndDpmGenerator("/home/katja/pandaroot/trunk/sim_katja/Background-micro.root"); primGen->AddGenerator(dpmGen); // now add the event filter FairEvtFilterOnCounts* fil= new FairEvtFilterOnCounts("PdgFilter"); //primGen->SetFilterMaxTries(100);//set the maximal number of attempts to find a suitable event to 150 Int_t evtNrsToAccept[] = {3,6,10,14,29,31,36,37,38,42,45,47,54,56,66,69,73,77,88,93}; Int_t nEntries = (sizeof(evtNrsToAccept)/sizeof(evtNrsToAccept[0])); fil->SetTestMode(evtNrsToAccept, nEntries); // for testing if event filter gives same results as before fil->AndMinMaxPdgCodes(3,5,211,-211);// min. 3 and max. 10 pions fil->AndMinAllParticles(4);// max. 4 particles in total //fil->AndMinMaxCharge(3,10,FairEvtFilterOnCounts::kPlus);// min. 3 and max. 10 pos. charged //fil->AndPzRange(0.,0.5);//transverse momentum within [0.,0.5] GeV/c primGen->SetVerbose();//highest commenting level of the FairPrimaryGenerator fil->SetVerbose();//highest commenting level of the FairEvtFilterOnCounts primGen->AndFilter(fil);//add filter to fFilterList if(UseBoxGenerator){ // Box Generator FairBoxGenerator* boxGenMuMinus = new FairBoxGenerator(13, 3); // 13 = muon; 1 = multipl. boxGenMuMinus->SetPRange(0.5,5.0); // GeV/c boxGenMuMinus->SetPhiRange(0., 360.); // Azimuth angle range [degree] boxGenMuMinus->SetThetaRange(0., 90.); // Polar angle in lab system range [degree] boxGenMuMinus->SetBoxXYZ(-1., -1., 1., 1., 10.); // cm primGen->AddGenerator(boxGenMuMinus); // FairBoxGenerator* boxGenMuPlus = new FairBoxGenerator(-13, 3); // 13 = muon; 3 = multipl. // boxGenMuPlus->SetPRange(0.5,5.0); // GeV/c // boxGenMuPlus->SetPhiRange(0., 360.); // Azimuth angle range [degree] // boxGenMuPlus->SetThetaRange(0., 90.); // Polar angle in lab system range [degree] // boxGenMuPlus->SetBoxXYZ(-1., -1., 1., 1., 10.); // cm // primGen->AddGenerator(boxGenMuPlus); } 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); } //---------------------Create and Set the Field(s)---------- PndMultiField *fField= new PndMultiField("FULL"); fRun->SetField(fField); // EMC Hit producer //------------------------------- PndEmcHitProducer* emcHitProd = new PndEmcHitProducer(); fRun->AddTask(emcHitProd); //------------------------- Initialize the RUN ----------------- fRun->Init(); //------------------------- Run the Simulation ----------------- fRun->Run(nEvents); //------------------------- 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); fil->ShowAcceptedEventNumbers(); // to generate the array of event numbers that were accepted by the filter with current settings fil->ShowEvtNrsToAccept(); if( kTRUE == fil->TestPassed() ){ std::cout << "\n\n\n QA test passed.\n\n\n"; }else{ std::cout << "\n\n\n QA test failed. \n\n\n"; } cout << " Test passed" << endl; cout << " All ok " << endl; //exit(0); }