ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.1
Committed: Wed Oct 15 06:05:00 2008 UTC (16 years, 6 months ago) by loizides
Content type: text/plain
Branch: MAIN
Log Message:
Added MitPhysics.

File Contents

# User Rev Content
1 loizides 1.1 // $Id: GeneratorMod.cc,v 1.1 2008/10/14 06:13:52 loizides Exp $
2    
3     #include "MitPhysics/Mods/interface/GeneratorMod.h"
4     #include "MitAna/DataTree/interface/Names.h"
5     #include "MitAna/DataCont/interface/ObjArray.h"
6     #include "MitCommon/MathTools/interface/MathUtils.h"
7     #include <TH1D.h>
8     #include <TH2D.h>
9    
10    
11     using namespace mithep;
12    
13     ClassImp(mithep::GeneratorMod)
14    
15     //--------------------------------------------------------------------------------------------------
16     GeneratorMod::GeneratorMod(const char *name, const char *title) :
17     BaseMod(name,title),
18     fPrintDebug(false),
19     fFillHist(false),
20     fIsMC(true),
21     fMCPartName(Names::gkMCPartBrn),
22     fMCLeptonsName(Names::gkMCLeptonsName),
23     fMCTausName(Names::gkMCTausName),
24     fMCNeutrinosName(Names::gkMCNeutrinosName),
25     fMCQuarksName(Names::gkMCQuarksName),
26     fMCqqHsName(Names::gkMCqqHsName),
27     fMCBosonsName(Names::gkMCBosonsName),
28     fParticles(0),
29     fNEventsProcessed(0)
30     {
31     // Constructor.
32     }
33    
34     //--------------------------------------------------------------------------------------------------
35     void GeneratorMod::Begin()
36     {
37     // Run startup code on the client machine. For this module, we dont do
38     // anything here.
39     }
40    
41     //--------------------------------------------------------------------------------------------------
42     void GeneratorMod::Process()
43     {
44     // Process entries of the tree
45     fNEventsProcessed++;
46    
47     if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
48     time_t systime;
49     systime = time(NULL);
50     cerr << endl << "GeneratorMod : Process Event " << fNEventsProcessed << " Time: "
51     << ctime(&systime) << endl;
52     }
53    
54     // These arrays will be filled in the loop of particles
55     ObjArray<MCParticle> *GenLeptons = new ObjArray<MCParticle>;
56     ObjArray<MCParticle> *GenTaus = new ObjArray<MCParticle>; GenTaus->SetOwner(true);
57     ObjArray<MCParticle> *GenNeutrinos = new ObjArray<MCParticle>;
58     ObjArray<MCParticle> *GenQuarks = new ObjArray<MCParticle>;
59     ObjArray<MCParticle> *GenqqHs = new ObjArray<MCParticle>;
60     ObjArray<MCParticle> *GenBosons = new ObjArray<MCParticle>;
61    
62     if(fIsMC == true){
63     // Get Generator Level information branch
64     LoadBranch(fMCPartName);
65     bool isqqH = false;
66     for (UInt_t i=0; i<fParticles->GetEntries(); ++i) {
67     MCParticle* p = fParticles->At(i);
68    
69     if(!p->IsGenerated()) continue;
70    
71     // muons/electrons from W/Z decays
72     if((p->AbsPdgId() == 11 || p->AbsPdgId() == 13) && p->Status() == 1){
73     bool isGoodLepton = false;
74     MCParticle* pm = p;
75     while (pm->HasMother() && isGoodLepton == false){
76     if (pm->Mother()->AbsPdgId() == 23 || pm->Mother()->AbsPdgId() == 24){
77     GenLeptons->Add(p);
78     isGoodLepton = true;
79     }
80     else if(pm->Mother()->AbsPdgId() == 111 || pm->Mother()->AbsPdgId() == 221){
81     // This is fake, but it is a trick to get rid of these cases
82     isGoodLepton = true;
83     }
84     else {
85     pm = (MCParticle *)pm->Mother();
86     }
87     }
88     }
89    
90     // taus
91     else if(p->AbsPdgId() == 16 && p->Status() == 1){
92     if(p->DistinctMother()){
93     MCParticle* pm = (mithep::MCParticle*)p->DistinctMother();
94     if(pm->AbsPdgId() == 15){
95     MCParticle* pm_f = new MCParticle(*pm);
96     pm_f->SetMom(pm->Px()-p->Px(), pm->Py()-p->Py(),
97     pm->Pz()-p->Pz(), pm->E()-p->E());
98     GenTaus->Add(pm_f);
99     }
100     }
101     }
102    
103     // neutrinos
104     else if(p->Status() == 1 &&
105     (p->AbsPdgId() == 12 || p->AbsPdgId() == 14 || p->AbsPdgId() == 16)){
106     GenNeutrinos->Add(p);
107     }
108    
109     // quarks from W/Z decays or top particles
110     else if(p->AbsPdgId() >=1 && p->AbsPdgId() <=6 && p->HasMother()){
111     if(p->Mother()->AbsPdgId() == 23 || p->Mother()->AbsPdgId() == 24 ||
112     p->AbsPdgId() == 6 || p->Mother()->AbsPdgId() == 6){
113     GenQuarks->Add(p);
114     }
115     }
116    
117     // qqH, information about the forward jets
118     else if(isqqH == false && p->AbsPdgId() == 25){
119     isqqH = true;
120     MCParticle* pq1 = fParticles->At(i-1);
121     MCParticle* pq2 = fParticles->At(i-2);
122    
123     if(pq1->HasMother() && pq2->HasMother() &&
124     pq1->Mother()->PdgId() == p->Mother()->PdgId() &&
125     pq2->Mother()->PdgId() == p->Mother()->PdgId() &&
126     pq1->AbsPdgId() < 7 && pq2->AbsPdgId() < 7 &&
127     pq1->AbsPdgId() > 0 && pq2->AbsPdgId() > 0){
128     GenqqHs->Add(pq1);
129     GenqqHs->Add(pq2);
130     }
131     }
132    
133     // information about bosons: W, Z, h, Z', W', H0, A0, H+
134     else if(p->Status() == 2 &&
135     (p->AbsPdgId() == 23 || p->AbsPdgId() == 24 || p->AbsPdgId() == 25 ||
136     p->AbsPdgId() == 32 || p->AbsPdgId() == 34 ||
137     p->AbsPdgId() == 35 || p->AbsPdgId() == 36 || p->AbsPdgId() == 37)){
138     GenBosons->Add(p);
139     }
140     }
141     } // IsMC?
142    
143     //Save Objects for Other Modules to use
144     AddObjThisEvt(GenLeptons, fMCLeptonsName.Data());
145     AddObjThisEvt(GenTaus, fMCTausName.Data());
146     AddObjThisEvt(GenNeutrinos,fMCNeutrinosName.Data());
147     AddObjThisEvt(GenQuarks, fMCQuarksName.Data());
148     AddObjThisEvt(GenqqHs, fMCqqHsName.Data());
149     AddObjThisEvt(GenBosons, fMCBosonsName.Data());
150    
151     // Fill histograms
152     if(fFillHist == true){
153     // Leptons
154     hDGenLeptons[0]->Fill(GenLeptons->GetEntries());
155     int idxMaxLep[2] = {-1, -1};
156     double ptMaxLep[2] = {-1.0, -1.0};
157     for(UInt_t i=0; i<GenLeptons->GetEntries(); i++){
158     hDGenLeptons[1]->Fill(GenLeptons->At(i)->Pt());
159     hDGenLeptons[2]->Fill(fabs(GenLeptons->At(i)->Eta()));
160     hDGenLeptons[3]->Fill(GenLeptons->At(i)->Phi() * 180. / TMath::Pi());
161     for(UInt_t j=i+1; j<GenLeptons->GetEntries(); j++){
162     CompositeParticle *dilepton = new CompositeParticle();
163     dilepton->AddDaughter(GenLeptons->At(i));
164     dilepton->AddDaughter(GenLeptons->At(j));
165     hDGenLeptons[4]->Fill(dilepton->Mass());
166     delete dilepton;
167     }
168     // Selecting the two highest Pt leptons
169     if (GenLeptons->At(i)->Pt() > ptMaxLep[0]){
170     ptMaxLep[1] = ptMaxLep[0];
171     idxMaxLep[1] = idxMaxLep[0];
172     ptMaxLep[0] = GenLeptons->At(i)->Pt();
173     idxMaxLep[0] = i;
174     }
175     else if(GenLeptons->At(i)->Pt() > ptMaxLep[1]){
176     ptMaxLep[1] = GenLeptons->At(i)->Pt();
177     idxMaxLep[1] = i;
178     }
179     }
180     // Looking at events with at least two leptons
181     if(ptMaxLep[0] > 0 && ptMaxLep[1] > 0){
182     hDGenLeptons[5]->Fill(TMath::Min(TMath::Max(fabs(GenLeptons->At(idxMaxLep[0])->Eta()),
183     fabs(GenLeptons->At(idxMaxLep[1])->Eta())),4.999));
184     hDGenLeptons[6]->Fill(TMath::Min(TMath::Min(fabs(GenLeptons->At(idxMaxLep[0])->Eta()),
185     fabs(GenLeptons->At(idxMaxLep[1])->Eta())),4.999));
186     if(fabs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
187     fabs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5){
188     hDGenLeptons[7]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[0])->Pt(),199.999));
189     if(GenLeptons->At(idxMaxLep[0])->Pt() > 20.0){
190     hDGenLeptons[8]->Fill(TMath::Min(GenLeptons->At(idxMaxLep[1])->Pt(),199.999));
191     if(GenLeptons->At(idxMaxLep[1])->Pt() > 10.0){
192     CompositeParticle *dilepton = new CompositeParticle();
193     dilepton->AddDaughter(GenLeptons->At(idxMaxLep[0]));
194     dilepton->AddDaughter(GenLeptons->At(idxMaxLep[1]));
195     hDGenLeptons[9]->Fill(TMath::Min(dilepton->Mass(),999.999));
196     hDGenLeptons[10]->Fill(MathUtils::DeltaPhi(GenLeptons->At(idxMaxLep[0])->Phi(),
197     GenLeptons->At(idxMaxLep[1])->Phi())
198     * 180./ TMath::Pi());
199     delete dilepton;
200     }
201     }
202     }
203     }
204    
205     // Taus
206     hDGenTaus[0]->Fill(GenTaus->GetEntries());
207     for(UInt_t i=0; i<GenTaus->GetEntries(); i++){
208     hDGenTaus[1]->Fill(GenTaus->At(i)->Pt());
209     hDGenTaus[2]->Fill(GenTaus->At(i)->Eta());
210     hDGenTaus[3]->Fill(GenTaus->At(i)->Phi() * 180. / TMath::Pi());
211     }
212    
213     // Neutrinos
214     hDGenNeutrinos[0]->Fill(GenNeutrinos->GetEntries());
215     CompositeParticle *neutrinoTotal = new CompositeParticle();
216     for(UInt_t i=0; i<GenNeutrinos->GetEntries(); i++){
217     if(GenNeutrinos->At(i)->HasMother())
218     neutrinoTotal->AddDaughter(GenNeutrinos->At(i));
219     }
220     if(GenNeutrinos->GetEntries() > 0){
221     hDGenNeutrinos[1]->Fill(neutrinoTotal->Pt());
222     hDGenNeutrinos[2]->Fill(neutrinoTotal->Eta());
223     hDGenNeutrinos[3]->Fill(neutrinoTotal->Phi() * 180./ TMath::Pi());
224     }
225     delete neutrinoTotal;
226    
227     // Quarks
228     hDGenQuarks[0]->Fill(GenQuarks->GetEntries());
229     for(UInt_t i=0; i<GenQuarks->GetEntries(); i++){
230     for(UInt_t j=i+1; j<GenQuarks->GetEntries(); j++){
231     CompositeParticle *dijet = new CompositeParticle();
232     dijet->AddDaughter(GenQuarks->At(i));
233     dijet->AddDaughter(GenQuarks->At(j));
234     hDGenQuarks[1]->Fill(dijet->Pt());
235     hDGenQuarks[2]->Fill(dijet->Mass());
236     if(fabs(GenQuarks->At(i)->Eta()) < 2.5 && fabs(GenQuarks->At(j)->Eta()) < 2.5){
237     hDGenQuarks[3]->Fill(dijet->Pt());
238     hDGenQuarks[4]->Fill(dijet->Mass());
239     }
240     delete dijet;
241     }
242     }
243    
244     // WBF
245     if(GenqqHs->GetEntries() == 2){
246     hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
247     GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
248     hDGenWBF[1]->Fill(fabs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
249     hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
250     hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
251     CompositeParticle *diqq = new CompositeParticle();
252     diqq->AddDaughter(GenqqHs->At(0));
253     diqq->AddDaughter(GenqqHs->At(1));
254     hDGenWBF[4]->Fill(diqq->Mass());
255     delete diqq;
256     }
257    
258     // Bosons
259     hDGenBosons[0]->Fill(GenBosons->GetEntries());
260     for(UInt_t i=0; i<GenBosons->GetEntries(); i++){
261     hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
262     hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
263     hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
264     }
265     }
266    
267     }
268    
269     //--------------------------------------------------------------------------------------------------
270     void GeneratorMod::SlaveBegin()
271     {
272     // Run startup code on the computer (slave) doing the actual analysis. Here,
273     // we typically initialize histograms and other analysis objects and request
274     // branches. For this module, we request a branch of the MitTree.
275     ReqBranch(fMCPartName, fParticles);
276    
277     // Fill histograms
278     if(fFillHist == true){
279     char sb[200];
280     sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
281     sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
282     sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
283     sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
284     sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
285     sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
286     sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
287     sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
288     sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
289     sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
290     sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
291     for(int i=0; i<11; i++) AddOutput(hDGenLeptons[i]);
292    
293     // Taus
294     sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
295     sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
296     sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
297     sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
298     for(int i=0; i<4; i++) AddOutput(hDGenTaus[i]);
299    
300     // Neutrinos
301     sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
302     sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
303     sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
304     sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
305     for(int i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
306    
307     // Quarks
308     sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
309     sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
310     sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
311     sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
312     sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
313     for(int i=0; i<5; i++) AddOutput(hDGenQuarks[i]);
314    
315     // qqH
316     sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
317     sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
318     sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
319     sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
320     sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
321     for(int i=0; i<5; i++) AddOutput(hDGenWBF[i]);
322    
323     // Bosons
324     sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
325     sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
326     sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
327     sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
328     for(int i=0; i<4; i++) AddOutput(hDGenBosons[i]);
329     }
330     }
331    
332     //--------------------------------------------------------------------------------------------------
333     void GeneratorMod::SlaveTerminate()
334     {
335     // Run finishing code on the computer (slave) that did the analysis. For this
336     // module, we dont do anything here.
337    
338     }
339    
340     //--------------------------------------------------------------------------------------------------
341     void GeneratorMod::Terminate()
342     {
343     // Run finishing code on the client computer. For this module, we dont do
344     // anything here.
345     }