00001 #include "NCPOTCounter.h"
00002
00003 #include "NCUtils/NCType.h"
00004 using namespace NCType;
00005
00006 #include "MessageService/MsgService.h"
00007
00008 CVSID("$Id: NCPOTCounter.cxx,v 1.14 2009/02/13 18:21:37 bckhouse Exp $");
00009
00010 #include "TFile.h"
00011 #include "TH2F.h"
00012 #include "TList.h"
00013 #include "TString.h"
00014 #include "TSystemDirectory.h"
00015 #include "TRegexp.h"
00016
00017 #include <cassert>
00018
00019 using std::map;
00020 using std::vector;
00021
00022
00023 NCPOTCounter::Cache NCPOTCounter::fgCache;
00024
00025 NCPOTCounter::NCPOTCounter(TString dataMCPath,
00026 TString mockDataSet,
00027 bool useMCAsData,
00028 bool useMockData,
00029 int fileLimitNearData,
00030 int fileLimitFarData,
00031 int fileLimitNearMC,
00032 int fileLimitFarMC,
00033 vector<BeamType::BeamType_t> beamIndex,
00034 int runToUse,
00035 int mockDataSubRun,
00036 bool scaleNearExposures,
00037 double farExposureScaleVal)
00038 :fDataMCPath(dataMCPath),
00039 fMockDataSet(mockDataSet),
00040 fBeamIndex(beamIndex),
00041 fRunToUse(runToUse),
00042 fMockDataSubRun(mockDataSubRun),
00043 fScaleNearExposures(scaleNearExposures),
00044 fFarExposureScaleVal(farExposureScaleVal)
00045 {
00046 fFileLimits[Detector::kNear][NCType::kData] = fileLimitNearData;
00047 fFileLimits[Detector::kFar] [NCType::kData] = fileLimitFarData;
00048 fFileLimits[Detector::kNear][NCType::kMC] = fileLimitNearMC;
00049 fFileLimits[Detector::kFar] [NCType::kMC] = fileLimitFarMC;
00050
00051 fMockOrMCAsData = useMCAsData || useMockData;
00052 }
00053
00054
00055 void NCPOTCounter::SetPOTValues(map<NCBeam::Info, double>& pot,
00056 map<NCBeam::Info, double>& scale,
00057 Detector::Detector_t det,
00058 NCType::EFileType fileType,
00059 NCType::EDataMC dataMC)
00060 {
00061 POTKey key(3);
00062 key[0] = det;
00063 key[1] = fileType;
00064 key[2] = dataMC;
00065
00066 Cache::iterator it = fgCache.find(key);
00067
00068 if(it != fgCache.end()){
00069 pot = it->second.first;
00070 scale = it->second.second;
00071 return;
00072 }
00073
00074
00075 pot.clear();
00076 scale.clear();
00077
00078 MSG("NCPOTCounter", Msg::kInfo)
00079 << "NCPOTCounter::SetPOTValues:"
00080 << " fDataMCPath = " << fDataMCPath << endl;
00081
00082
00083 vector<vector<double> > POT;
00084 POT.resize(NCType::kMaxRun);
00085 for(int i = 0; i < NCType::kMaxRun; ++i){
00086 POT[i].resize(fBeamIndex.size());
00087 for(unsigned int j = 0; j < fBeamIndex.size(); ++j){
00088 POT[i][j] = 0.;
00089 }
00090 }
00091
00092
00093
00094 for(UInt_t bt = 0; bt < fBeamIndex.size(); ++bt){
00095
00096 const vector<TString> fileList = GetListOfFiles(bt, fileType, det, dataMC);
00097 for(unsigned int nf = 0; nf < fileList.size(); ++nf){
00098 const TString fileName = fileList[nf];
00099
00100 const double potToAdd = GetPOTForSingleFile(fileName, fileType,
00101 det, dataMC, bt);
00102 assert(potToAdd >= 0);
00103
00104 const int runType = NCType::FindRunType(fileName);
00105 POT[runType-1].at(bt) += potToAdd;
00106
00107 MSG("NCPOTCounter", Msg::kInfo) << " Adding " << potToAdd << " POT to "
00108 << FileDesignator(fileType, det, dataMC)
00109 << endl;
00110 }
00111
00112
00113 const TString beamType = BeamType::AsString(fBeamIndex[bt]);
00114 MSG("NCPOTCounter", Msg::kInfo) << "POT INFORMATION("
00115 << (det == Detector::kFar ? "far" : "near")
00116 << "): " << endl
00117 << FileDesignator(fileType, det, dataMC)
00118 << " run I POT[" << beamType << "] = "
00119 << POT[0].at(bt)
00120 << " run II POT[" << beamType << "] = "
00121 << POT[1].at(bt)
00122 << endl;
00123
00124 for(int i = 0; i < NCType::kMaxRun; ++i){
00125 const NCBeam::Info beamInfo(fBeamIndex[bt], NCType::ERunType(i+1));
00126 if(POT[i].at(bt) > 0.){
00127 pot[beamInfo] = POT[i][bt];
00128 scale[beamInfo] = 1.;
00129 }
00130 }
00131 }
00132
00133 DoScaling(pot, scale, det, dataMC);
00134
00135 fgCache[key] = POTVal(pot, scale);
00136 }
00137
00138
00139 double NCPOTCounter::GetPOTValue(Detector::Detector_t det,
00140 NCType::EFileType fileType,
00141 NCType::EDataMC dataMC,
00142 BeamType::BeamType_t beamType,
00143 NCType::ERunType run)
00144 {
00145 map<NCBeam::Info, double> pot;
00146 map<NCBeam::Info, double> scale;
00147 SetPOTValues(pot, scale, det, fileType, dataMC);
00148
00149 const NCBeam::Info beamInfo(beamType, run);
00150
00151 return pot[beamInfo];
00152 }
00153
00154
00155 double NCPOTCounter::GetScaleValue(Detector::Detector_t det,
00156 NCType::EFileType fileType,
00157 NCType::EDataMC dataMC,
00158 BeamType::BeamType_t beamType,
00159 NCType::ERunType run)
00160 {
00161 POTMap pot;
00162 ScaleMap scale;
00163 SetPOTValues(pot, scale, det, fileType, dataMC);
00164
00165 const NCBeam::Info beamInfo(beamType, run);
00166
00167 return scale[beamInfo];
00168 }
00169
00170
00171 void NCPOTCounter::DoScaling(POTMap& pot,
00172 ScaleMap& scale,
00173 Detector::Detector_t det,
00174 NCType::EDataMC dataMC) const
00175 {
00176 if(fScaleNearExposures &&
00177 det == Detector::kNear &&
00178 dataMC == NCType::kData){
00179
00180 double smallest = 1e25;
00181 for(POTMap::iterator it = pot.begin(); it != pot.end(); ++it){
00182 if(it->second < smallest)
00183 smallest = it->second;
00184 }
00185
00186 ToFixedScale(pot, scale, smallest);
00187 }
00188
00189
00190 if(fFarExposureScaleVal != 1 &&
00191 det == Detector::kFar &&
00192 dataMC == NCType::kData){
00193
00194 ToFixedScale(pot, scale, fFarExposureScaleVal);
00195 }
00196 }
00197
00198
00199 void NCPOTCounter::ToFixedScale(POTMap& pot, ScaleMap& scale,
00200 double fixedScale) const
00201 {
00202 for(POTMap::iterator itr = pot.begin(); itr != pot.end(); ++itr){
00203 scale[itr->first] = fixedScale/itr->second;
00204 itr->second = fixedScale;
00205
00206
00207 MSG("NCPOTCounter", Msg::kInfo) << "set "
00208 << itr->first.GetDescription()
00209 << " scale to "
00210 << scale[itr->first]
00211 << " - " << itr->second << " POT"
00212 << endl;
00213 }
00214 }
00215
00216
00217 vector<TString> NCPOTCounter::GetListOfFiles(int bt,
00218 NCType::EFileType fileType,
00219 Detector::Detector_t det,
00220 NCType::EDataMC dataMC) const
00221 {
00222
00223
00224 TDirectory* saveDir = gDirectory;
00225
00226 vector<TString> ret;
00227
00228 const TString fileDesignator = FileDesignator(fileType, det, dataMC);
00229
00230 const TString beamType = BeamType::AsString(fBeamIndex[bt]);
00231 const TString detName = (det == Detector::kFar) ? "far" : "near";
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 TSystemDirectory dataDir(beamType, fDataMCPath+beamType+"/");
00242 TList* files = dataDir.GetListOfFiles();
00243 assert(files);
00244
00245 for(int de = 0; de < files->GetEntries(); ++de){
00246 TString fileName = files->At(de)->GetName();
00247
00248
00249 if(fileName == "." || fileName == "..") continue;
00250
00251
00252 if(fileName.Contains(beamType) &&
00253 fileName.Contains("strip") &&
00254 fileName.Contains(detName) &&
00255 fileName.Contains(fileDesignator)){
00256
00257
00258 const int runType = NCType::FindRunType(fileName);
00259 if(dataMC == NCType::kData
00260 && fRunToUse != NCType::kRunAll
00261 && runType != fRunToUse) continue;
00262
00263 if(fileType == NCType::kMockFile && fMockDataSubRun != -1){
00264
00265 int subRun = MockFileSubRun(fileName);
00266
00267
00268
00269
00270
00271 if(subRun != -1 && subRun != fMockDataSubRun) continue;
00272 }
00273 ret.push_back(fileName);
00274 }
00275 }
00276
00277 delete files;
00278
00279 saveDir->cd();
00280
00281 return ret;
00282 }
00283
00284
00285
00286 double NCPOTCounter::GetPOTForSingleFile(TString fileName,
00287 NCType::EFileType fileType,
00288 Detector::Detector_t det,
00289 NCType::EDataMC dataMC,
00290 int bt) const
00291 {
00292
00293
00294 TFile tf(fileName);
00295
00296 const TString histName = HistName(fileType, dataMC);
00297 const TString fileDesignator = FileDesignator(fileType, det, dataMC);
00298
00299 TH1F* potHist = (TH1F*)(tf.Get(histName));
00300
00301 if(fileType==NCType::kMockFile) {
00302
00303
00304
00305 if (MockFileSubRun(fileName) != -1) return potHist->Integral()*1e12;
00306 else return POTPerFile(bt, det, fileType);
00307 }
00308
00309 MSG("NCPOTCounter", Msg::kInfo) << " Running on "
00310 << fileName << " "
00311 << potHist->GetEntries()
00312 << " files" << endl;
00313
00314
00315
00316 if(fileDesignator == "data" || fileDesignator == "blind"){
00317 double potToAdd = potHist->Integral()*1e12;
00318
00319
00320
00321
00322
00323 if(fileDesignator == "blind" && fileName.Contains("2005-06"))
00324 potToAdd += 2.56e18;
00325
00326 return potToAdd;
00327 }
00328
00329 const double potPerFile = POTPerFile(bt, det, fileType);
00330
00331 const int fileLimit = fFileLimits[det][dataMC];
00332 if(fileLimit <= potHist->GetEntries()){
00333 return potPerFile*fileLimit;
00334 }
00335 else{
00336 return potHist->Integral()*1e12;
00337 }
00338 }
00339
00340
00341 TString NCPOTCounter::FileDesignator(NCType::EFileType fileType,
00342 Detector::Detector_t det,
00343 NCType::EDataMC dataMC) const
00344 {
00345 if(fileType == NCType::kMockFile) return fMockDataSet;
00346
00347 if(dataMC == NCType::kData && !fMockOrMCAsData){
00348 if(det == Detector::kNear) return "data";
00349 else if(det == Detector::kFar) return "blind";
00350 }
00351
00352 if(fileType == NCType::kTauFile) return "tau";
00353
00354 if(fileType == NCType::kElectronFile) return "electron";
00355
00356 return "mc";
00357 }
00358
00359
00360 TString NCPOTCounter::HistName(NCType::EFileType fileType,
00361 NCType::EDataMC dataMC) const
00362 {
00363 if(fileType == NCType::kMockFile) return "dataPOT";
00364
00365 if(dataMC == NCType::kData && !fMockOrMCAsData) return "dataPOT";
00366
00367 return "mcPOT";
00368 }
00369
00370
00371 double NCPOTCounter::POTPerFile(int bt, Detector::Detector_t det,
00372 NCType::EFileType fileType) const
00373 {
00374
00375
00376
00377 if(det == Detector::kFar){
00378 if(fileType == NCType::kMockFile) {
00379 const double kPotFarMock = 2.5e20;
00380 return kPotFarMock;
00381 }
00382
00383
00384 const double kPotFarMC = 6.5e20;
00385 return kPotFarMC;
00386 }
00387
00388 const double ret = kPOTNear[fBeamIndex[bt]];
00389 assert(ret >= 0);
00390 return ret;
00391 }
00392
00393 int NCPOTCounter::MockFileSubRun(TString fileName) const
00394 {
00395
00396
00397
00398
00399
00400
00401
00402 TRegexp re("_[0-9][0-9][0-9][0-9]\\.uDST");
00403 int stupidROOT;
00404 int index = re.Index(fileName, &stupidROOT);
00405 if(index != kNPOS) {
00406
00407 TString subRunStr(fileName(index+1, 4));
00408 return subRunStr.Atoi();
00409 } else {
00410
00411 return -1;
00412 }
00413 }