00001
00002
00003
00004
00005
00006
00008 #include "TFile.h"
00009 #include "TRandom.h"
00010 #include "TString.h"
00011 #include "NueAna/NueHandScan.h"
00012 #include "NueAna/NuePID.h"
00013 #include "MessageService/MsgService.h"
00014 #include "MinosObjectMap/MomNavigator.h"
00015 #include "JobControl/JobCModuleRegistry.h"
00016 #include "JobControl/JobCEnv.h"
00017 #include "MCNtuple/NtpMCRecord.h"
00018 #include "MCNtuple/NtpMCTruth.h"
00019 #include "TruthHelperNtuple/NtpTHRecord.h"
00020 #include "StandardNtuple/NtpStRecord.h"
00021 #include "CandNtupleSR/NtpSRRecord.h"
00022 #include "CandNtupleSR/NtpSREvent.h"
00023 #include "CandNtupleSR/NtpSRTrack.h"
00024 #include "CandNtupleSR/NtpSRShower.h"
00025 #include "NueAna/SntpHelpers.h"
00026 #include <iostream>
00027 #include <fstream>
00028
00029 JOBMODULE(NueHandScan, "NueHandScan",
00030 "Makes Random Hand Scan file from ntuples");
00031 CVSID("$Id:");
00032
00033
00034 NueHandScan::NueHandScan() :
00035 allEntries(0),randomSeed(0),fracPassed(0.8),
00036 preScaleFactorNue(1),preScaleFactorNuMu(1),preScaleFactorNuTau(1),
00037 preScaleFactorBNue(1),preScaleFactorNC(1),
00038 nNueCC(0),nNueNC(0),nNuMuCC(0),nNuMuNC(0),
00039 nNuTauCC(0),nNuTauNC(0),nBeamNueCC(0),nBeamNueNC(0),nPassNueCC(0),
00040 nPassNueNC(0),nPassNuMuCC(0),nPassNuMuNC(0),nPassNuTauCC(0),
00041 nPassNuTauNC(0),nPassBeamNueCC(0),nPassBeamNueNC(0)
00042 {
00043 firstPass = true;
00044 }
00045
00046
00047 NueHandScan::~NueHandScan()
00048 {}
00049
00050
00051 JobCResult NueHandScan::Ana(const MomNavigator* mom)
00052 {
00053
00054 std::string curFileName = this->GetCurrentFile();
00055 if(firstPass) {
00056 ntupleFileNames.push_back(curFileName);
00057 MSG("NueHandScan",Msg::kInfo) << curFileName << endl;
00058 gRandom->SetSeed(randomSeed);
00059 randomSeed = gRandom->GetSeed();
00060 firstPass = false;
00061 }
00062
00063 std::vector<std::string>::iterator iter = ntupleFileNames.begin();
00064 std::vector<std::string>::iterator end = ntupleFileNames.end();
00065 Bool_t fileFlag = false;
00066 while(iter!=end){
00067 if(strcmp(iter->c_str(),curFileName.c_str())==0) fileFlag = true;
00068 iter++;
00069 }
00070 if(!fileFlag) {
00071 ntupleFileNames.push_back(curFileName);
00072 MSG("NueHandScan",Msg::kInfo) << curFileName << endl;
00073 }
00074
00075 TObject *obj=0;
00076 TIter objiter = mom->FragmentIter();
00077 NtpSRRecord *sr=0;
00078 NtpStRecord *st=0;
00079 NtpMCRecord *mc=0;
00080 NtpTHRecord *th=0;
00081 vector<NuePID *> vpid;
00082 while((obj=objiter.Next())){
00083 const char *cn=obj->ClassName();
00084 if(strcmp(cn,"NtpSRRecord")==0){
00085 sr = dynamic_cast<NtpSRRecord *>(obj);
00086 }
00087 else if(strcmp(cn,"NtpStRecord")==0){
00088 st = dynamic_cast<NtpStRecord *>(obj);
00089 }
00090 else if(strcmp(cn,"NtpMCRecord")==0){
00091 mc = dynamic_cast<NtpMCRecord *>(obj);
00092 }
00093 else if(strcmp(cn,"NtpTHRecord")==0){
00094 th = dynamic_cast<NtpTHRecord *>(obj);
00095 }
00096 else if(strcmp(cn,"NuePID")==0){
00097 NuePID *npid = dynamic_cast<NuePID *>(obj);
00098 vpid.push_back(npid);
00099 }
00100 else{
00101 continue;
00102 }
00103 }
00104 if(!sr&&!st){
00105 MSG("NueHandScan",Msg::kError)<<"No NtpSR(t)Record in Mom"<<endl;
00106 return JobCResult::kFailed;
00107 }
00108 if(!st&&!mc){
00109 MSG("NueHandScan",Msg::kError)<<"No NtpMCRecord in Mom"<<endl;
00110 return JobCResult::kFailed;
00111 }
00112 if(!st&&!th){
00113 MSG("NueHandScan",Msg::kError)<<"No NtpTHRecord in Mom"<<endl;
00114 return JobCResult::kFailed;
00115 }
00116
00117
00118
00119
00120 Int_t nEvent = 0;
00121 if(st) nEvent = st->evthdr.nevent;
00122 else if(sr) nEvent = sr->evthdr.nevent;
00123
00124 for(int i=0;i<nEvent;i++){
00125
00126
00127 Int_t pidevno = -1;
00128 if(vpid.size()>0) {
00129 for(unsigned int j=0;j<vpid.size();j++){
00130 if(vpid[j]->GetHeader().GetEventNo()==i){
00131 pidevno = j;
00132 }
00133 }
00134 }
00135
00136 NtpSREvent *event = 0;
00137 Int_t mcindex = -1;
00138 NtpMCTruth *mcth = 0;
00139 Int_t detector = 0;
00140 if(st) {
00141 if(st->GetHeader().GetVldContext().
00142 GetDetector()==Detector::kNear) detector = 1;
00143 else if(st->GetHeader().GetVldContext().
00144 GetDetector()==Detector::kFar) detector = 2;
00145 event = SntpHelpers::GetEvent(i,st);
00146 mcindex = SntpHelpers::GetEvent2MCIndex(i,st);
00147 mcth = SntpHelpers::GetMCTruth(mcindex,st);
00148 runMap[allEntries] = st->GetHeader().GetRun();
00149 subrunMap[allEntries] = st->GetHeader().GetSubRun();
00150 snarlMap[allEntries] = st->GetHeader().GetSnarl();
00151 eventMap[allEntries] = i;
00152 }
00153 else if(sr){
00154 if(sr->GetHeader().GetVldContext().
00155 GetDetector()==Detector::kNear) detector = 1;
00156 else if(sr->GetHeader().GetVldContext().
00157 GetDetector()==Detector::kFar) detector = 2;
00158 event = SntpHelpers::GetEvent(i,sr);
00159 mcindex = SntpHelpers::GetEvent2MCIndex(i,th);
00160 mcth = SntpHelpers::GetMCTruth(mcindex,mc);
00161 runMap[allEntries] = sr->GetHeader().GetRun();
00162 subrunMap[allEntries] = sr->GetHeader().GetSubRun();
00163 snarlMap[allEntries] = sr->GetHeader().GetSnarl();
00164 eventMap[allEntries] = i;
00165 }
00166
00167 Double_t preScale = 1;
00168
00169 if(mcth->iaction == 1){
00170 if(TMath::Abs(mcth->inu)==14) { nNuMuCC++; preScale = preScaleFactorNuMu; }
00171 else if(TMath::Abs(mcth->inu)==16) { nNuTauCC++; preScale = preScaleFactorNuTau; }
00172 else if(TMath::Abs(mcth->inu)==12) {
00173 if(TMath::Abs(mcth->inunoosc)==12) { nBeamNueCC++; preScale = preScaleFactorBNue; }
00174 else { nNueCC++; preScale = preScaleFactorNue; }
00175 }
00176 }
00177 else {
00178 preScale = preScaleFactorNC;
00179 if(TMath::Abs(mcth->inu)==14) nNuMuNC++;
00180 else if(TMath::Abs(mcth->inu)==16) nNuTauNC++;
00181 else if(TMath::Abs(mcth->inu)==12) {
00182 if(TMath::Abs(mcth->inunoosc)==12) nBeamNueNC++;
00183 else nNueNC++;
00184 }
00185 }
00186
00187 if(gRandom->Uniform()<=preScale &&
00188 this->PassCuts(event,detector)) {
00189 passMap[allEntries] = 1;
00190 if(mcth->iaction == 1){
00191 if(TMath::Abs(mcth->inu)==14) nPassNuMuCC++;
00192 else if(TMath::Abs(mcth->inu)==16) nPassNuTauCC++;
00193 else if(TMath::Abs(mcth->inu)==12) {
00194 if(TMath::Abs(mcth->inunoosc)==12) nPassBeamNueCC++;
00195 else nPassNueCC++;
00196 }
00197 }
00198 else {
00199 if(TMath::Abs(mcth->inu)==14) nPassNuMuNC++;
00200 else if(TMath::Abs(mcth->inu)==16) nPassNuTauNC++;
00201 else if(TMath::Abs(mcth->inu)==12) {
00202 if(TMath::Abs(mcth->inunoosc)==12) nPassBeamNueNC++;
00203 else nPassNueNC++;
00204 }
00205 }
00206 }
00207 else passMap[allEntries] = 0;
00208 allEntries+=1;
00209 }
00210 return JobCResult::kPassed;
00211 }
00212
00213
00214 Bool_t NueHandScan::PassCuts(NtpSREvent *event,Int_t detector)
00215 {
00216 if(!event) return false;
00217 if(event->plane.n>16) return false;
00218 if(detector==1){
00219 if(event->vtx.z<1||event->vtx.z>5) return false;
00220 if(TMath::Sqrt(TMath::Power(event->vtx.x-1.4885,2) +
00221 TMath::Power(event->vtx.y-0.1397,2)) > 1) return false;
00222 }
00223 else if(detector==2){
00224 if(event->vtx.z<1||event->vtx.z>29) return false;
00225 if(event->vtx.z>14&&event->vtx.z<17) return false;
00226 if(TMath::Sqrt(TMath::Power(event->vtx.x,2) +
00227 TMath::Power(event->vtx.y,2)) > 3.87) return false;
00228 }
00229 return true;
00230 }
00231
00232
00233 void NueHandScan::EndJob()
00234 {
00235
00236 char filename[256];
00237 sprintf(filename,"randomEventFile_%s.txt",fileTag.c_str());
00238 ofstream outFile(filename);
00239
00240 std::vector<Int_t> alreadyGot(allEntries,0);
00241 for(int i=0;i<allEntries;i++) alreadyGot[i] = 0;
00242
00243 Int_t counter = 0;
00244 Int_t targetNumber = Int_t(fracPassed*(nPassNueCC + nPassNueNC +
00245 nPassNuMuCC + nPassNuMuNC +
00246 nPassNuTauCC + nPassNuTauNC +
00247 nPassBeamNueCC + nPassBeamNueNC));
00248 while(counter<targetNumber){
00249 Int_t ev_no = Int_t(Double_t(allEntries)*gRandom->Uniform());
00250 if(passMap[ev_no]==1&&alreadyGot[ev_no]==0) {
00251 outFile << runMap[ev_no] << " "
00252 << subrunMap[ev_no] << " "
00253 << snarlMap[ev_no] << " "
00254 << eventMap[ev_no] << " 0 0" << endl;
00255 alreadyGot[ev_no] = 1;
00256 counter++;
00257 }
00258 }
00259 outFile.close();
00260
00261 sprintf(filename,"EfficiencySummary_%s.txt",fileTag.c_str());
00262 ofstream outSum(filename);
00263 outSum << "The following ntuple files are required" << endl;
00264 outSum << "to use the associated random event file:" << endl;
00265 std::vector<std::string>::iterator iter = ntupleFileNames.begin();
00266 std::vector<std::string>::iterator end = ntupleFileNames.end();
00267 while(iter!=end) {
00268 TString ts(iter->c_str());
00269 char c[2] = "/";
00270 ts.Remove(0,ts.Last(c[0])+1);
00271 outSum << ts.Data() << endl;
00272 iter++;
00273 }
00274 outSum << "------------------------------------------------------------------" << endl;
00275 outSum << "Total number of reconstructed events in files = " << allEntries << endl;
00276 outSum << "Total number of events passing cuts = "
00277 << (nPassNueCC + nPassNueNC + nPassNuMuCC + nPassNuMuNC +
00278 nPassNuTauCC + nPassNuTauNC + nPassBeamNueCC + nPassBeamNueNC) << endl;
00279 outSum << "Fraction of events passing cuts in random file = " << fracPassed << endl;
00280 outSum << "RandomSeed = " << randomSeed << endl;
00281 outSum << "Pre-Scale Factors (Nue,NuMu,NuTau,BNue,NC): "
00282 << preScaleFactorNue << ", " << preScaleFactorNuMu << ", "
00283 << preScaleFactorNuTau << ", " << preScaleFactorBNue << ", "
00284 << preScaleFactorNC << endl;
00285 outSum << "(pre-scale factors are applied in addition to cuts in table below)" << endl;
00286 outSum << "------------------------------------------------------------------" << endl;
00287 outSum << " \t" << "Total\tPassed\tEff.(%)\tError(%)" << endl;
00288 outSum << "NueCC: \t" << nNueCC << "\t" << nPassNueCC << "\t";
00289 if(nNueCC>0) outSum << 100.*nPassNueCC/nNueCC
00290 << "\t" << 100.*sqrt(nPassNueCC)/nNueCC << endl;
00291 else outSum << "0\t0" << endl;
00292 outSum << "NuMuCC: \t" << nNuMuCC << "\t" << nPassNuMuCC << "\t";
00293 if(nNuMuCC>0) outSum << 100.*nPassNuMuCC/nNuMuCC
00294 << "\t" << 100.*sqrt(nPassNuMuCC)/nNuMuCC << endl;
00295 else outSum << "0\t0" << endl;
00296 outSum << "NuTauCC: \t" << nNuTauCC << "\t" << nPassNuTauCC << "\t";
00297 if(nNuTauCC>0) outSum << 100.*nPassNuTauCC/nNuTauCC
00298 << "\t" << 100.*sqrt(nPassNuTauCC)/nNuTauCC << endl;
00299 else outSum << "0\t0" << endl;
00300 outSum << "BeamNueCC:\t" << nBeamNueCC << "\t" << nPassBeamNueCC << "\t";
00301 if(nBeamNueCC>0) outSum << 100.*nPassBeamNueCC/nBeamNueCC
00302 <<"\t"<<100.*sqrt(nPassBeamNueCC)/nBeamNueCC << endl;
00303 else outSum << "0\t0" << endl;
00304 outSum << "NueNC: \t" << nNueNC << "\t" << nPassNueNC << "\t";
00305 if(nNueNC>0) outSum << 100.*nPassNueNC/nNueNC
00306 << "\t" << 100.*sqrt(nPassNueNC)/nNueNC << endl;
00307 else outSum << "0\t0" << endl;
00308 outSum << "NuMuNC: \t" << nNuMuNC << "\t" << nPassNuMuNC << "\t";
00309 if(nNuMuNC>0) outSum << 100.*nPassNuMuNC/nNuMuNC
00310 << "\t" << 100.*sqrt(nPassNuMuNC)/nNuMuNC << endl;
00311 else outSum << "0\t0" << endl;
00312 outSum << "NuTauNC: \t" << nNuTauNC << "\t" << nPassNuTauNC << "\t";
00313 if(nNuTauNC>0) outSum << 100.*nPassNuTauNC/nNuTauNC
00314 << "\t" << 100.*sqrt(nPassNuTauNC)/nNuTauNC << endl;
00315 else outSum << "0\t0" << endl;
00316 outSum << "BeamNueNC:\t" << nBeamNueNC << "\t" << nPassBeamNueNC << "\t";
00317 if(nBeamNueNC>0) outSum << 100.*nPassBeamNueNC/nBeamNueNC
00318 << "\t" << 100.*sqrt(nPassBeamNueNC)/nBeamNueNC<<endl;
00319 else outSum << "0\t0" << endl;
00320 outSum << "------------------------------------------------------------------" << endl;
00321 outSum.close();
00322 }
00323
00324
00325 const Registry& NueHandScan::DefaultConfig() const
00326 {
00327 static Registry r;
00328 std::string name = this->GetName();
00329 name += ".config.default";
00330 r.SetName(name.c_str());
00331 r.UnLockValues();
00332 r.Set("RandomSeed",0);
00333 r.Set("FracPassed",0.8);
00334 r.Set("FileTag","test");
00335 r.Set("PreScaleFactorNue",1);
00336 r.Set("PreScaleFactorNuMu",1);
00337 r.Set("PreScaleFactorNuTau",1);
00338 r.Set("PreScaleFactorBNue",1);
00339 r.Set("PreScaleFactorNC",1);
00340 r.LockValues();
00341 return r;
00342 }
00343
00344
00345 void NueHandScan::Config(const Registry& r)
00346 {
00347 Int_t tmpi = 0;
00348 Double_t tmpd = 0;
00349 const char *tmpc = 0;
00350 if(r.Get("RandomSeed",tmpi)) { randomSeed = tmpi; }
00351 if(r.Get("FracPassed",tmpd)) { fracPassed = tmpd; }
00352 if(r.Get("FileTag",tmpc)) { fileTag = tmpc; }
00353 if(r.Get("PreScaleFactorNue",tmpd)) { preScaleFactorNue = tmpd; }
00354 if(r.Get("PreScaleFactorNuMu",tmpd)) { preScaleFactorNuMu = tmpd; }
00355 if(r.Get("PreScaleFactorNuTau",tmpd)) { preScaleFactorNuTau = tmpd; }
00356 if(r.Get("PreScaleFactorBNue",tmpd)) { preScaleFactorBNue = tmpd; }
00357 if(r.Get("PreScaleFactorNC",tmpd)) { preScaleFactorNC = tmpd; }
00358 }
00359
