00001 //Contrary to its name, this code does not make histograms. Rather, it skims the trees looking for electrons 00002 //and writes them to a slimmer tree 00003 #include <iostream> 00004 #include <fstream> 00005 #include <set> 00006 #include <pair> 00007 #include <map> 00008 00009 using namespace std; 00010 00011 void electron_ntuple_maker(const char* file_list="",const char* skimfile="electronskimfile.root") 00012 { 00013 gROOT->Macro("LoadLogger.C"); 00014 gROOT->Macro("loadMuDst.C"); 00015 gSystem->Load("StTpcDb"); 00016 gSystem->Load("StDaqLib"); 00017 gSystem->Load("StDetectorDbMaker"); 00018 gSystem->Load("St_db_Maker"); 00019 gSystem->Load("StDbUtilities"); 00020 gSystem->Load("StEmcRawMaker"); 00021 gSystem->Load("StMcEvent"); 00022 gSystem->Load("StMcEventMaker");//*** 00023 gSystem->Load("StEmcSimulatorMaker");//*** 00024 gSystem->Load("StEmcADCtoEMaker"); 00025 gSystem->Load("StEpcMaker"); 00026 gSystem->Load("StDbBroker"); 00027 gSystem->Load("StEmcUtil"); 00028 gSystem->Load("StAssociationMaker"); 00029 gSystem->Load("StEmcTriggerMaker"); 00030 gSystem->Load("StTriggerUtilities"); 00031 gSystem->Load("StEmcOfflineCalibrationMaker"); 00032 00033 //chain all input files together 00034 char file[300]; 00035 TChain* calib_tree = new TChain("calibTree"); 00036 ifstream filelist(file_list); 00037 while(1){ 00038 filelist >> file; 00039 if(!filelist.good()) break; 00040 cout<<file<<endl; 00041 calib_tree->Add(file); 00042 } 00043 00044 char* dbtime = "2008-02-28 00:00:00"; 00045 StEmcOfflineCalibrationElectronAnalyzer* ana = new StEmcOfflineCalibrationElectronAnalyzer; 00046 00047 ana->HTtrigs.push_back(220500); 00048 ana->HTtrigs.push_back(220510); 00049 ana->HTtrigs.push_back(220520); 00050 00051 /*2006 00052 ana->HTtrigs.push_back(117211); 00053 ana->HTtrigs.push_back(117212); 00054 ana->HTtrigs.push_back(117611); 00055 ana->HTtrigs.push_back(117821); 00056 ana->HTtrigs.push_back(127212); 00057 ana->HTtrigs.push_back(127213); 00058 ana->HTtrigs.push_back(127611); 00059 ana->HTtrigs.push_back(127821); 00060 ana->HTtrigs.push_back(137213); 00061 ana->HTtrigs.push_back(137611); 00062 ana->HTtrigs.push_back(137821); 00063 */ 00064 ana->analyze(calib_tree,skimfile,dbtime); 00065 00066 /* 00067 StEmcGeom* emcgeom = StEmcGeom::instance("bemc"); 00068 StBemcTablesWriter* bemctables = new StBemcTablesWriter(); 00069 bemctables->loadTables("2006-05-31 00:00:00","ofl"); 00070 cout<<"input filelist: "<<file_list<<endl; 00071 cout<<"skimmed tree: "<<skimfile<<endl; 00072 00073 StEmcOfflineCalibrationEvent* myEvent = new StEmcOfflineCalibrationEvent(); 00074 calib_tree->SetBranchAddress("event_branch",&myEvent); 00075 StEmcOfflineCalibrationTrack* track = new StEmcOfflineCalibrationTrack(); 00076 StEmcOfflineCalibrationTrack* dumtrack = new StEmcOfflineCalibrationTrack(); 00077 00078 TFile* skim_file = new TFile(skimfile,"RECREATE"); 00079 TNtuple *ntuple = new TNtuple("nt","electron ntuple","p:teta:tphi:dedx:np:id:idx:eta:energy:r:vz:ht:en3x3:enHN:nTracks"); 00080 00081 //keep track of all hit towers and exclude any with >1 track/tower 00082 set<int> track_towers; 00083 set<int> excluded_towers; 00084 00085 unsigned int nAccept = 0; 00086 unsigned int nGoodEvents = 0; 00087 unsigned int nentries = calib_tree->GetEntries(); 00088 for(unsigned int i=0; i<nentries; i++){ 00089 if(i%100000 == 0) cout<<"processing "<<i<<" of "<<nentries<<endl; 00090 //cout<<i<<endl; 00091 calib_tree->GetEntry(i); 00092 00093 //event level cuts 00094 float vz = myEvent->vz[0]; 00095 00096 int ht = 0; 00097 for(map<int, int>::iterator iter = (map<int, int>::iterator)myEvent->triggerResult.begin(); iter != (map<int, int>::iterator)myEvent->triggerResult.end(); ++iter){ 00098 if((*iter).first < 200000 && (*iter).second == 1)ht = 1; 00099 } 00100 00101 //do a quick loop over tracks to get excluded towers 00102 track_towers.clear(); 00103 excluded_towers.clear(); 00104 00105 for(int j=0; j<myEvent->tracks->GetEntries(); j++){ 00106 track = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j); 00107 int id = track->tower_id[0]; 00108 00109 if(track_towers.find(id) != track_towers.end()){ 00110 excluded_towers.insert(id); 00111 } 00112 else{ 00113 track_towers.insert(id); 00114 } 00115 } 00116 00117 //now we loop again and look for electrons / hadrons 00118 for(int j=0; j<myEvent->tracks->GetEntries(); j++){ 00119 track = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j); 00120 00121 if(j==0) nGoodEvents++; 00122 00123 double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi); 00124 //if(dR > 0.03) continue; 00125 //cout<<track->nSigmaElectron<<" "<<track->p<<" "<<track->nHits<<" "<<track->tower_id[0]<<" "<<track->tower_id_exit<<endl; 00126 float squarefid = 0.03/TMath::Sqrt(2.0); 00127 //if(TMath::Abs(track->deta) > squarefid || TMath::Abs(track->dphi) > squarefid)continue; 00128 00129 float p = track->p; 00130 float teta = track->track.pseudoRapidity(); 00131 float tphi = track->track.phi(); 00132 float np = track->nHits; 00133 float dedx = track->dEdx; 00134 float energy = (track->tower_adc[0]-track->tower_pedestal[0])*bemctables->calib(1,track->tower_id[0]); 00135 if(energy < 0)energy = 0; 00136 if(track->tower_status[0] != 1) continue; 00137 float id = track->tower_id[0]; 00138 float idx = track->tower_id_exit; 00139 float eta; 00140 emcgeom->getEta(id,eta); 00141 if(excluded_towers.find(track->tower_id[0]) != excluded_towers.end()) continue; 00142 00143 //looking for tracks at surrounding towers 00144 //cout<<"passed basic cuts"<<endl; 00145 00146 float nnt = 0; 00147 float en3x3 = energy; 00148 float enHN = 0; 00149 if(track->nSigmaElectron > -5.){ 00150 //create start using cluster 00151 00152 for (int k = 1; k < 9; k++)//loop over surrounding towers 00153 { 00154 float tenergy = (track->tower_adc[k]-track->tower_pedestal[k])*bemctables->calib(1,track->tower_id[k]); 00155 if(tenergy > 0)en3x3+=tenergy; 00156 if(tenergy > enHN)enHN = tenergy; 00157 if(track_towers.find(track->tower_id[k]) != track_towers.end())//if there's a tower with a track 00158 { 00159 nnt++; 00160 } 00161 } 00162 nt->Fill(p,teta,tphi,dedx,np,id,idx,eta,energy,dR,vz,ht,en3x3,enHN,nnt); 00163 nAccept++; 00164 track->Clear(); 00165 } 00166 } 00167 myEvent->Clear(); 00168 } 00169 00170 00171 cout<<"found "<<nGoodEvents<<" events with at least one good track"<<endl; 00172 cout<<"accepted "<<nAccept<<" electrons"<<endl; 00173 00174 skim_file->cd(); 00175 skim_file->Write(); 00176 skim_file->Close(); 00177 */ 00178 }
1.5.9