StRoot: StJetMaker/macros/RunJetFinder2009sim.C Source File

00001 //
00002 // Pibero Djawotho <pibero@tamu.edu>
00003 // Texas A&M University
00004 // 24 Jan 2011
00005 //
00006 
00007 void RunJetFinder2009sim(int nevents = 1e6,
00008                          const char* mudstfile = "/star/data47/reco/pp200/pythia6_410/15_25gev/cdf_a/y2006c/gheisha_on/p07ic/rcf1307_01_2000evts.MuDst.root",
00009                          const char* geantfile = "/star/data47/reco/pp200/pythia6_410/15_25gev/cdf_a/y2006c/gheisha_on/p07ic/rcf1307_01_2000evts.geant.root",
00010                          const char* jetfile   = "rcf1307_01_2000evts.jets.root",
00011                          const char* skimfile  = "rcf1307_01_2000evts.skim.root",
00012                          bool useL2 = true)
00013 {
00014   cout << "Read MuDst file:\t" << mudstfile << endl;
00015   cout << "Read geant file:\t" << geantfile << endl;
00016   cout << "Write jet file:\t" << jetfile << endl;
00017   cout << "Write skim file:\t" << skimfile << endl;
00018 
00019   gROOT->Macro("loadMuDst.C");
00020   gROOT->Macro("LoadLogger.C");
00021 
00022   gSystem->Load("StDetectorDbMaker");
00023   gSystem->Load("StTpcDb");
00024   gSystem->Load("StDbUtilities");
00025   gSystem->Load("StMcEvent");
00026   gSystem->Load("StMcEventMaker");
00027   gSystem->Load("StDaqLib");
00028   gSystem->Load("StEmcRawMaker");
00029   gSystem->Load("StEmcADCtoEMaker");
00030   gSystem->Load("StPreEclMaker");
00031   gSystem->Load("StEmcSimulatorMaker");
00032   gSystem->Load("StDbBroker");
00033   gSystem->Load("St_db_Maker");
00034   gSystem->Load("StEEmcUtil");
00035   gSystem->Load("StEEmcDbMaker");
00036   gSystem->Load("StSpinDbMaker");
00037   gSystem->Load("StEmcTriggerMaker");
00038   gSystem->Load("StTriggerUtilities");
00039   gSystem->Load("StMCAsymMaker");
00040   gSystem->Load("StRandomSelector");
00041   gSystem->Load("libfastjet.so");
00042   gSystem->Load("libCDFConesPlugin.so");
00043   gSystem->Load("libEECambridgePlugin.so");
00044   gSystem->Load("libJadePlugin.so");
00045   gSystem->Load("libNestedDefsPlugin.so");
00046   gSystem->Load("libSISConePlugin.so");
00047   gSystem->Load("StJetFinder");
00048   gSystem->Load("StJetSkimEvent");
00049   gSystem->Load("StJets");
00050   gSystem->Load("StJetEvent");
00051   gSystem->Load("StJetMaker");
00052   gSystem->Load("StEEmcSimulatorMaker");
00053 
00054   // Create chain
00055   StChain* chain = new StChain;
00056 
00057   // I/O maker
00058   StIOMaker* ioMaker = new StIOMaker;
00059   ioMaker->SetFile(geantfile);
00060   ioMaker->SetIOMode("r");
00061   ioMaker->SetBranch("*",0,"0");             // Deactivate all branches
00062   ioMaker->SetBranch("geantBranch",0,"r");   // Activate geant Branch
00063 
00064   // StMcEvent maker
00065   StMcEventMaker* mcEventMaker = new StMcEventMaker;
00066   mcEventMaker->doPrintEventInfo = false;
00067   mcEventMaker->doPrintMemoryInfo = false;
00068 
00069   // MuDst reader
00070   StMuDstMaker* muDstMaker = new StMuDstMaker(0,0,"",mudstfile,"",100000,"MuDst");
00071 
00072   // MuDst DB
00073   StMuDbReader* muDstDb = StMuDbReader::instance();
00074 
00075   // star database
00076   St_db_Maker* starDb = new St_db_Maker("StarDb","MySQL:StarDb");
00077   starDb->SetDateTime(20090628,53220); // Run 10179006
00078   //starDb->SetDateTime(20060516,110349); // Run 7136022
00079 
00080   // Endcap database
00081   StEEmcDbMaker* eemcDb = new StEEmcDbMaker;
00082 
00083   // EEMC slow simulator
00084   StEEmcSlowMaker* eess = new StEEmcSlowMaker;
00085   eess->setSamplingFraction(0.0384);
00086   eess->setAddPed(true);
00087   eess->setSmearPed(true);
00088 
00089   // BEMC simulator
00090   StEmcSimulatorMaker* emcSim = new StEmcSimulatorMaker;
00091   emcSim->setCalibSpread(kBarrelEmcTowerId,0.15);
00092   emcSim->setCheckStatus(kBarrelEmcTowerId,false);
00093   emcSim->setMakeFullDetector(kBarrelEmcTowerId,true);
00094   emcSim->setDoZeroSuppression(kBarrelEmcTowerId,false);
00095 
00096   StPreEclMaker* preEcl = new StPreEclMaker; // need this to fill new StEvent information
00097 
00098   // Barrel ADC to energy maker
00099   StEmcADCtoEMaker* adc = new StEmcADCtoEMaker;
00100   adc->saveAllStEvent(true);
00101 
00102   // Trigger simulator
00103   // -- StL2_2009EmulatorMaker must run before StTriggerSimuMaker?
00104   StL2_2009EmulatorMaker* simL2Mk = 0;
00105   if (useL2) {
00106     simL2Mk = new StL2_2009EmulatorMaker;
00107     //simL2Mk->setSetupPath("/star/u/pibero/public/StarTrigSimuSetup/");
00108     simL2Mk->setSetupPath("/star/institutions/mit/corliss/L2setup/");
00109     simL2Mk->setOutPath("./");
00110   }
00111   StTriggerSimuMaker* simuTrig = new StTriggerSimuMaker;
00112   simuTrig->setMC(true); // Must be before individual detectors, to be passed
00113   //simuTrig->useBbc(); // No BBC in Run 9
00114   simuTrig->useBemc();
00115   simuTrig->useEemc(); // No EEMC in Run 6
00116   simuTrig->bemc->setConfig(StBemcTriggerSimu::kOffline);
00117   if (useL2) simuTrig->useL2(simL2Mk);
00118 
00119   // Get Pythia record
00120   StMCAsymMaker* asym = new StMCAsymMaker;
00121 
00122   // Skim event maker
00123   StJetSkimEventMaker* skimEventMaker = new StJetSkimEventMaker("StJetSkimEventMaker",muDstMaker,skimfile);
00124   skimEventMaker->addSimuTrigger(240530); // BHT3
00125   skimEventMaker->addSimuTrigger(240652); // L2JetHigh
00126   skimEventMaker->addSimuTrigger(240411); // JP1
00127 
00128 #if 0
00129   // Run 6 triggers
00130   skimEventMaker->addSimuTrigger(127611); //HTTP
00131   skimEventMaker->addSimuTrigger(137611);
00132   skimEventMaker->addSimuTrigger(5);
00133   skimEventMaker->addSimuTrigger(127821); //HTTP-fast
00134   skimEventMaker->addSimuTrigger(137821);
00135   skimEventMaker->addSimuTrigger(137822);
00136   skimEventMaker->addSimuTrigger(127212); //HT2
00137   skimEventMaker->addSimuTrigger(137213);
00138   skimEventMaker->addSimuTrigger(127501); //JP0
00139   skimEventMaker->addSimuTrigger(137501);
00140   skimEventMaker->addSimuTrigger(127622); //JP0-etot
00141   skimEventMaker->addSimuTrigger(137622);
00142   skimEventMaker->addSimuTrigger(127221); //JP1
00143   skimEventMaker->addSimuTrigger(137221);
00144   skimEventMaker->addSimuTrigger(137222);
00145 #endif
00146 
00147   // Jet maker
00148   StJetMaker2009* jetmaker = new StJetMaker2009;
00149   jetmaker->setJetFile(jetfile);
00150 
00151   //------------------------------------------------------------------------------------
00152 
00153   // Set analysis cuts for 12-point branch with 100% probability of accepting TPC tracks before jet finding
00154   StAnaPars* anapars12_100 = new StAnaPars;
00155   anapars12_100->useTpc  = true;
00156   anapars12_100->useBemc = true;
00157   anapars12_100->useEemc = true;
00158   anapars12_100->randomSelectorProb = 1.00;
00159 
00160   // The classes available for correcting tower energy for tracks are:
00161   // 1. StjTowerEnergyCorrectionForTracksMip
00162   // 2. StjTowerEnergyCorrectionForTracksFraction
00163   // 3. StjTowerEnergyCorrectionForTracksNull (default: no correction)
00164   anapars12_100->setTowerEnergyCorrection(new StjTowerEnergyCorrectionForTracksFraction(1.00));
00165 
00166   // TPC cuts
00167   anapars12_100->addTpcCut(new StjTrackCutFlag(0));
00168   anapars12_100->addTpcCut(new StjTrackCutNHits(12));
00169   anapars12_100->addTpcCut(new StjTrackCutPossibleHitRatio(0.51));
00170   anapars12_100->addTpcCut(new StjTrackCutDca(3));
00171   anapars12_100->addTpcCut(new StjTrackCutDcaPtDependent);
00172   anapars12_100->addTpcCut(new StjTrackCutChi2(0,4));
00173   anapars12_100->addTpcCut(new StjTrackCutPt(0.2,200));
00174   anapars12_100->addTpcCut(new StjTrackCutEta(-2.5,2.5));
00175   anapars12_100->addTpcCut(new StjTrackCutLastPoint(125));
00176 
00177   // BEMC cuts
00178   anapars12_100->addBemcCut(new StjTowerEnergyCutBemcStatus(1));
00179   anapars12_100->addBemcCut(new StjTowerEnergyCutAdc(4,3));     // ADC-ped>4 AND ADC-ped>3*RMS
00180   anapars12_100->addBemcCut(new StjTowerEnergyCutEt(0.2));
00181 
00182   // EEMC cuts
00183   anapars12_100->addEemcCut(new StjTowerEnergyCutBemcStatus(1));
00184   anapars12_100->addEemcCut(new StjTowerEnergyCutAdc(4,3));     // ADC-ped>4 AND ADC-ped>3*RMS
00185   anapars12_100->addEemcCut(new StjTowerEnergyCutEt(0.2));
00186 
00187   // Jet cuts
00188   anapars12_100->addJetCut(new StProtoJetCutPt(5,200));
00189   anapars12_100->addJetCut(new StProtoJetCutEta(-100,100));
00190 
00191   //------------------------------------------------------------------------------------
00192 
00193   // Set analysis cuts for 12-point branch with 93% probability of accepting TPC tracks before jet finding
00194   StAnaPars* anapars12_093 = new StAnaPars;
00195   anapars12_093->useTpc  = true;
00196   anapars12_093->useBemc = true;
00197   anapars12_093->useEemc = true;
00198   anapars12_093->randomSelectorProb = 0.93;
00199 
00200   // The classes available for correcting tower energy for tracks are:
00201   // 1. StjTowerEnergyCorrectionForTracksMip
00202   // 2. StjTowerEnergyCorrectionForTracksFraction
00203   // 3. StjTowerEnergyCorrectionForTracksNull (default: no correction)
00204   anapars12_093->setTowerEnergyCorrection(new StjTowerEnergyCorrectionForTracksFraction(1.00));
00205 
00206   // TPC cuts
00207   anapars12_093->addTpcCut(new StjTrackCutFlag(0));
00208   anapars12_093->addTpcCut(new StjTrackCutNHits(12));
00209   anapars12_093->addTpcCut(new StjTrackCutPossibleHitRatio(0.51));
00210   anapars12_093->addTpcCut(new StjTrackCutDca(3));
00211   anapars12_093->addTpcCut(new StjTrackCutDcaPtDependent);
00212   anapars12_093->addTpcCut(new StjTrackCutChi2(0,4));
00213   anapars12_093->addTpcCut(new StjTrackCutPt(0.2,200));
00214   anapars12_093->addTpcCut(new StjTrackCutEta(-2.5,2.5));
00215   anapars12_093->addTpcCut(new StjTrackCutLastPoint(125));
00216 
00217   // BEMC cuts
00218   anapars12_093->addBemcCut(new StjTowerEnergyCutBemcStatus(1));
00219   anapars12_093->addBemcCut(new StjTowerEnergyCutAdc(4,3));     // ADC-ped>4 AND ADC-ped>3*RMS
00220   anapars12_093->addBemcCut(new StjTowerEnergyCutEt(0.2));
00221 
00222   // EEMC cuts
00223   anapars12_093->addEemcCut(new StjTowerEnergyCutBemcStatus(1));
00224   anapars12_093->addEemcCut(new StjTowerEnergyCutAdc(4,3));     // ADC-ped>4 AND ADC-ped>3*RMS
00225   anapars12_093->addEemcCut(new StjTowerEnergyCutEt(0.2));
00226 
00227   // Jet cuts
00228   anapars12_093->addJetCut(new StProtoJetCutPt(5,200));
00229   anapars12_093->addJetCut(new StProtoJetCutEta(-100,100));
00230 
00231   //------------------------------------------------------------------------------------
00232 
00233   // Set analysis cuts for 5-point branch
00234   StAnaPars* anapars5 = new StAnaPars;
00235   anapars5->useTpc  = true;
00236   anapars5->useBemc = true;
00237   anapars5->useEemc = true;
00238   anapars5->randomSelectorProb = 0.93;
00239 
00240   // The classes available for correcting tower energy for tracks are:
00241   // 1. StjTowerEnergyCorrectionForTracksMip
00242   // 2. StjTowerEnergyCorrectionForTracksFraction
00243   // 3. StjTowerEnergyCorrectionForTracksNull (default: no correction)
00244   anapars5->setTowerEnergyCorrection(new StjTowerEnergyCorrectionForTracksFraction(1.00));
00245 
00246   // TPC cuts
00247   anapars5->addTpcCut(new StjTrackCutFlag(0));
00248   anapars5->addTpcCut(new StjTrackCutNHits(5));
00249   anapars5->addTpcCut(new StjTrackCutPossibleHitRatio(0.51));
00250   anapars5->addTpcCut(new StjTrackCutDca(3));
00251   anapars5->addTpcCut(new StjTrackCutDcaPtDependent);
00252   anapars5->addTpcCut(new StjTrackCutChi2(0,4));
00253   anapars5->addTpcCut(new StjTrackCutPt(0.2,200));
00254   anapars5->addTpcCut(new StjTrackCutEta(-2.5,2.5));
00255 
00256   // BEMC cuts
00257   anapars5->addBemcCut(new StjTowerEnergyCutBemcStatus(1));
00258   anapars5->addBemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
00259   anapars5->addBemcCut(new StjTowerEnergyCutEt(0.2));
00260 
00261   // EEMC cuts
00262   anapars5->addEemcCut(new StjTowerEnergyCutBemcStatus(1));
00263   anapars5->addEemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
00264   anapars5->addEemcCut(new StjTowerEnergyCutEt(0.2));
00265 
00266   // Jet cuts
00267   anapars5->addJetCut(new StProtoJetCutPt(5,200));
00268   anapars5->addJetCut(new StProtoJetCutEta(0.8,2.5));
00269 
00270   //------------------------------------------------------------------------------------
00271 
00272   // Set analysis cuts for EMC branch
00273   StAnaPars* anaparsEMC = new StAnaPars;
00274   anaparsEMC->useTpc  = true;
00275   anaparsEMC->useBemc = true;
00276   anaparsEMC->useEemc  = true;
00277 
00278   // TPC cuts
00279   anaparsEMC->addTpcCut(new StjTrackCutFlag(0));
00280   anaparsEMC->addTpcCut(new StjTrackCutNHits(1000000));
00281 
00282   // BEMC cuts
00283   anaparsEMC->addBemcCut(new StjTowerEnergyCutBemcStatus(1));
00284   anaparsEMC->addBemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
00285   anaparsEMC->addBemcCut(new StjTowerEnergyCutEt(0.2));
00286 
00287   // EEMC cuts
00288   anaparsEMC->addEemcCut(new StjTowerEnergyCutBemcStatus(1));
00289   anaparsEMC->addEemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
00290   anaparsEMC->addEemcCut(new StjTowerEnergyCutEt(0.2));
00291 
00292   // Jet cuts
00293   anaparsEMC->addJetCut(new StProtoJetCutPt(5,200));
00294   anaparsEMC->addJetCut(new StProtoJetCutEta(-100,100));
00295 
00296   //------------------------------------------------------------------------------------
00297 
00298   // Set analysis cuts for particle jets branch
00299   StAnaPars* anaparsParticle = new StAnaPars;
00300   anaparsParticle->useMonteCarlo = true;
00301 
00302   // MC cuts
00303   anaparsParticle->addMcCut(new StjMCParticleCutStatus(1)); // final state particles
00304 
00305   // Jet cuts
00306   anaparsParticle->addJetCut(new StProtoJetCutPt(3,200));
00307   anaparsParticle->addJetCut(new StProtoJetCutEta(-100,100));
00308     
00309   //------------------------------------------------------------------------------------
00310 
00311   // Set analysis cuts for parton jets branch
00312   StAnaPars* anaparsParton = new StAnaPars;
00313   anaparsParton->useMonteCarlo = true;
00314 
00315   // MC cuts
00316   anaparsParton->addMcCut(new StjMCParticleCutParton);
00317 
00318   // Jet cuts
00319   anaparsParton->addJetCut(new StProtoJetCutPt(3,200));
00320   anaparsParton->addJetCut(new StProtoJetCutEta(-100,100));
00321     
00322   // Set cone jet finder parameters
00323   StConePars* conepars = new StConePars;
00324   conepars->setGridSpacing(105,-3.0,3.0,120,-TMath::Pi(),TMath::Pi());
00325   conepars->setConeRadius(0.7);
00326   conepars->setSeedEtMin(0.5);
00327   conepars->setAssocEtMin(0.1);
00328   conepars->setSplitFraction(0.5);
00329   conepars->setPerformMinimization(true);
00330   conepars->setAddMidpoints(true);
00331   conepars->setRequireStableMidpoints(true);
00332   conepars->setDoSplitMerge(true);
00333   conepars->setDebug(false);
00334 
00335   //------------------------------------------------------------------------------------
00336 
00337   jetmaker->addBranch("ConeJets12_100",anapars12_100,conepars);
00338   jetmaker->addBranch("ConeJets12_093",anapars12_093,conepars);
00339   jetmaker->addBranch("ConeJets5",anapars5,conepars);
00340   jetmaker->addBranch("ConeJetsEMC",anaparsEMC,conepars);
00341   jetmaker->addBranch("ParticleConeJets",anaparsParticle,conepars);
00342   jetmaker->addBranch("PartonConeJets",anaparsParton,conepars);
00343 
00344   //------------------------------------------------------------------------------------
00345 
00346   chain->Init();
00347   chain->EventLoop(nevents);
00348 }