Attachment 'eventSelection.c'
Download 1 #include <math.h>
2 #include <TFile.h>
3 #include <TTree.h>
4 #include <TGraph.h>
5 #include<TString.h>
6 #include<iostream>
7
8 /******** eventSelection processor ********/
9 /* */
10 /* reads nFiles root trees, and applies */
11 /* cuts defined in selection criteria */
12 /* result is written to outTree */
13 /* */
14 /******************************************/
15
16 const int nFiles = 1;
17 const int MAXCLUSTER = 30;
18
19 const float E_cms = 500.0; //center of mass energy [GeV]
20
21 //selection criteria
22 const int _minNhitsEcal = 1500; //minimum # hits in ECAL: 1500
23 const int _maxNhitsEcal = 6000; //maximum # hits in ECAL: 6000
24 const float _EminEcal = 0; //minimum energy in ECAL: 80 GeV
25 const float _EmaxEcal = 700; //maximum energy in ECAL: 450 GeV
26 const float _minEt_miss = 0; //minimum transverse energy: 30 GeV
27 const float _minNgamma = 0; //minimum # of cluster: 2
28 const float _minEcluster = 0; //minimum energy of cluster: 20 GeV
29 const float _maxCosTheta = 1.00; //maximum cluster theta position: 0.75
30
31 using namespace std;
32
33 void eventSelection()
34 {
35 TFile *inFile;
36 TString inFileNames[nFiles] = {
37 "/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa.root"
38 };
39
40 // float nExpected[nFiles] = {216728,3113697,2001046,70319,70674};
41 float nExpected[nFiles] = {70319};
42 //500fb = 216728 ee_neu1neu1
43 // 2001046 ee_n1n1aa
44 // 70319 ee_n2n2aa
45 // 70674 ee_n3n3aa
46 // 3113697 ee_aa
47 // 180211 ee_eaea
48 // 73715739 ee_ee
49 float polWeight[nFiles] = {1.};//,1.,1.,1.};
50 //0.02 ep+1.0em-1.0, 0.72 ep-1.0em+1.0
51 //0.40 ep-1.0em+0.0, 0.10 ep+1.0em+0.0, 0.05 ep+0.0em-1.0, 0.45 ep+0.0em+1.0
52 float weight[nFiles]; // nExpected*polWeight/nTrue;
53
54 /******** create output tree **********/
55 struct {
56 int iEvt; // event number
57 int nCluster; // # cluster per event
58 int nGamma; // # photons per event
59 float etMiss; // missing transverse energy [GeV]
60 float esumCalos; // energy sum in calorimeters [GeV]
61 int nHitsCalos; // # hits in calorimeters
62 float esumECAL; // energy sum in ECAL [GeV]
63 int nHitsECAL; // # hits in ECAL
64 float esumHCAL; // energy sum in HCAL [GeV]
65 int nHitsHCAL; // # hits in HCAL
66 float energy[MAXCLUSTER]; // energy per cluster [GeV]
67 float clusterPos[MAXCLUSTER][3]; // cluster center of gravity
68 float clusterAxis[MAXCLUSTER][3]; // cluster axis of inertia
69 float thetaDir[MAXCLUSTER]; // cluster angle theta [degree]
70 float phiDir[MAXCLUSTER]; // cluster angle phi [degree]
71 float cosThetaPos[MAXCLUSTER]; // cosine of cluster position theta
72 float phiPos[MAXCLUSTER]; // cluster position phi [degree]
73 float dca[MAXCLUSTER]; // distance of closest approach to IP [mm]
74 int nHitsPerCluster[MAXCLUSTER]; // #hits per cluster
75 int trueNcluster; // true # cluster / event
76 float trueDist[MAXCLUSTER]; // true distance of MCparticle to IP [mm]
77 float trueGammaEnergy[MAXCLUSTER]; // true energy of MCparticle gamma [GeV]
78 float trueGammaDirX[MAXCLUSTER]; // true x direction of MCparticle gamma [mm]
79 float trueGammaDirY[MAXCLUSTER]; // true y direction of MCparticle gamma [mm]
80 float trueGammaDirZ[MAXCLUSTER]; // true z direction of MCparticle gamma [mm]
81 float trueGammaPosX[MAXCLUSTER]; // true x position of MCparticle gamma [mm]
82 float trueGammaPosY[MAXCLUSTER]; // true y position of MCparticle gamma [mm]
83 float trueGammaPosZ[MAXCLUSTER]; // true z position of MCparticle gamma [mm]
84 float trueNeutralinoEnergy[MAXCLUSTER]; // true energy of MCparticle neu1 [GeV]
85 float trueNeutralinoDirX[MAXCLUSTER]; // true x direction of MCparticle neu1 [mm]
86 float trueNeutralinoDirY[MAXCLUSTER]; // true y direction of MCparticle neu1 [mm]
87 float trueNeutralinoDirZ[MAXCLUSTER]; // true z direction of MCparticle neu1 [mm]
88 float trueThetaDir[MAXCLUSTER]; // true gamma angle theta [degree]
89 float truePhiDir[MAXCLUSTER]; // true gamma angle phi [degree]
90 } _clusterFill;
91
92 TFile *outFile = new TFile("/data/nwattime/LDC_optimization/rootTrees/ILD00/ee_n3n3aa_Nhits.root","RECREATE");
93 TTree *outTree = new TTree("selectedTree","selectedTree");
94
95 outTree->Branch("eventNumber", &_clusterFill.iEvt, "eventNumber/I");
96 outTree->Branch("ECAL_energySum", &_clusterFill.esumECAL, "ECAL_energySum/F");
97 outTree->Branch("ECAL_nHits", &_clusterFill.nHitsECAL, "ECAL_nHits/I");
98 outTree->Branch("HCAL_energySum", &_clusterFill.esumHCAL, "HCAL_energySum/F");
99 outTree->Branch("HCAL_nHits", &_clusterFill.nHitsHCAL, "HCAL_nHits/I");
100 outTree->Branch("calos_energySum", &_clusterFill.esumCalos, "calos_energySum/F");
101 outTree->Branch("calos_nHits", &_clusterFill.nHitsCalos, "calos_nHits/I");
102 outTree->Branch("cluster_nCluster", &_clusterFill.nCluster, "cluster_nCluster/I");
103 outTree->Branch("cluster_nGamma", &_clusterFill.nGamma, "cluster_nGamma/I");
104 outTree->Branch("cluster_cosThetaPos", &_clusterFill.cosThetaPos,"cluster_cosThetaPos[cluster_nCluster]/F");
105 outTree->Branch("cluster_dca", &_clusterFill.dca, "cluster_dca[cluster_nCluster]/F");
106 outTree->Branch("cluster_energy", &_clusterFill.energy, "cluster_energy[cluster_nCluster]/F");
107 outTree->Branch("cluster_etMiss", &_clusterFill.etMiss, "cluster_etMiss/F");
108 outTree->Branch("cluster_nHits", &_clusterFill.nHitsPerCluster, "cluster_nHits[cluster_nCluster]/I");
109 outTree->Branch("cluster_phiDir", &_clusterFill.phiDir, "cluster_phiDir[cluster_nCluster]/F");
110 outTree->Branch("cluster_phiPos", &_clusterFill.phiPos, "cluster_phiPos[cluster_nCluster]/F");
111 outTree->Branch("cluster_position", &_clusterFill.clusterPos, "cluster_position[cluster_nCluster][3]/F");
112 outTree->Branch("cluster_clusterAxis", &_clusterFill.clusterAxis,"cluster_axisOfInertia[cluster_nCluster][3]/F");
113 outTree->Branch("cluster_thetaDir", &_clusterFill.thetaDir, "cluster_thetaDir[cluster_nCluster]/F");
114 outTree->Branch("cluster_trueNcluster", &_clusterFill.trueNcluster, "cluster_trueNcluster/I");
115 outTree->Branch("cluster_trueDist", &_clusterFill.trueDist, "cluster_trueDist[cluster_nCluster]/F");
116 outTree->Branch("cluster_trueGammaEnergy", &_clusterFill.trueGammaEnergy, "cluster_trueGammaEnergy[cluster_nCluster]/F");
117 outTree->Branch("cluster_trueGammaDirX", &_clusterFill.trueGammaDirX, "cluster_trueGammaDirX[cluster_nCluster]/F");
118 outTree->Branch("cluster_trueGammaDirY", &_clusterFill.trueGammaDirY, "cluster_trueGammaDirY[cluster_nCluster]/F");
119 outTree->Branch("cluster_trueGammaDirZ", &_clusterFill.trueGammaDirZ, "cluster_trueGammaDirZ[cluster_nCluster]/F");
120 outTree->Branch("cluster_trueGammaPosX", &_clusterFill.trueGammaPosX, "cluster_trueGammaPosX[cluster_nCluster]/F");
121 outTree->Branch("cluster_trueGammaPosY", &_clusterFill.trueGammaPosY, "cluster_trueGammaPosY[cluster_nCluster]/F");
122 outTree->Branch("cluster_trueGammaPosZ", &_clusterFill.trueGammaPosZ, "cluster_trueGammaPosZ[cluster_nCluster]/F");
123 outTree->Branch("cluster_trueNeutralinoEnergy", &_clusterFill.trueNeutralinoEnergy, "cluster_trueNeutralinoEnergy[cluster_nCluster]/F");
124 outTree->Branch("cluster_trueNeutralinoDirX", &_clusterFill.trueNeutralinoDirX, "cluster_trueNeutralinoDirX[cluster_nCluster]/F");
125 outTree->Branch("cluster_trueNeutralinoDirY", &_clusterFill.trueNeutralinoDirY, "cluster_trueNeutralinoDirY[cluster_nCluster]/F");
126 outTree->Branch("cluster_trueNeutralinoDirZ", &_clusterFill.trueNeutralinoDirZ, "cluster_trueNeutralinoDirZ[cluster_nCluster]/F");
127 outTree->Branch("cluster_truePhiDir", &_clusterFill.truePhiDir, "cluster_truePhiDir[cluster_nCluster]/F");
128 outTree->Branch("cluster_trueThetaDir", &_clusterFill.trueThetaDir, "cluster_trueThetaDir[cluster_nCluster]/F");
129
130 for (int fileCounter = 0; fileCounter<nFiles; ++fileCounter) {
131
132 inFile = new TFile(inFileNames[fileCounter],"READ");
133
134 if ( inFile->IsOpen() ) {
135
136 cout<<"opened file "<<fileCounter<<": "<<inFileNames[fileCounter]<<endl;
137
138 /******** read input tree **********/
139 int _iEvt;
140 int _nCluster; int _nGamma; int _trueNcluster;
141 float _etMiss;
142 float _esumCalos; float _esumECAL; float _esumHCAL;
143 int _nHitsCalos; int _nHitsECAL; int _nHitsHCAL;
144 float _energy[MAXCLUSTER]; int _nHitsPerCluster[MAXCLUSTER];
145 float _clusterPos[MAXCLUSTER][3]; float _clusterAxis[MAXCLUSTER][3];
146 float _thetaDir[MAXCLUSTER]; float _cosThetaPos[MAXCLUSTER];
147 float _phiDir[MAXCLUSTER]; float _phiPos[MAXCLUSTER];
148 float _dca[MAXCLUSTER];
149 float _trueDist[MAXCLUSTER];
150 float _trueGammaEnergy[MAXCLUSTER]; float _trueNeutralinoEnergy[MAXCLUSTER];
151 float _trueGammaDirX[MAXCLUSTER]; float _trueGammaDirY[MAXCLUSTER]; float _trueGammaDirZ[MAXCLUSTER];
152 float _trueGammaPosX[MAXCLUSTER]; float _trueGammaPosY[MAXCLUSTER]; float _trueGammaPosZ[MAXCLUSTER];
153 float _trueNeutralinoDirX[MAXCLUSTER]; float _trueNeutralinoDirY[MAXCLUSTER]; float _trueNeutralinoDirZ[MAXCLUSTER];
154 float _trueThetaDir[MAXCLUSTER]; float _truePhiDir[MAXCLUSTER];
155
156 TTree* inTree = (TTree*)inFile->Get("tree");
157 inTree->SetBranchAddress("eventNumber", &_iEvt);
158 inTree->SetBranchAddress("ECAL_energySum", &_esumECAL);
159 inTree->SetBranchAddress("ECAL_nHits", &_nHitsECAL);
160 inTree->SetBranchAddress("HCAL_energySum", &_esumHCAL);
161 inTree->SetBranchAddress("HCAL_nHits", &_nHitsHCAL);
162 inTree->SetBranchAddress("calos_energySum", &_esumCalos);
163 inTree->SetBranchAddress("calos_nHits", &_nHitsCalos);
164 inTree->SetBranchAddress("cluster_cosThetaPos", &_cosThetaPos);
165 inTree->SetBranchAddress("cluster_dca", &_dca);
166 inTree->SetBranchAddress("cluster_energy", &_energy);
167 inTree->SetBranchAddress("cluster_etMiss", &_etMiss);
168 inTree->SetBranchAddress("cluster_nCluster", &_nCluster);
169 inTree->SetBranchAddress("cluster_nGamma", &_nGamma);
170 inTree->SetBranchAddress("cluster_nHits", &_nHitsPerCluster);
171 inTree->SetBranchAddress("cluster_phiDir", &_phiDir);
172 inTree->SetBranchAddress("cluster_phiPos", &_phiPos);
173 inTree->SetBranchAddress("cluster_position", &_clusterPos);
174 inTree->SetBranchAddress("cluster_clusterAxis", &_clusterAxis);
175 inTree->SetBranchAddress("cluster_thetaDir", &_thetaDir);
176 inTree->SetBranchAddress("cluster_trueNcluster", &_trueNcluster);
177 inTree->SetBranchAddress("cluster_trueDist", &_trueDist);
178 inTree->SetBranchAddress("cluster_trueGammaEnergy", &_trueGammaEnergy);
179 inTree->SetBranchAddress("cluster_trueGammaDirX", &_trueGammaDirX);
180 inTree->SetBranchAddress("cluster_trueGammaDirY", &_trueGammaDirY);
181 inTree->SetBranchAddress("cluster_trueGammaDirZ", &_trueGammaDirZ);
182 inTree->SetBranchAddress("cluster_trueGammaPosX", &_trueGammaPosX);
183 inTree->SetBranchAddress("cluster_trueGammaPosY", &_trueGammaPosY);
184 inTree->SetBranchAddress("cluster_trueGammaPosZ", &_trueGammaPosZ);
185 inTree->SetBranchAddress("cluster_trueNeutralinoEnergy", &_trueNeutralinoEnergy);
186 inTree->SetBranchAddress("cluster_trueNeutralinoDirX", &_trueNeutralinoDirX);
187 inTree->SetBranchAddress("cluster_trueNeutralinoDirY", &_trueNeutralinoDirY);
188 inTree->SetBranchAddress("cluster_trueNeutralinoDirZ", &_trueNeutralinoDirZ);
189 inTree->SetBranchAddress("cluster_truePhiDir", &_truePhiDir);
190 inTree->SetBranchAddress("cluster_trueThetaDir", &_trueThetaDir);
191
192 int nEntries = inTree->GetEntriesFast();
193 cout<<" found input tree with "<<nEntries<<" entries"<<endl;
194
195 //weight events according to expected crosssection
196 // polWeight[fileCounter] = 1;
197 weight[fileCounter] = nExpected[fileCounter]*polWeight[fileCounter]/nEntries;
198 cout<<" and weight "<<weight[fileCounter]<<endl;
199
200 /******** apply selection criteria chosen above **********/
201 int selectedEvents = 0;
202 int notEcluster = 0;
203 int notCosTheta = 0;
204 int notNhits = 0;
205 int notEtmiss = 0;
206 int notEsum = 0;
207 int notNgamma = 0;
208 int isGood[nEntries];
209 int isCosTheta[nEntries];
210
211 for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
212 inTree->GetEntry(entryCounter);
213 isGood[entryCounter] = 0;
214 isCosTheta[entryCounter] = 0;
215 if ( _esumECAL > _EminEcal && _esumECAL < _EmaxEcal) {
216 if ( _etMiss > _minEt_miss) {
217 if ( _nHitsECAL > _minNhitsEcal && _nHitsECAL < _maxNhitsEcal) {
218 if ( _nCluster >= _minNgamma ) {
219 for (int clusterCounter=0; clusterCounter<_nCluster; clusterCounter++) {
220 if ( (_cosThetaPos[clusterCounter] < _maxCosTheta) && (_cosThetaPos[clusterCounter] > -1*_maxCosTheta) ) {
221 isCosTheta[entryCounter]++;
222 if ( _energy[clusterCounter] > _minEcluster ) {
223 isGood[entryCounter]++;
224 }
225 else notEcluster++;
226 }
227 else notCosTheta++;
228 }//all cluster
229 }
230 else notNgamma++;
231 }
232 else notNhits++;
233 }
234 else notEtmiss++;
235 }
236 else notEsum++;
237
238 }//all entries
239 cout<<" looped over input tree "<<endl;
240
241 /******** fill output tree **********/
242 for (int entryCounter = 0; entryCounter < nEntries; ++entryCounter) {
243 inTree->GetEntry(entryCounter);
244
245 if ( isGood[entryCounter] >= _minNgamma && isGood[entryCounter] > 0 && isCosTheta[entryCounter] == _nCluster ) {
246
247 // cout<<isCosTheta[entryCounter]<<" of "<<_nCluster<<" cluster with |cos(theta)| > "<<_maxCosTheta<<endl;
248
249 _clusterFill.iEvt = _iEvt;
250 _clusterFill.esumECAL = _esumECAL;
251 _clusterFill.esumHCAL = _esumHCAL;
252 _clusterFill.esumCalos = _esumCalos;
253 _clusterFill.nHitsECAL = _nHitsECAL;
254 _clusterFill.nHitsHCAL = _nHitsHCAL;
255 _clusterFill.nHitsCalos = _nHitsCalos;
256 _clusterFill.etMiss = _etMiss;
257 _clusterFill.nCluster = _nCluster;
258 _clusterFill.nGamma = _nGamma;
259 _clusterFill.trueNcluster = _trueNcluster;
260
261 for (int clusterCounter=0; clusterCounter<_nCluster; ++clusterCounter) {
262 _clusterFill.energy[clusterCounter] = _energy[clusterCounter];
263 _clusterFill.nHitsPerCluster[clusterCounter] = _nHitsPerCluster[clusterCounter];
264 _clusterFill.clusterPos[clusterCounter][0] = _clusterPos[clusterCounter][0];
265 _clusterFill.clusterPos[clusterCounter][1] = _clusterPos[clusterCounter][1];
266 _clusterFill.clusterPos[clusterCounter][2] = _clusterPos[clusterCounter][2];
267 _clusterFill.clusterAxis[clusterCounter][0] = _clusterAxis[clusterCounter][0];
268 _clusterFill.clusterAxis[clusterCounter][1] = _clusterAxis[clusterCounter][1];
269 _clusterFill.clusterAxis[clusterCounter][2] = _clusterAxis[clusterCounter][2];
270 _clusterFill.thetaDir[clusterCounter] = _thetaDir[clusterCounter];
271 _clusterFill.phiDir[clusterCounter] = _phiDir[clusterCounter];
272 _clusterFill.cosThetaPos[clusterCounter] = _cosThetaPos[clusterCounter];
273 _clusterFill.phiPos[clusterCounter] = _phiPos[clusterCounter];
274 _clusterFill.dca[clusterCounter] = _dca[clusterCounter];
275 _clusterFill.trueDist[clusterCounter] = _trueDist[clusterCounter];
276 _clusterFill.trueGammaEnergy[clusterCounter] = _trueGammaEnergy[clusterCounter];
277 _clusterFill.trueGammaDirX[clusterCounter] = _trueGammaDirX[clusterCounter];
278 _clusterFill.trueGammaDirY[clusterCounter] = _trueGammaDirY[clusterCounter];
279 _clusterFill.trueGammaDirZ[clusterCounter] = _trueGammaDirZ[clusterCounter];
280 _clusterFill.trueGammaPosX[clusterCounter] = _trueGammaPosX[clusterCounter];
281 _clusterFill.trueGammaPosY[clusterCounter] = _trueGammaPosY[clusterCounter];
282 _clusterFill.trueGammaPosZ[clusterCounter] = _trueGammaPosZ[clusterCounter];
283 _clusterFill.trueNeutralinoEnergy[clusterCounter] = _trueNeutralinoEnergy[clusterCounter];
284 _clusterFill.trueNeutralinoDirX[clusterCounter] = _trueNeutralinoDirX[clusterCounter];
285 _clusterFill.trueNeutralinoDirY[clusterCounter] = _trueNeutralinoDirY[clusterCounter];
286 _clusterFill.trueNeutralinoDirZ[clusterCounter] = _trueNeutralinoDirZ[clusterCounter];
287 _clusterFill.trueThetaDir[clusterCounter] = _trueThetaDir[clusterCounter];
288 _clusterFill.truePhiDir[clusterCounter] = _truePhiDir[clusterCounter];
289 }//all cluster
290
291 outTree->SetWeight(weight[fileCounter]);
292 outTree->Fill();
293 selectedEvents++;
294 }//is good
295 }//all entries
296
297 cout<<"writing output tree with "<<selectedEvents<<" selected events"<<endl;
298 cout<<"\t found "<<notEcluster<<" cluster with less then "<<_minEcluster<<" GeV per cluster"<<endl;
299 cout<<"\t found "<<notCosTheta<<" cluster with cos(theta) more then "<<_maxCosTheta<<endl;
300 cout<<"\t found "<<notNhits<<" events with less then "<<_minNhitsEcal<<" hits in ECAL"<<endl;
301 cout<<"\t found "<<notEtmiss<<" events with less then "<<_minEt_miss<<" GeV transverse energy"<<endl;
302 cout<<"\t found "<<notEsum<<" events with less then "<<_EminEcal<<" GeV in ECAL"<<endl;
303 cout<<"\t found "<<notNgamma<<" events with less then "<<_minNgamma<<" photons"<<endl;
304
305 }//file is open
306 inFile->Close();
307 }//all files
308
309 outFile->Write();
310
311 }//eventSelection
312
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.