StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
electron_ntuple_maker.C
1 //Contrary to its name, this code does not make histograms. Rather, it skims the trees looking for electrons
2 //and writes them to a slimmer tree
3 #include <iostream>
4 #include <fstream>
5 #include <set>
6 #include <pair>
7 #include <map>
8 
9 using namespace std;
10 
11 void electron_ntuple_maker(const char* file_list="",const char* skimfile="electronskimfile.root")
12 {
13  gROOT->Macro("LoadLogger.C");
14  gROOT->Macro("loadMuDst.C");
15  gSystem->Load("StTpcDb");
16  gSystem->Load("StDaqLib");
17  gSystem->Load("StDetectorDbMaker");
18  gSystem->Load("St_db_Maker");
19  gSystem->Load("StDbUtilities");
20  gSystem->Load("StEmcRawMaker");
21  gSystem->Load("StMcEvent");
22  gSystem->Load("StMcEventMaker");//***
23  gSystem->Load("StEmcSimulatorMaker");//***
24  gSystem->Load("StEmcADCtoEMaker");
25  gSystem->Load("StEpcMaker");
26  gSystem->Load("StDbBroker");
27  gSystem->Load("StEmcUtil");
28  gSystem->Load("StAssociationMaker");
29  gSystem->Load("StEmcTriggerMaker");
30  gSystem->Load("StTriggerUtilities");
31  gSystem->Load("StEmcOfflineCalibrationMaker");
32 
33  //chain all input files together
34  char file[300];
35  TChain* calib_tree = new TChain("calibTree");
36  ifstream filelist(file_list);
37  while(1){
38  filelist >> file;
39  if(!filelist.good()) break;
40  cout<<file<<endl;
41  calib_tree->Add(file);
42  }
43 
44  char* dbtime = "2008-02-28 00:00:00";
45  StEmcOfflineCalibrationElectronAnalyzer* ana = new StEmcOfflineCalibrationElectronAnalyzer;
46 
47  ana->HTtrigs.push_back(220500);
48  ana->HTtrigs.push_back(220510);
49  ana->HTtrigs.push_back(220520);
50 
51  /*2006
52  ana->HTtrigs.push_back(117211);
53  ana->HTtrigs.push_back(117212);
54  ana->HTtrigs.push_back(117611);
55  ana->HTtrigs.push_back(117821);
56  ana->HTtrigs.push_back(127212);
57  ana->HTtrigs.push_back(127213);
58  ana->HTtrigs.push_back(127611);
59  ana->HTtrigs.push_back(127821);
60  ana->HTtrigs.push_back(137213);
61  ana->HTtrigs.push_back(137611);
62  ana->HTtrigs.push_back(137821);
63  */
64  ana->analyze(calib_tree,skimfile,dbtime);
65 
66  /*
67  StEmcGeom* emcgeom = StEmcGeom::instance("bemc");
68  StBemcTablesWriter* bemctables = new StBemcTablesWriter();
69  bemctables->loadTables("2006-05-31 00:00:00","ofl");
70  cout<<"input filelist: "<<file_list<<endl;
71  cout<<"skimmed tree: "<<skimfile<<endl;
72 
73  StEmcOfflineCalibrationEvent* myEvent = new StEmcOfflineCalibrationEvent();
74  calib_tree->SetBranchAddress("event_branch",&myEvent);
75  StEmcOfflineCalibrationTrack* track = new StEmcOfflineCalibrationTrack();
76  StEmcOfflineCalibrationTrack* dumtrack = new StEmcOfflineCalibrationTrack();
77 
78  TFile* skim_file = new TFile(skimfile,"RECREATE");
79  TNtuple *ntuple = new TNtuple("nt","electron ntuple","p:teta:tphi:dedx:np:id:idx:eta:energy:r:vz:ht:en3x3:enHN:nTracks");
80 
81  //keep track of all hit towers and exclude any with >1 track/tower
82  set<int> track_towers;
83  set<int> excluded_towers;
84 
85  unsigned int nAccept = 0;
86  unsigned int nGoodEvents = 0;
87  unsigned int nentries = calib_tree->GetEntries();
88  for(unsigned int i=0; i<nentries; i++){
89  if(i%100000 == 0) cout<<"processing "<<i<<" of "<<nentries<<endl;
90  //cout<<i<<endl;
91  calib_tree->GetEntry(i);
92 
93  //event level cuts
94  float vz = myEvent->vz[0];
95 
96  int ht = 0;
97  for(map<int, int>::iterator iter = (map<int, int>::iterator)myEvent->triggerResult.begin(); iter != (map<int, int>::iterator)myEvent->triggerResult.end(); ++iter){
98  if((*iter).first < 200000 && (*iter).second == 1)ht = 1;
99  }
100 
101  //do a quick loop over tracks to get excluded towers
102  track_towers.clear();
103  excluded_towers.clear();
104 
105  for(int j=0; j<myEvent->tracks->GetEntries(); j++){
106  track = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j);
107  int id = track->tower_id[0];
108 
109  if(track_towers.find(id) != track_towers.end()){
110  excluded_towers.insert(id);
111  }
112  else{
113  track_towers.insert(id);
114  }
115  }
116 
117  //now we loop again and look for electrons / hadrons
118  for(int j=0; j<myEvent->tracks->GetEntries(); j++){
119  track = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j);
120 
121  if(j==0) nGoodEvents++;
122 
123  double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);
124  //if(dR > 0.03) continue;
125  //cout<<track->nSigmaElectron<<" "<<track->p<<" "<<track->nHits<<" "<<track->tower_id[0]<<" "<<track->tower_id_exit<<endl;
126  float squarefid = 0.03/TMath::Sqrt(2.0);
127  //if(TMath::Abs(track->deta) > squarefid || TMath::Abs(track->dphi) > squarefid)continue;
128 
129  float p = track->p;
130  float teta = track->track.pseudoRapidity();
131  float tphi = track->track.phi();
132  float np = track->nHits;
133  float dedx = track->dEdx;
134  float energy = (track->tower_adc[0]-track->tower_pedestal[0])*bemctables->calib(1,track->tower_id[0]);
135  if(energy < 0)energy = 0;
136  if(track->tower_status[0] != 1) continue;
137  float id = track->tower_id[0];
138  float idx = track->tower_id_exit;
139  float eta;
140  emcgeom->getEta(id,eta);
141  if(excluded_towers.find(track->tower_id[0]) != excluded_towers.end()) continue;
142 
143  //looking for tracks at surrounding towers
144  //cout<<"passed basic cuts"<<endl;
145 
146  float nnt = 0;
147  float en3x3 = energy;
148  float enHN = 0;
149  if(track->nSigmaElectron > -5.){
150  //create start using cluster
151 
152  for (int k = 1; k < 9; k++)//loop over surrounding towers
153  {
154  float tenergy = (track->tower_adc[k]-track->tower_pedestal[k])*bemctables->calib(1,track->tower_id[k]);
155  if(tenergy > 0)en3x3+=tenergy;
156  if(tenergy > enHN)enHN = tenergy;
157  if(track_towers.find(track->tower_id[k]) != track_towers.end())//if there's a tower with a track
158  {
159  nnt++;
160  }
161  }
162  nt->Fill(p,teta,tphi,dedx,np,id,idx,eta,energy,dR,vz,ht,en3x3,enHN,nnt);
163  nAccept++;
164  track->Clear();
165  }
166  }
167  myEvent->Clear();
168  }
169 
170 
171  cout<<"found "<<nGoodEvents<<" events with at least one good track"<<endl;
172  cout<<"accepted "<<nAccept<<" electrons"<<endl;
173 
174  skim_file->cd();
175  skim_file->Write();
176  skim_file->Close();
177 */
178 }
Generated by   doxygen 1.8.5