-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRunJetFinder2009emb_ue.C
297 lines (239 loc) · 11.7 KB
/
RunJetFinder2009emb_ue.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
//
// Pibero Djawotho <[email protected]>
// Texas A&M University
// 4 Dec 2011
//
void RunJetFinder2009emb_ue(int nevents = 1e6,
const char* mudstfile = "/gpfs01/star/scratch/gdwebb/Embedding_inclu/10103041/pt9_11_10103041_6.MuDst.root",
const char* geantfile = "/gpfs01/star/scratch/gdwebb/Embedding_inclu/10103041/pt9_11_10103041_6.geant.root",
const char* jetfile = "pt9_11_10103041_6.jets.root",
const char* skimfile = "pt9_11_10103041_6.skim.root",
const char* uefile = "pt9_11_10103041_6.ue_region.root",
const char* uefile2 = "pt9_11_10103041_6.ue_offcone.root",
bool useL2 = false)
{
cout << "Read MuDst file:\t" << mudstfile << endl;
cout << "Read geant file:\t" << geantfile << endl;
cout << "Write jet file:\t" << jetfile << endl;
cout << "Write skim file:\t" << skimfile << endl;
cout << "Write ue region file:\t" << uefile << endl;
cout << "Write ue off cone file:\t" << uefile2 << endl;
gROOT->Macro("loadMuDst.C");
gROOT->Macro("LoadLogger.C");
gSystem->Load("StDetectorDbMaker");
gSystem->Load("StTpcDb");
gSystem->Load("StDbUtilities");
gSystem->Load("StMcEvent");
gSystem->Load("StMcEventMaker");
gSystem->Load("StDaqLib");
gSystem->Load("StEmcRawMaker");
gSystem->Load("StEmcADCtoEMaker");
gSystem->Load("StPreEclMaker");
gSystem->Load("StEmcSimulatorMaker");
gSystem->Load("StDbBroker");
gSystem->Load("St_db_Maker");
gSystem->Load("StEEmcUtil");
gSystem->Load("StEEmcDbMaker");
gSystem->Load("StSpinDbMaker");
gSystem->Load("StEmcTriggerMaker");
gSystem->Load("StTriggerUtilities");
gSystem->Load("StMCAsymMaker");
gSystem->Load("StRandomSelector");
gSystem->Load("libfastjet.so");
gSystem->Load("libsiscone.so");
gSystem->Load("libsiscone_spherical.so");
gSystem->Load("libfastjetplugins.so");
gSystem->Load("StJetFinder");
gSystem->Load("StJetSkimEvent");
gSystem->Load("StJets");
gSystem->Load("StJetEvent");
gSystem->Load("StUeEvent");
gSystem->Load("StJetMaker");
gSystem->Load("StEEmcSimulatorMaker");
// Create chain
StChain* chain = new StChain;
// I/O maker
StIOMaker* ioMaker = new StIOMaker;
ioMaker->SetFile(geantfile);
ioMaker->SetIOMode("r");
ioMaker->SetBranch("*",0,"0"); // Deactivate all branches
ioMaker->SetBranch("geantBranch",0,"r"); // Activate geant Branch
// StMcEvent maker
StMcEventMaker* mcEventMaker = new StMcEventMaker;
mcEventMaker->doPrintEventInfo = false;
mcEventMaker->doPrintMemoryInfo = false;
// MuDst reader
StMuDstMaker* muDstMaker = new StMuDstMaker(0,0,"",mudstfile,"",100000,"MuDst");
// MuDst DB
StMuDbReader* muDstDb = StMuDbReader::instance();
// star database
St_db_Maker* starDb = new St_db_Maker("StarDb","MySQL:StarDb");
// Endcap database
StEEmcDbMaker* eemcDb = new StEEmcDbMaker;
// Barrel ADC to energy maker
StEmcADCtoEMaker* adc2e = new StEmcADCtoEMaker;
adc2e->saveAllStEvent(true);
// Trigger simulator
// -- StL2_2009EmulatorMaker must run before StTriggerSimuMaker?
StL2_2009EmulatorMaker* simL2Mk = 0;
if (useL2) {
simL2Mk = new StL2_2009EmulatorMaker;
simL2Mk->setSetupPath("/star/u/pibero/public/StarTrigSimuSetup/");
simL2Mk->setOutPath("/star/institutions/uky/gdwebb/2009/500EmbedJets/L2/");
}
StTriggerSimuMaker* simuTrig = new StTriggerSimuMaker;
// Pibero, I changed the next line from 2 to 1.
simuTrig->useOnlineDB();
simuTrig->setMC(2); // 0=data, 1=simulation, 2=embedding
simuTrig->useBemc();
simuTrig->useEemc();
simuTrig->bemc->setConfig(StBemcTriggerSimu::kOffline);
if (useL2) simuTrig->useL2(simL2Mk);
// Set trigger thresholds using run 10102039
simuTrig->setBarrelJetPatchTh(0,32);
simuTrig->setBarrelJetPatchTh(1,40);
simuTrig->setBarrelJetPatchTh(2,60);
simuTrig->setOverlapJetPatchTh(0,32);
simuTrig->setOverlapJetPatchTh(1,40);
simuTrig->setOverlapJetPatchTh(2,60);
simuTrig->setEndcapJetPatchTh(0,32);
simuTrig->setEndcapJetPatchTh(1,40);
simuTrig->setEndcapJetPatchTh(2,60);
simuTrig->setBarrelHighTowerTh(0,-1);
simuTrig->setBarrelHighTowerTh(1,-1);
simuTrig->setBarrelHighTowerTh(2,25);
simuTrig->setBarrelHighTowerTh(3,31);
simuTrig->setEndcapHighTowerTh(0,25);
simuTrig->setEndcapHighTowerTh(1,31);
// Get Pythia record
StMCAsymMaker* asym = new StMCAsymMaker;
// Skim event maker
StJetSkimEventMaker* skimEventMaker = new StJetSkimEventMaker("StJetSkimEventMaker",muDstMaker,skimfile);
skimEventMaker->addSimuTrigger(230531); // BHT3
skimEventMaker->addSimuTrigger(230410); // JP1
skimEventMaker->addSimuTrigger(230411); // JP2
skimEventMaker->addSimuTrigger(230420); // AJP
// Jet maker
StJetMaker2012* jetmaker = new StJetMaker2012;
jetmaker->setJetFile(jetfile);
jetmaker->setJetFileUe(uefile2);
// UE maker
StUEMaker2009* uemaker = new StUEMaker2009;
uemaker->setUeFile(uefile);
//------------------------------------------------------------------------------------
// Set analysis cuts for 12-point branch
StAnaPars* anapars12 = new StAnaPars;
anapars12->useTpc = true;
anapars12->useBemc = true;
anapars12->useEemc = false;
anapars12->randomSelectorProb = 1.00;
// The classes available for correcting tower energy for tracks are:
// 1. StjTowerEnergyCorrectionForTracksMip
// 2. StjTowerEnergyCorrectionForTracksFraction
// 3. StjTowerEnergyCorrectionForTracksNull (default: no correction
//The classes available for shifting track and tower parameters prior to jet reconstruction are:
// 1. StjTrackPtFraction
// 2. StjTowerEnergyFraction
// The input parameter for is the fraction which you want to increase or decrease the track/tower quantity.
// i.e. If you want to decrease the track pT by 5%, you should put in -0.05 as the input parameter for StjTrackPtFraction
anapars12->setTrackShift(new StjTrackPtFraction(0.0));
anapars12->setTowerShift(new StjTowerEnergyFraction(0.0));
anapars12->setTowerEnergyCorrection(new StjTowerEnergyCorrectionForTracksFraction(1.00));
// TPC cuts
anapars12->addTpcCut(new StjTrackCutFlag(0));
anapars12->addTpcCut(new StjTrackCutNHits(12));
anapars12->addTpcCut(new StjTrackCutPossibleHitRatio(0.51));
anapars12->addTpcCut(new StjTrackCutDca(3));
anapars12->addTpcCut(new StjTrackCutTdcaPtDependent);
anapars12->addTpcCut(new StjTrackCutPt(0.2,200));
anapars12->addTpcCut(new StjTrackCutEta(-2.5,2.5));
anapars12->addTpcCut(new StjTrackCutLastPoint(125));
// BEMC cuts
anapars12->addBemcCut(new StjTowerEnergyCutBemcStatus(1));
anapars12->addBemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
anapars12->addBemcCut(new StjTowerEnergyCutEt(0.2));
// EEMC cuts
anapars12->addEemcCut(new StjTowerEnergyCutBemcStatus(1));
anapars12->addEemcCut(new StjTowerEnergyCutAdc(4,3)); // ADC-ped>4 AND ADC-ped>3*RMS
anapars12->addEemcCut(new StjTowerEnergyCutEt(0.2));
// Jet cuts
anapars12->addJetCut(new StProtoJetCutPt(5,200));
anapars12->addJetCut(new StProtoJetCutEta(-100,100));
// _________________ Regional Criteria_____________________
StAnaPars *anapars12_tow = new StAnaPars(*anapars12);
anapars12_tow->setTrackRegion(new StjTrackRegion(60.0,-60.0,1.0));
anapars12_tow->setTowerRegion(new StjTowerRegion(60.0,-60.0,1.0));
StAnaPars *anapars12_away = new StAnaPars(*anapars12);
anapars12_away->setTrackRegion(new StjTrackRegion(120.0,-120.0,1.0));
anapars12_away->setTowerRegion(new StjTowerRegion(120.0,-120.0,1.0));
StAnaPars *anapars12_transPlus = new StAnaPars(*anapars12);
anapars12_transPlus->setTrackRegion(new StjTrackRegion(120.0,60.0,1.0));
anapars12_transPlus->setTowerRegion(new StjTowerRegion(120.0,60.0,1.0));
StAnaPars *anapars12_transMinus = new StAnaPars(*anapars12);
anapars12_transMinus->setTrackRegion(new StjTrackRegion(-60.0,-120.0,1.0));
anapars12_transMinus->setTowerRegion(new StjTowerRegion(-60.0,-120.0,1.0));
// -----------Set analysis cuts for particle jets branch---------------------
StAnaPars* anaparsParticle = new StAnaPars;
anaparsParticle->useMonteCarlo = true;
// MC cuts
anaparsParticle->addMcCut(new StjMCParticleCutStatus(1)); // final state particles
// Jet cuts
anaparsParticle->addJetCut(new StProtoJetCutPt(3,200));
anaparsParticle->addJetCut(new StProtoJetCutEta(-100,100));
StAnaPars* anaparsParticle_toward = new StAnaPars(*anaparsParticle);
anaparsParticle_toward->setParticleRegion(new StjMCParticleRegion(60.0,-60.0,1.0));
StAnaPars* anaparsParticle_away = new StAnaPars(*anaparsParticle);
anaparsParticle_away->setParticleRegion(new StjMCParticleRegion(120.0,-120.0,1.0));
StAnaPars* anaparsParticle_transPlus = new StAnaPars(*anaparsParticle);
anaparsParticle_transPlus->setParticleRegion(new StjMCParticleRegion(120.0,60.0,1.0));
StAnaPars* anaparsParticle_transMinus = new StAnaPars(*anaparsParticle);
anaparsParticle_transMinus->setParticleRegion(new StjMCParticleRegion(-60.0,-120.0,1.0));
//------------------------------------------------------------------------------------
// Set analysis cuts for parton jets branch
StAnaPars* anaparsParton = new StAnaPars;
anaparsParton->useMonteCarlo = true;
// MC cuts
anaparsParton->addMcCut(new StjMCParticleCutParton);
// Jet cuts
anaparsParton->addJetCut(new StProtoJetCutPt(3,200));
anaparsParton->addJetCut(new StProtoJetCutEta(-100,100));
StAnaPars* anaparsParton_toward = new StAnaPars(*anaparsParton);
anaparsParton_toward->setParticleRegion(new StjMCParticleRegion(60.0,-60.0,1.0));
StAnaPars* anaparsParton_away = new StAnaPars(*anaparsParton);
anaparsParton_away->setParticleRegion(new StjMCParticleRegion(120.0,-120.0,1.0));
StAnaPars* anaparsParton_transPlus = new StAnaPars(*anaparsParton);
anaparsParton_transPlus->setParticleRegion(new StjMCParticleRegion(120.0,60.0,1.0));
StAnaPars* anaparsParton_transMinus = new StAnaPars(*anaparsParton);
anaparsParton_transMinus->setParticleRegion(new StjMCParticleRegion(-60.0,-120.0,1.0));
//--------------------------------------------------------------------------------------------
StFastJetAreaPars *JetAreaPars = new StFastJetAreaPars;
JetAreaPars->setGhostArea(0.01);
// Set anti-kt R=0.6 parameters
StFastJetPars* AntiKtR060Pars = new StFastJetPars;
AntiKtR060Pars->setJetAlgorithm(StFastJetPars::antikt_algorithm);
AntiKtR060Pars->setRparam(0.6);
AntiKtR060Pars->setRecombinationScheme(StFastJetPars::E_scheme);
AntiKtR060Pars->setStrategy(StFastJetPars::Best);
AntiKtR060Pars->setPtMin(5.0);
AntiKtR060Pars->setJetArea(JetAreaPars);
jetmaker->addBranch("AntiKtR060Particle",anaparsParticle,AntiKtR060Pars);
jetmaker->addBranch("AntiKtR060Parton",anaparsParton,AntiKtR060Pars);
jetmaker->addBranch("AntiKtR060NHits12",anapars12,AntiKtR060Pars);
uemaker->addBranch("toward",anapars12_tow,"AntiKtR060NHits12");
uemaker->addBranch("away",anapars12_away,"AntiKtR060NHits12");
uemaker->addBranch("transP",anapars12_transPlus,"AntiKtR060NHits12");
uemaker->addBranch("transM",anapars12_transMinus,"AntiKtR060NHits12");
uemaker->addBranch("toward",anaparsParticle_toward,"AntiKtR060Particle");
uemaker->addBranch("away",anaparsParticle_away,"AntiKtR060Particle");
uemaker->addBranch("transP",anaparsParticle_transPlus,"AntiKtR060Particle");
uemaker->addBranch("transM",anaparsParticle_transMinus,"AntiKtR060Particle");
uemaker->addBranch("toward",anaparsParton_toward,"AntiKtR060Parton");
uemaker->addBranch("away",anaparsParton_away,"AntiKtR060Parton");
uemaker->addBranch("transP",anaparsParton_transPlus,"AntiKtR060Parton");
uemaker->addBranch("transM",anaparsParton_transMinus,"AntiKtR060Parton");
StOffAxisConesPars *off060 = new StOffAxisConesPars(0.6);
jetmaker->addUeBranch("OffAxisConesR060", off060);
//------------------------------------------------------------------------------------
chain->Init();
chain->EventLoop(nevents);
}