WCSim
nuPRISMconvert.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <TH1F.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include <TFile.h>
6 #include <TTree.h>
7 #include <TKey.h>
8 #include <WCSimRootEvent.hh>
9 
10 /* converts position and direction branches in the fiTQun tree to nuPRISM coordinate system
11  * and copies the rest of the branches into a new fiTQun tree
12  */
13 void convertfiTQun(TTree* fiTQun, TTree* fiTQunCv) {
14 
15  // copies fiTQun tree to a new tree without any position and direction branches
16  fiTQun->SetBranchStatus("*", 1);
17  fiTQun->SetBranchStatus("fq*pos*",0);
18  fiTQun->SetBranchStatus("fq*dir*",0);
19  fiTQunCv = fiTQun->CloneTree();
20  fiTQun->CopyAddresses(fiTQunCv, kTRUE);
21 
22  // sets branches back to be processed
23  fiTQun->SetBranchStatus("fq*pos*",1);
24  fiTQun->SetBranchStatus("fq*dir*",1);
25 
26  // initialize position and direction variables and their converted counterparts
27  const Int_t MAX_TEMP = 100;
28  Int_t fqntwnd;
29  Float_t fqtwnd_prftpos[MAX_TEMP][3]; // pre-fitter vertex position
30  Float_t fqtwnd_prftposCv[MAX_TEMP][3];
31  Int_t fqnse;
32  Float_t fq1rpos[MAX_TEMP][7][3]; // 1 ring fit vertex
33  Float_t fq1rposCv[MAX_TEMP][7][3];
34  Float_t fq1rdir[MAX_TEMP][7][3]; // 1 ring fit direction
35  Float_t fq1rdirCv[MAX_TEMP][7][3];
36  Float_t fqpi0pos[2][3]; // Pi0 fit vertex position
37  Float_t fqpi0posCv[2][3];
38  Float_t fqpi0dir1[2][3]; // Pi0 fit direction of the first photon
39  Float_t fqpi0dir1Cv[2][3];
40  Float_t fqpi0dir2[2][3]; // Pi0 fit direction of the second photon
41  Float_t fqpi0dir2Cv[2][3];
42  Float_t fqpi0dirtot[2][3]; // Pi0 fit direction of the Pi0
43  Float_t fqpi0dirtotCv[2][3];
44  Int_t fqnmrfit;
45  Float_t fqmrpos[MAX_TEMP][6][3]; // Multi-Ring fit vertex position of each ring
46  Float_t fqmrposCv[MAX_TEMP][6][3];
47  Float_t fqmrdir[MAX_TEMP][6][3]; // Multi-Ring fit direction of each ring
48  Float_t fqmrdirCv[MAX_TEMP][6][3];
49  Int_t fqmsnfit;
50  Float_t fqmspos[MAX_TEMP][20][3]; // Multi-Segment Muon fit vertex position of each segment
51  Float_t fqmsposCv[MAX_TEMP][20][3];
52  Float_t fqmsdir[MAX_TEMP][20][3]; // Multi-Segment Muon fit direction of each segment
53  Float_t fqmsdirCv[MAX_TEMP][20][3];
54  Int_t fqtestn1r;
55  Float_t fqtest1rpos[MAX_TEMP][3]; // test 1 ring fit vertex
56  Float_t fqtest1rposCv[MAX_TEMP][3];
57  Float_t fqtest1rdir[MAX_TEMP][3]; // test 1 ring fit directio
58  Float_t fqtest1rdirCv[MAX_TEMP][3];
59  Int_t fqtestnpi0;
60  Float_t fqtestpi0pos[MAX_TEMP][3]; // test Pi0 fit vertex position
61  Float_t fqtestpi0posCv[MAX_TEMP][3];
62  Float_t fqtestpi0dir1[MAX_TEMP][3]; // test Pi0 fit direction of the first photon
63  Float_t fqtestpi0dir1Cv[MAX_TEMP][3];
64  Float_t fqtestpi0dir2[MAX_TEMP][3]; // test Pi0 fit direction of the second photon
65  Float_t fqtestpi0dir2Cv[MAX_TEMP][3];
66  Float_t fqtestpi0dirtot[MAX_TEMP][3]; // test Pi0 fit direction of the Pi0
67  Float_t fqtestpi0dirtotCv[MAX_TEMP][3];
68 
69  // initializes branches in fiTQunCv output
70  TBranch *bfqtwnd_prftpos = fiTQunCv->Branch("fqtwnd_prftpos",fqtwnd_prftposCv,"fqtwnd_prftpos[fqntwnd][3]/F");
71  TBranch *bfq1rpos = fiTQunCv->Branch("fq1rpos",fq1rposCv,"fq1rpos[fqnse][7][3]/F");
72  TBranch *bfq1rdir = fiTQunCv->Branch("fq1rdir",fq1rdirCv,"fq1rdir[fqnse][7][3]/F");
73  TBranch *bfqpi0pos = fiTQunCv->Branch("fqpi0pos",fqpi0posCv,"fqpi0pos[2][3]/F");
74  TBranch *bfqpi0dir1 = fiTQunCv->Branch("fqpi0dir1",fqpi0dir1Cv,"fqpi0dir1[2][3]/F");
75  TBranch *bfqpi0dir2 = fiTQunCv->Branch("fqpi0dir2",fqpi0dir2Cv,"fqpi0dir2[2][3]/F");
76  TBranch *bfqpi0dirtot = fiTQunCv->Branch("fqpi0dirtot",fqpi0dirtotCv,"fqpi0dirtot[2][3]/F");
77  TBranch *bfqmrpos = fiTQunCv->Branch("fqmrpos",fqmrposCv,"fqmrpos[fqnmrfit][6][3]/F");
78  TBranch *bfqmrdir = fiTQunCv->Branch("fqmrdir",fqmrdirCv,"fqmrdir[fqnmrfit][6][3]/F");
79  TBranch *bfqmspos = fiTQunCv->Branch("fqmspos",fqmsposCv,"fqmspos[fqmsnfit][20][3]/F");
80  TBranch *bfqmsdir = fiTQunCv->Branch("fqmsdir",fqmsdirCv,"fqmsdir[fqmsnfit][20][3]/F");
81  TBranch *bfqtest1rpos = fiTQunCv->Branch("fqtest1rpos",fqtest1rposCv,"fqtest1rpos[fqtestn1r][3]/F");
82  TBranch *bfqtest1rdir = fiTQunCv->Branch("fqtest1rdir",fqtest1rdirCv,"fqtest1rdir[fqtestn1r][3]/F");
83  TBranch *bfqtestpi0pos = fiTQunCv->Branch("fqtestpi0pos",fqtestpi0posCv,"fqtestpi0pos[fqtestnpi0][3]/F");
84  TBranch *bfqtestpi0dir1 = fiTQunCv->Branch("fqtestpi0dir1",fqtestpi0dir1Cv,"fqtestpi0dir1[fqtestnpi0][3]/F");
85  TBranch *bfqtestpi0dir2 = fiTQunCv->Branch("fqtestpi0dir2",fqtestpi0dir2Cv,"fqtestpi0dir2[fqtestnpi0][3]/F");
86  TBranch *bfqtestpi0dirtot = fiTQunCv->Branch("fqtestpi0dirtot",fqtestpi0dirtotCv,"fqtestpi0dirtot[fqtestnpi0][3]/F");
87  fiTQun->SetBranchAddress("fqntwnd", &fqntwnd);
88  fiTQun->SetBranchAddress("fqtwnd_prftpos", &fqtwnd_prftpos);
89  fiTQun->SetBranchAddress("fqnse", &fqnse);
90  fiTQun->SetBranchAddress("fq1rpos", &fq1rpos);
91  fiTQun->SetBranchAddress("fq1rdir", &fq1rdir);
92  fiTQun->SetBranchAddress("fqpi0pos", &fqpi0pos);
93  fiTQun->SetBranchAddress("fqpi0dir1", &fqpi0dir1);
94  fiTQun->SetBranchAddress("fqpi0dir2", &fqpi0dir2);
95  fiTQun->SetBranchAddress("fqpi0dirtot", &fqpi0dirtot);
96  fiTQun->SetBranchAddress("fqnmrfit", &fqnmrfit);
97  fiTQun->SetBranchAddress("fqmrpos", &fqmrpos);
98  fiTQun->SetBranchAddress("fqmrdir", &fqmrdir);
99  fiTQun->SetBranchAddress("fqmsnfit", &fqmsnfit);
100  fiTQun->SetBranchAddress("fqmspos", &fqmspos);
101  fiTQun->SetBranchAddress("fqmsdir", &fqmsdir);
102  fiTQun->SetBranchAddress("fqtestn1r", &fqtestn1r);
103  fiTQun->SetBranchAddress("fqtest1rpos", &fqtest1rpos);
104  fiTQun->SetBranchAddress("fqtest1rdir", &fqtest1rdir);
105  fiTQun->SetBranchAddress("fqtestnpi0", &fqtestnpi0);
106  fiTQun->SetBranchAddress("fqtestpi0pos", &fqtestpi0pos);
107  fiTQun->SetBranchAddress("fqtestpi0dir1", &fqtestpi0dir1);
108  fiTQun->SetBranchAddress("fqtestpi0dir2", &fqtestpi0dir2);
109  fiTQun->SetBranchAddress("fqtestpi0dirtot", &fqtestpi0dirtot);
110 
111  // gets the total number of entries
112  int fqEntries = fiTQun->GetEntries();
113 
114  double offset = 2151.035; // particle gun y-offset of position in centimeters
115 
116  // for each entry in the fiTQun input
117  // changes the fiTQun output's WCSim coordinates to the global nuPRISM coordinate geometry
118  for (int entry=0;entry<fqEntries;entry++) {
119  fiTQun->GetEntry(entry);
120  // time-window information
121  for (int i=0;i<fqntwnd;i++) {
122  for (int axis=0;axis<3;axis++) {
123  if (axis == 0) { // set the x-axis to the z-axis
124  fqtwnd_prftposCv[i][0] = fqtwnd_prftpos[i][axis];
125  } else if (axis == 1) { // set the y-axis to the x-axis
126  fqtwnd_prftposCv[i][1] = fqtwnd_prftpos[i][axis];
127  } else if (axis == 2) { // set the z-axis to the y-axis
128  fqtwnd_prftposCv[i][2] = fqtwnd_prftpos[i][axis];
129  }
130  }
131  }
132  // 1-ring fits
133  for (int i=0;i<fqnse;i++) {
134  for (int j=0;j<7;j++) {
135  for (int axis=0;axis<3;axis++) {
136  if (axis == 0) {
137  fq1rposCv[i][j][0] = fq1rpos[i][j][axis];
138  fq1rdirCv[i][j][0] = fq1rdir[i][j][axis];
139  } else if (axis == 1) {
140  fq1rposCv[i][j][1] = fq1rpos[i][j][axis];
141  fq1rdirCv[i][j][1] = fq1rdir[i][j][axis];
142  } else if (axis == 2) {
143  fq1rposCv[i][j][2] = fq1rpos[i][j][axis];
144  fq1rdirCv[i][j][2] = fq1rdir[i][j][axis];
145  }
146  }
147  }
148  }
149  // Pi0 fits
150  for (int i=0;i<2;i++) {
151  for (int axis=0;axis<3;axis++) {
152  if (axis == 0) {
153  fqpi0posCv[i][0] = fqpi0pos[i][axis];
154  fqpi0dir1Cv[i][0] = fqpi0dir1[i][axis];
155  fqpi0dir2Cv[i][0] = fqpi0dir2[i][axis];
156  fqpi0dirtotCv[i][0] = fqpi0dirtot[i][axis];
157  } else if (axis == 1) {
158  fqpi0posCv[i][1] = fqpi0pos[i][axis];
159  fqpi0dir1Cv[i][1] = fqpi0dir1[i][axis];
160  fqpi0dir2Cv[i][1] = fqpi0dir2[i][axis];
161  fqpi0dirtotCv[i][1] = fqpi0dirtot[i][axis];
162  } else if (axis == 2) {
163  fqpi0posCv[i][2] = fqpi0pos[i][axis];
164  fqpi0dir1Cv[i][2] = fqpi0dir1[i][axis];
165  fqpi0dir2Cv[i][2] = fqpi0dir2[i][axis];
166  fqpi0dirtotCv[i][2] = fqpi0dirtot[i][axis];
167  }
168  }
169  }
170  // multi-ring fits
171  for (int i=0;i<fqnmrfit;i++) {
172  for (int j=0;j<6;j++) {
173  for (int axis=0;axis<3;axis++) {
174  if (axis == 0) {
175  fqmrposCv[i][j][0] = fqmrpos[i][j][axis];
176  fqmrdirCv[i][j][0] = fqmrdir[i][j][axis];
177  } else if (axis == 1) {
178  fqmrposCv[i][j][1] = fqmrpos[i][j][axis];
179  fqmrdirCv[i][j][1] = fqmrdir[i][j][axis];
180  } else if (axis == 2) {
181  fqmrposCv[i][j][2] = fqmrpos[i][j][axis];
182  fqmrdirCv[i][j][2] = fqmrdir[i][j][axis];
183  }
184  }
185  }
186  }
187  // multi-segment muon fits
188  for (int i=0;i<fqmsnfit;i++) {
189  for (int j=0;j<20;j++) {
190  for (int axis=0;axis<3;axis++) {
191  if (axis == 0) {
192  fqmsposCv[i][j][0] = fqmspos[i][j][axis];
193  fqmsdirCv[i][j][0] = fqmsdir[i][j][axis];
194  } else if (axis == 1) {
195  fqmsposCv[i][j][1] = fqmspos[i][j][axis];
196  fqmsdirCv[i][j][1] = fqmsdir[i][j][axis];
197  } else if (axis == 2) {
198  fqmsposCv[i][j][2] = fqmspos[i][j][axis];
199  fqmsdirCv[i][j][2] = fqmsdir[i][j][axis];
200  }
201  }
202  }
203  }
204  // test 1-ring fits
205  for (int i=0;i<fqtestn1r;i++) {
206  for (int axis=0;axis<3;axis++) {
207  if (axis == 0) {
208  fqtest1rposCv[i][0] = fqtest1rpos[i][axis];
209  fqtest1rdirCv[i][0] = fqtest1rdir[i][axis];
210  } else if (axis == 1) {
211  fqtest1rposCv[i][1] = fqtest1rpos[i][axis];
212  fqtest1rdirCv[i][1] = fqtest1rdir[i][axis];
213  } else if (axis == 2) {
214  fqtest1rposCv[i][2] = fqtest1rpos[i][axis];
215  fqtest1rdirCv[i][2] = fqtest1rdir[i][axis];
216  }
217  }
218  }
219  // test Pi0 fits
220  for (int i=0;i<fqtestnpi0;i++) {
221  for (int axis=0;axis<3;axis++) {
222  if (axis == 0) {
223  fqtestpi0posCv[i][0] = fqtestpi0pos[i][axis];
224  fqtestpi0dir1Cv[i][0] = fqtestpi0dir1[i][axis];
225  fqtestpi0dir2Cv[i][0] = fqtestpi0dir2[i][axis];
226  fqtestpi0dirtotCv[i][0] = fqtestpi0dirtot[i][axis];
227  } else if (axis == 1) {
228  fqtestpi0posCv[i][1] = fqtestpi0pos[i][axis];
229  fqtestpi0dir1Cv[i][1] = fqtestpi0dir1[i][axis];
230  fqtestpi0dir2Cv[i][1] = fqtestpi0dir2[i][axis];
231  fqtestpi0dirtotCv[i][1] = fqtestpi0dirtot[i][axis];
232  } else if (axis == 2) {
233  fqtestpi0posCv[i][2] = fqtestpi0pos[i][axis];
234  fqtestpi0dir1Cv[i][2] = fqtestpi0dir1[i][axis];
235  fqtestpi0dir2Cv[i][2] = fqtestpi0dir2[i][axis];
236  fqtestpi0dirtotCv[i][2] = fqtestpi0dirtot[i][axis];
237  }
238  }
239  }
240  // fill all branches for each entry
241  bfqtwnd_prftpos->Fill();
242  bfq1rpos->Fill();
243  bfq1rdir->Fill();
244  bfqpi0pos->Fill();
245  bfqpi0dir1->Fill();
246  bfqpi0dir2->Fill();
247  bfqpi0dirtot->Fill();
248  bfqmrpos->Fill();
249  bfqmrdir->Fill();
250  bfqmspos->Fill();
251  bfqmsdir->Fill();
252  bfqtest1rpos->Fill();
253  bfqtest1rdir->Fill();
254  bfqtestpi0pos->Fill();
255  bfqtestpi0dir1->Fill();
256  bfqtestpi0dir2->Fill();
257  bfqtestpi0dirtot->Fill();
258  }
259 
260  fiTQunCv->Write("",TObject::kOverwrite);
261 
262 }
263 
264 /* extracts the true information from the wcsimT tree
265  * the position and direction information is converted to nuPRISM geometry
266  * and then all the information are saved in the new wcsimEx tree
267  */
268 void extractTruth(TTree* wcsimT, TTree* wcsimEx) {
269 
270  WCSimRootTrigger * wcsimrootTrigger;
271  WCSimRootEvent * wcsimrootEvent;
272  wcsimrootEvent = new WCSimRootEvent();
273  wcsimT->SetBranchAddress("wcsimrootevent" ,&wcsimrootEvent);
274 
275  // initializes variables
276  const Int_t maxtrack = 200;
277  Int_t ntrack;
278  Int_t parent[maxtrack];
279  Int_t pid[maxtrack];
280  Double_t dir[maxtrack][4];
281  Double_t start[maxtrack][4];
282  Double_t stop[maxtrack][4];
283  Double_t mom[maxtrack];
284  Double_t mass[maxtrack];
285  Double_t energy[maxtrack];
286  Double_t startvol[maxtrack];
287  Double_t stopvol[maxtrack];
288  Double_t time[maxtrack];
289 
290  // sets output Ttree to above variables
291  wcsimEx->Branch("ntrack",&ntrack,"ntrack/I");
292  wcsimEx->Branch("parent",parent,"parent[ntrack]/I");
293  wcsimEx->Branch("pid",pid,"pid[ntrack]/I");
294  wcsimEx->Branch("dir",dir,"dir[ntrack][4]/D");
295  wcsimEx->Branch("start",start,"start[ntrack][4]/D");
296  wcsimEx->Branch("stop",stop,"stop[ntrack][4]/D");
297  wcsimEx->Branch("mom",mom,"mom[ntrack]/D");
298  wcsimEx->Branch("mass",mass,"mass[ntrack]/D");
299  wcsimEx->Branch("energy",energy,"energy[ntrack]/D");
300  wcsimEx->Branch("startvol",startvol,"startvol[ntrack]/D");
301  wcsimEx->Branch("stopvol",stopvol,"stopvol[ntrack]/D");
302  wcsimEx->Branch("time",time,"time[ntrack]/D");
303 
304  // for each event in the wcsimT input
305  int Entries = wcsimT->GetEntries();
306  for (int ev=0; ev<Entries; ev++) {
307  wcsimT->GetEvent(ev);
308  wcsimrootTrigger = wcsimrootEvent->GetTrigger(0);
309  ntrack = wcsimrootTrigger->GetNtrack();
310  if(ntrack>0) {
311  // Loop through elements in the TClonesArray of WCSimTracks
312  // for each particle in the event set the respective variable's information
313  for (int par=0; par<ntrack; par++) {
314  // wcsimroottrack holds all info about particle
315  TObject *element = (wcsimrootTrigger->GetTracks())->At(par);
316  WCSimRootTrack *wcsimroottrack = dynamic_cast<WCSimRootTrack*>(element);
317  for (int axis=0;axis<4;axis++) {
318  if (axis == 0) { // set wcsim x-axis to nuPRISM x-axis
319  start[par][0] = wcsimroottrack->GetStart(axis)/100;
320  stop[par][0] = wcsimroottrack->GetStop(axis)/100;
321  dir[par][0] = wcsimroottrack->GetDir(axis);
322  } else if (axis == 1) { // set wcsim y-axis to nuPRISM y-axis
323  start[par][1] = wcsimroottrack->GetStart(axis)/100;
324  stop[par][1] = wcsimroottrack->GetStop(axis)/100;
325  dir[par][1] = wcsimroottrack->GetDir(axis);
326  } else if (axis == 2) { // set wcsim z-axis to nuPRISM z-axis
327  start[par][2] = wcsimroottrack->GetStart(axis)/100;
328  stop[par][2] = wcsimroottrack->GetStop(axis)/100;
329  dir[par][2] = wcsimroottrack->GetDir(axis);
330  } else if (axis == 3) {
331  start[par][3] = wcsimroottrack->GetStart(axis)/100;
332  stop[par][3] = wcsimroottrack->GetStop(axis)/100;
333  dir[par][3] = wcsimroottrack->GetDir(axis);
334  }
335  }
336  parent[par]=wcsimroottrack->GetParenttype();
337  pid[par]=wcsimroottrack->GetIpnu();
338  mom[par]=wcsimroottrack->GetP();
339  mass[par]=wcsimroottrack->GetM();
340  energy[par]=wcsimroottrack->GetE();
341  startvol[par]=wcsimroottrack->GetStartvol();
342  stopvol[par]=wcsimroottrack->GetStopvol();
343  time[par]=wcsimroottrack->GetTime();
344  }
345  // fill the wcsimEx with all the current event's information in the current loop
346  wcsimEx->Fill();
347  }
348  }
349 
350  wcsimEx->Write("",TObject::kOverwrite);
351 
352 }
353 
354 int convertNuPrism(char* input, char* output, bool dofiTQun=0, bool copy=0){
355 
356  char *filename = input;
357  char *outfilename = output;
358  bool toCopy = copy;
359  bool useFiTQun = dofiTQun;
360 
361  TFile *file = new TFile(filename,"read");
362  if (!file->IsOpen()){
363  std::cout << "Error, could not open input file: " << filename << std::endl;
364  exit(-1);
365  }
366 
367  // creates new trees for output file
368  TFile* outFile = TFile::Open(outfilename, "RECREATE");
369  TTree* fiTQunCv = new TTree("fiTQunCv", "fiTQun Converted Geometry");
370  TTree* wcsimEx = new TTree("wcsimEx","Converted Particle Info");
371 
372  // gets the old trees from the fiTQun file
373  TTree* wcsimT = (TTree*) file->Get("wcsimT");
374  TTree* wcsimGeoT = (TTree*) file->Get("wcsimGeoT");
375  TTree* fiTQun = (TTree*) file->Get("fiTQun"); // gets the old fiTQun information
376  TTree* Settings = (TTree*) file->Get("Settings");
377 
378  // prepares trees to be copied
379  TTree* newEventTree;
380  TTree* newGeomTree;
381  TTree* newSettings;
382  if (toCopy) {
383  newEventTree = wcsimT->CloneTree();
384  newGeomTree = wcsimGeoT->CloneTree();
385  if (Settings) {
386  newSettings = Settings->CloneTree();
387  }
388  }
389 
390  if(useFiTQun) convertfiTQun(fiTQun, fiTQunCv);
391  extractTruth(wcsimT, wcsimEx);
392 
393  // copies the other trees
394  if (toCopy){
395  newEventTree->Write();
396  newGeomTree->Write();
397  if (Settings) {
398  newSettings->Write();
399  }
400  }
401 
402  outFile->Close();
403  return 0;
404 }
405 
void convertfiTQun(TTree *fiTQun, TTree *fiTQunCv)
Float_t GetStop(Int_t i=0) const
int convertNuPrism(char *input, char *output, bool dofiTQun=0, bool copy=0)
Int_t GetIpnu() const
dictionary pid
Definition: MakeKin.py:11
TClonesArray * GetTracks() const
Float_t GetDir(Int_t i=0) const
Int_t GetParenttype() const
Float_t GetStart(Int_t i=0) const
Float_t GetM() const
Class storing trigger information.
Float_t GetE() const
WCSimRootTrigger * GetTrigger(int number)
Int_t GetStartvol() const
string filename
Definition: MakeKin.py:235
Class containing event information.
Int_t GetNtrack() const
Double_t GetTime() const
void extractTruth(TTree *wcsimT, TTree *wcsimEx)
Int_t GetStopvol() const
Class holding true track information.
Float_t GetP() const