isoCluster/macro/Fun4All_testEtIso_Pythia.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 #include <fun4all/SubsysReco.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllInputManager.h>
0006 #include <fun4all/Fun4AllDstInputManager.h>
0007
0008 #include <fun4all/Fun4AllDstOutputManager.h>
0009 #include <fun4all/Fun4AllOutputManager.h>
0010 #include <fun4all/Fun4AllServer.h>
0011
0012 #include <phool/PHRandomSeed.h>
0013 #include <phool/recoConsts.h>
0014
0015 #include <g4jets/FastJetAlgo.h>
0016 #include <g4jets/JetReco.h>
0017 #include <g4jets/TowerJetInput.h>
0018 #include <g4jets/TruthJetInput.h>
0019
0020 #include <jetbackground/CopyAndSubtractJets.h>
0021 #include <jetbackground/DetermineTowerBackground.h>
0022 #include <jetbackground/FastJetAlgoSub.h>
0023 #include <jetbackground/RetowerCEMC.h>
0024 #include <jetbackground/SubtractTowers.h>
0025 #include <jetbackground/SubtractTowersCS.h>
0026
0027 #include <ClusterIso.h>
0028
0029 #include <isocluster/isoCluster.h>
0030
0031 R__LOAD_LIBRARY(libfun4all.so)
0032 R__LOAD_LIBRARY(libisoCluster.so)
0033 R__LOAD_LIBRARY(libg4jets.so)
0034 R__LOAD_LIBRARY(libjetbackground.so)
0035 R__LOAD_LIBRARY(libclusteriso.so)
0036
0037 namespace Enable
0038 {
0039   bool HIJETS = false;
0040   int HIJETS_VERBOSITY = 1;
0041 }  // namespace Enable
0042
0043 namespace G4HIJETS
0044 {
0045   bool do_flow = true;
0046   bool do_CS = false;
0047 }  // namespace G4HIJETS
0048 #endif
0049
0050
0051 void Fun4All_testEtIso_Pythia(
0052                      int nEvents = 1,
0053                      const char *filelist1 = "dst_calo_cluster_run40_PhotonJet.list",
0054                      const char *filelist2 = "dst_truth_run40_PhotonJet.list",
0055              const string outname = "test.root")
0056 {
0057
0058
0059
0060   Fun4AllServer *se = Fun4AllServer::instance();
0061   int verbosity = 0;
0062
0063   se->Verbosity(verbosity);
0064   recoConsts *rc = recoConsts::instance();
0065   //rc->set_IntFlag("RUNNUMBER",62);
0066
0067
0068   RetowerCEMC *rcemc = new RetowerCEMC();
0069   //rcemc->SetEnergyDistribution(1);
0070   se->registerSubsystem(rcemc);
0071
0072   //seeds
0073   JetReco *towerjetreco = new JetReco();
0074   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWER_RETOWER));
0075   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWER));
0076   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWER));
0077   towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2), "AntiKt_Tower_HIRecoSeedsRaw_r02");
0078   towerjetreco->set_algo_node("ANTIKT");
0079   towerjetreco->set_input_node("TOWER");
0080
0081   se->registerSubsystem(towerjetreco);
0082
0083
0084
0085   //first UE subtraction pass
0086   DetermineTowerBackground *dtb = new DetermineTowerBackground();
0087   dtb->SetBackgroundOutputName("TowerBackground_Sub1");
0088   dtb->SetFlow(G4HIJETS::do_flow);
0089   dtb->SetSeedType(0);
0090   dtb->SetSeedJetD(3);
0091
0092   se->registerSubsystem(dtb);
0093
0094   CopyAndSubtractJets *casj = new CopyAndSubtractJets();
0095   casj->SetFlowModulation(G4HIJETS::do_flow);
0096   casj->Verbosity(verbosity);
0097   se->registerSubsystem(casj);
0098
0099
0100
0101   DetermineTowerBackground *dtb2 = new DetermineTowerBackground();
0102   dtb2->SetBackgroundOutputName("TowerBackground_Sub2");
0103   dtb2->SetFlow(G4HIJETS::do_flow);
0104   dtb2->SetSeedType(1);
0105   dtb2->SetSeedJetPt(7);
0106   dtb2->Verbosity(verbosity);
0107   se->registerSubsystem(dtb2);
0108
0109   //subtract UE from towers
0110   SubtractTowers *st = new SubtractTowers();
0111   st->SetFlowModulation(G4HIJETS::do_flow);
0112   st->Verbosity(verbosity);
0113   se->registerSubsystem(st);
0114
0115    //I think this goes here?
0116   ClusterIso *makeClusterEt = new ClusterIso("isoCluster",0,3,1,0);
0117   makeClusterEt -> Verbosity(0);
0118   se->registerSubsystem(makeClusterEt);
0119   ///???
0120   //substract UE from jets
0121   /*CopyAndSubtractJets *casj = new CopyAndSubtractJets();
0122   casj->SetFlowModulation(G4HIJETS::do_flow);
0123
0124   se->registerSubsystem(casj);
0125
0126   //second pass of UE subtraction
0127   DetermineTowerBackground *dtb2 = new DetermineTowerBackground();
0128   dtb2->SetBackgroundOutputName("TowerBackground_Sub2");
0129   dtb2->SetFlow(G4HIJETS::do_flow);
0130   dtb2->SetSeedType(1);
0131   dtb2->SetSeedJetPt(7);
0132
0133   se->registerSubsystem(dtb2);*/
0134
0135   //now we build our regular jets. Not used for eT analysis right now
0136   //but could add in function easily.
0137   /*towerjetreco = new JetReco();
0138   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWER_SUB1));
0139   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWER_SUB1));
0140   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWER_SUB1));
0141   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2, 1), "AntiKt_Tower_r02_Sub1");
0142   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.3, 1), "AntiKt_Tower_r03_Sub1");
0143   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4, 1), "AntiKt_Tower_r04_Sub1");
0144   towerjetreco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.5, 1), "AntiKt_Tower_r05_Sub1");
0145   towerjetreco->set_algo_node("ANTIKT");
0146   towerjetreco->set_input_node("TOWER");
0147   towerjetreco->Verbosity(verbosity);
0148   se->registerSubsystem(towerjetreco);*/
0149
0150
0151
0152   Fun4AllInputManager *in = new Fun4AllDstInputManager("DSTcalo");
0153   in->AddListFile(filelist1,1);
0154   se->registerInputManager(in);
0155
0156   Fun4AllInputManager *in2 = new Fun4AllDstInputManager("DSTTruth");
0157   in2 -> AddListFile(filelist2,1);
0158   se -> registerInputManager(in2);
0159
0160   isoCluster *eval = new isoCluster(outname.c_str());
0161   eval -> setGenEvent(1);
0162   se -> registerSubsystem(eval);
0163
0164   se -> run(nEvents);
0165   se -> End();
0166
0167 }