ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/GeneratorMod.cc
Revision: 1.2
Committed: Thu Oct 23 12:23:44 2008 UTC (16 years, 6 months ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.1: +42 -2 lines
Log Message:
changes after new life

File Contents

# User Rev Content
1 ceballos 1.2 // $Id: GeneratorMod.cc,v 1.5 2008/10/10 10:54:13 ceballos Exp $
2 loizides 1.1
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 ceballos 1.2 // b quark info
243     if (GenQuarks->At(i)->AbsPdgId() == 5){
244     hDGenQuarks[5]->Fill(GenQuarks->At(i)->Pt());
245     hDGenQuarks[6]->Fill(GenQuarks->At(i)->Eta());
246     hDGenQuarks[7]->Fill(GenQuarks->At(i)->Phi());
247     if(ptMaxLep[0] > 20 && ptMaxLep[1] > 15){
248     if(fabs(GenLeptons->At(idxMaxLep[0])->Eta()) < 2.5 &&
249     fabs(GenLeptons->At(idxMaxLep[1])->Eta()) < 2.5){
250     hDGenQuarks[8]->Fill(GenQuarks->At(i)->Pt());
251     hDGenQuarks[9]->Fill(GenQuarks->At(i)->Eta());
252     hDGenQuarks[10]->Fill(GenQuarks->At(i)->Phi());
253     }
254     }
255     }
256     // t quark info
257     else if(GenQuarks->At(i)->AbsPdgId() == 6){
258     hDGenQuarks[11]->Fill(GenQuarks->At(i)->Pt());
259     hDGenQuarks[12]->Fill(GenQuarks->At(i)->Eta());
260     hDGenQuarks[13]->Fill(GenQuarks->At(i)->Phi());
261     }
262     // light quark info
263     else {
264     hDGenQuarks[14]->Fill(GenQuarks->At(i)->Pt());
265     hDGenQuarks[15]->Fill(GenQuarks->At(i)->Eta());
266     hDGenQuarks[16]->Fill(GenQuarks->At(i)->Phi());
267     }
268 loizides 1.1 }
269    
270     // WBF
271     if(GenqqHs->GetEntries() == 2){
272     hDGenWBF[0]->Fill(MathUtils::DeltaPhi(GenqqHs->At(0)->Phi(),
273     GenqqHs->At(1)->Phi()) * 180./ TMath::Pi());
274     hDGenWBF[1]->Fill(fabs(GenqqHs->At(0)->Eta()-GenqqHs->At(1)->Eta()));
275     hDGenWBF[2]->Fill(TMath::Max(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
276     hDGenWBF[3]->Fill(TMath::Min(GenqqHs->At(0)->Pt(),GenqqHs->At(1)->Pt()));
277     CompositeParticle *diqq = new CompositeParticle();
278     diqq->AddDaughter(GenqqHs->At(0));
279     diqq->AddDaughter(GenqqHs->At(1));
280     hDGenWBF[4]->Fill(diqq->Mass());
281     delete diqq;
282     }
283    
284     // Bosons
285     hDGenBosons[0]->Fill(GenBosons->GetEntries());
286     for(UInt_t i=0; i<GenBosons->GetEntries(); i++){
287     hDGenBosons[1]->Fill(GenBosons->At(i)->Pt());
288     hDGenBosons[2]->Fill(GenBosons->At(i)->Eta());
289     hDGenBosons[3]->Fill(GenBosons->At(i)->Mass());
290     }
291     }
292    
293     }
294    
295     //--------------------------------------------------------------------------------------------------
296     void GeneratorMod::SlaveBegin()
297     {
298     // Run startup code on the computer (slave) doing the actual analysis. Here,
299     // we typically initialize histograms and other analysis objects and request
300     // branches. For this module, we request a branch of the MitTree.
301     ReqBranch(fMCPartName, fParticles);
302    
303     // Fill histograms
304     if(fFillHist == true){
305     char sb[200];
306     sprintf(sb,"hDGenLeptons_%d", 0); hDGenLeptons[0] = new TH1D(sb,sb,10,-0.5,9.5);
307     sprintf(sb,"hDGenLeptons_%d", 1); hDGenLeptons[1] = new TH1D(sb,sb,100,0.0,200.0);
308     sprintf(sb,"hDGenLeptons_%d", 2); hDGenLeptons[2] = new TH1D(sb,sb,50,0.0,5.0);
309     sprintf(sb,"hDGenLeptons_%d", 3); hDGenLeptons[3] = new TH1D(sb,sb,90,0.0,180.0);
310     sprintf(sb,"hDGenLeptons_%d", 4); hDGenLeptons[4] = new TH1D(sb,sb,1000,0.0,1000.0);
311     sprintf(sb,"hDGenLeptons_%d", 5); hDGenLeptons[5] = new TH1D(sb,sb,50,0.0,5.0);
312     sprintf(sb,"hDGenLeptons_%d", 6); hDGenLeptons[6] = new TH1D(sb,sb,50,0.0,5.0);
313     sprintf(sb,"hDGenLeptons_%d", 7); hDGenLeptons[7] = new TH1D(sb,sb,100,0.0,200.0);
314     sprintf(sb,"hDGenLeptons_%d", 8); hDGenLeptons[8] = new TH1D(sb,sb,100,0.0,200.0);
315     sprintf(sb,"hDGenLeptons_%d", 9); hDGenLeptons[9] = new TH1D(sb,sb,1000,0.0,1000.0);
316     sprintf(sb,"hDGenLeptons_%d",10); hDGenLeptons[10] = new TH1D(sb,sb,90,0.0,180.0);
317     for(int i=0; i<11; i++) AddOutput(hDGenLeptons[i]);
318    
319     // Taus
320     sprintf(sb,"hDGenTaus_%d", 0); hDGenTaus[0] = new TH1D(sb,sb,10,-0.5,9.5);
321     sprintf(sb,"hDGenTaus_%d", 1); hDGenTaus[1] = new TH1D(sb,sb,100,0.0,200.0);
322     sprintf(sb,"hDGenTaus_%d", 2); hDGenTaus[2] = new TH1D(sb,sb,100,-5.0,5.0);
323     sprintf(sb,"hDGenTaus_%d", 3); hDGenTaus[3] = new TH1D(sb,sb,90,0.0,180.0);
324     for(int i=0; i<4; i++) AddOutput(hDGenTaus[i]);
325    
326     // Neutrinos
327     sprintf(sb,"hDGenNeutrinos_%d", 0); hDGenNeutrinos[0] = new TH1D(sb,sb,10,-0.5,9.5);
328     sprintf(sb,"hDGenNeutrinos_%d", 1); hDGenNeutrinos[1] = new TH1D(sb,sb,100,0.0,200.0);
329     sprintf(sb,"hDGenNeutrinos_%d", 2); hDGenNeutrinos[2] = new TH1D(sb,sb,100,-5.0,5.0);
330     sprintf(sb,"hDGenNeutrinos_%d", 3); hDGenNeutrinos[3] = new TH1D(sb,sb,90,0.0,180.0);
331     for(int i=0; i<4; i++) AddOutput(hDGenNeutrinos[i]);
332    
333     // Quarks
334     sprintf(sb,"hDGenQuarks_%d", 0); hDGenQuarks[0] = new TH1D(sb,sb,10,-0.5,9.5);
335     sprintf(sb,"hDGenQuarks_%d", 1); hDGenQuarks[1] = new TH1D(sb,sb,200,0.0,400.);
336     sprintf(sb,"hDGenQuarks_%d", 2); hDGenQuarks[2] = new TH1D(sb,sb,2000,0.0,2000.);
337     sprintf(sb,"hDGenQuarks_%d", 3); hDGenQuarks[3] = new TH1D(sb,sb,200,0.0,400.);
338     sprintf(sb,"hDGenQuarks_%d", 4); hDGenQuarks[4] = new TH1D(sb,sb,2000,0.0,2000.);
339 ceballos 1.2 sprintf(sb,"hDGenQuarks_%d", 5); hDGenQuarks[5] = new TH1D(sb,sb,200,0.0,400.);
340     sprintf(sb,"hDGenQuarks_%d", 6); hDGenQuarks[6] = new TH1D(sb,sb,200,-10.0,10.);
341     sprintf(sb,"hDGenQuarks_%d", 7); hDGenQuarks[7] = new TH1D(sb,sb,200,-TMath::Pi(),TMath::Pi());
342     sprintf(sb,"hDGenQuarks_%d", 8); hDGenQuarks[8] = new TH1D(sb,sb,200,0.0,400.);
343     sprintf(sb,"hDGenQuarks_%d", 9); hDGenQuarks[9] = new TH1D(sb,sb,200,-10.0,10.);
344     sprintf(sb,"hDGenQuarks_%d",10); hDGenQuarks[10] = new TH1D(sb,sb,200,-TMath::Pi(),TMath::Pi());
345     sprintf(sb,"hDGenQuarks_%d",11); hDGenQuarks[11] = new TH1D(sb,sb,200,0.0,400.);
346     sprintf(sb,"hDGenQuarks_%d",12); hDGenQuarks[12] = new TH1D(sb,sb,200,-10.0,10.);
347     sprintf(sb,"hDGenQuarks_%d",13); hDGenQuarks[13] = new TH1D(sb,sb,200,-TMath::Pi(),TMath::Pi());
348     sprintf(sb,"hDGenQuarks_%d",14); hDGenQuarks[14] = new TH1D(sb,sb,200,0.0,400.);
349     sprintf(sb,"hDGenQuarks_%d",15); hDGenQuarks[15] = new TH1D(sb,sb,200,-10.0,10.);
350     sprintf(sb,"hDGenQuarks_%d",16); hDGenQuarks[16] = new TH1D(sb,sb,200,-TMath::Pi(),TMath::Pi());
351     for(int i=0; i<17; i++) AddOutput(hDGenQuarks[i]);
352 loizides 1.1
353     // qqH
354     sprintf(sb,"hDGenWBF_%d", 0); hDGenWBF[0] = new TH1D(sb,sb,90,0.0,180.);
355     sprintf(sb,"hDGenWBF_%d", 1); hDGenWBF[1] = new TH1D(sb,sb,100,0.0,10.);
356     sprintf(sb,"hDGenWBF_%d", 2); hDGenWBF[2] = new TH1D(sb,sb,200,0.0,400.);
357     sprintf(sb,"hDGenWBF_%d", 3); hDGenWBF[3] = new TH1D(sb,sb,200,0.0,400.);
358     sprintf(sb,"hDGenWBF_%d", 4); hDGenWBF[4] = new TH1D(sb,sb,200,0.0,4000.);
359     for(int i=0; i<5; i++) AddOutput(hDGenWBF[i]);
360    
361     // Bosons
362     sprintf(sb,"hDGenBosons_%d", 0); hDGenBosons[0] = new TH1D(sb,sb,10,-0.5,9.5);
363     sprintf(sb,"hDGenBosons_%d", 1); hDGenBosons[1] = new TH1D(sb,sb,200,0.0,400.0);
364     sprintf(sb,"hDGenBosons_%d", 2); hDGenBosons[2] = new TH1D(sb,sb,100,-5.0,5.0);
365     sprintf(sb,"hDGenBosons_%d", 3); hDGenBosons[3] = new TH1D(sb,sb,2000,0.0,2000.0);
366     for(int i=0; i<4; i++) AddOutput(hDGenBosons[i]);
367     }
368     }
369    
370     //--------------------------------------------------------------------------------------------------
371     void GeneratorMod::SlaveTerminate()
372     {
373     // Run finishing code on the computer (slave) that did the analysis. For this
374     // module, we dont do anything here.
375    
376 ceballos 1.2 cout << "Total Number of Processed Events: " << fNEventsProcessed << endl;
377    
378 loizides 1.1 }
379    
380     //--------------------------------------------------------------------------------------------------
381     void GeneratorMod::Terminate()
382     {
383     // Run finishing code on the client computer. For this module, we dont do
384     // anything here.
385     }