00001 #include "AtNuOscFitTemplateMaker.h"
00002
00003 #include "AtNuAna/Oscillation/AtNuOscillate.h"
00004 #include "AtNuAna/SolarCycleRW/SolarCycleRW.h"
00005 #include "AtNuAna/FluxRW/AtNuFluxRW.h"
00006
00007 #include <iostream>
00008 #include <cmath>
00009
00010 ClassImp(AtNuOscFitTemplateMaker)
00011
00012 AtNuOscFitTemplateMaker::AtNuOscFitTemplateMaker()
00013 {
00014 std::cout << " *** AtNuOscFitTemplateMaker::AtNuOscFitTemplateMaker() *** " << std::endl;
00015
00016
00017
00018
00019
00020 fNsystematics = 0;
00021
00022
00023 fEnuSpec = 3.0;
00024 fSigmaSpecNu = 0.06;
00025 fSigmaSpecNuBar = 0.06;
00026 fSigmaEmuRange = 0.03;
00027 fSigmaEmuCurve = 0.05;
00028 fSigmaEshw = 0.15;
00029
00030
00031 fUseDeweightedEnergy = 1;
00032 fUseUpwardMuonMomentum = 1;
00033
00034
00035 fZenithBinningCV = 0;
00036 fZenithBinningUPMU = 0;
00037
00038
00039 fContainmentBinningCV = 0;
00040 fContainmentBinningUPMU = 0;
00041
00042
00043 fRecombineResolutionBins = 0;
00044
00045
00046 fDoProductionHeights = 1;
00047 fReDoProductionHeights = 0;
00048
00049
00050 fDoSolarCycleWeights = 1;
00051 fReDoSolarCycleWeights = 0;
00052
00053
00054 fDoFluxWeights = 1;
00055 fReDoFluxWeights = 0;
00056
00057
00058 fReweightToFluka = 0;
00059 fReweightToHonda = 0;
00060
00061
00062 fShieldSigRate = 1.0;
00063 fShieldBkgRate = 0.0;
00064
00065
00066 fNE = 25;
00067 fE1 = -0.5;
00068 fE2 = +4.5;
00069
00070
00071 fNX = 0;
00072 fX1 = 0.0;
00073 fX2 = 0.0;
00074 fX = 0;
00075 fXedge = 0;
00076 fDmsqSubBins = 0;
00077 fDmsqOverLap = 0.0;
00078
00079
00080 fChainData = new TChain("ntuple","chainData");
00081
00082 this->Reset();
00083 }
00084
00085 AtNuOscFitTemplateMaker::~AtNuOscFitTemplateMaker()
00086 {
00087 std::cout << " *** AtNuOscFitTemplateMaker::~AtNuOscFitTemplateMaker() *** " << std::endl;
00088
00089 if( fX ) delete [] fX;
00090 if( fXedge ) delete [] fXedge;
00091 if( fTemplateArray ) delete [] fTemplateArray;
00092 }
00093
00094 void AtNuOscFitTemplateMaker::Reset()
00095 {
00096 std::cout << " *** AtNuOscFitTemplateMaker::Reset() *** " << std::endl;
00097
00098
00099
00100
00101 fKtYrsData = 0.0;
00102 fTemplateType = AtNuTemplateType::kNone;
00103 fDataType = AtNuDataType::kUnknown;
00104
00105 fNsystematics = 0;
00106
00107 if( fTemplateArray ) delete [] fTemplateArray;
00108 fTemplateArray = 0;
00109
00110 fTemplateFileName = "templates.root";
00111
00112 fChainData->Reset();
00113 fChainData->SetBranchAddress("unixtime",&unixtime);
00114 fChainData->SetBranchAddress("simflag",&simflag);
00115 fChainData->SetBranchAddress("simflag",&simflag);
00116 fChainData->SetBranchAddress("index",&index);
00117 fChainData->SetBranchAddress("mc.inu",&idnu);
00118 fChainData->SetBranchAddress("mc.iact",&idact);
00119 fChainData->SetBranchAddress("mc.ires",&idres);
00120 fChainData->SetBranchAddress("mc.enu",&enu);
00121 fChainData->SetBranchAddress("mc.pnux",&pnux);
00122 fChainData->SetBranchAddress("mc.pnuy",&pnuy);
00123 fChainData->SetBranchAddress("mc.pnuz",&pnuz);
00124 fChainData->SetBranchAddress("evt.goodslice",&goodslice);
00125 fChainData->SetBranchAddress("evt.goodevent",&goodevent);
00126 fChainData->SetBranchAddress("evt.goodtrack",&goodtrack);
00127 fChainData->SetBranchAddress("evt.goodshower",&goodshower);
00128 fChainData->SetBranchAddress("evt.cv",&cv);
00129 fChainData->SetBranchAddress("evt.ce",&ce);
00130 fChainData->SetBranchAddress("evt.fc",&fc);
00131 fChainData->SetBranchAddress("evt.pc",&pc);
00132 fChainData->SetBranchAddress("evt.upmu",&upmu);
00133 fChainData->SetBranchAddress("evt.veto",&veto);
00134 fChainData->SetBranchAddress("evt.spill",&spill);
00135 fChainData->SetBranchAddress("evt.goodshield",&goodshield);
00136 fChainData->SetBranchAddress("evt.gooddirection",&gooddirection);
00137 fChainData->SetBranchAddress("evt.goodenergy",&goodenergy);
00138 fChainData->SetBranchAddress("evt.goodcharge",&goodcharge);
00139 fChainData->SetBranchAddress("evt.positivecharge",&positivecharge);
00140 fChainData->SetBranchAddress("evt.negativecharge",&negativecharge);
00141 fChainData->SetBranchAddress("evt.atmosnumu",&atmosnumu);
00142 fChainData->SetBranchAddress("evt.atmosnumucv",&atmosnumucv);
00143 fChainData->SetBranchAddress("evt.atmosnumuup",&atmosnumuup);
00144 fChainData->SetBranchAddress("evt.atmosnue",&atmosnue);
00145 fChainData->SetBranchAddress("evt.trkreco",&evttrkreco);
00146 fChainData->SetBranchAddress("evt.trkemcharge",&evttrkemcharge);
00147 fChainData->SetBranchAddress("evt.trkprange",&evttrkprange);
00148 fChainData->SetBranchAddress("evt.trkpcurve",&evttrkpcurve);
00149 fChainData->SetBranchAddress("evt.trkdiru",&evttrkdiru);
00150 fChainData->SetBranchAddress("evt.trkdirv",&evttrkdirv);
00151 fChainData->SetBranchAddress("evt.trkdirx",&evttrkdirx);
00152 fChainData->SetBranchAddress("evt.trkdiry",&evttrkdiry);
00153 fChainData->SetBranchAddress("evt.trkdirz",&evttrkdirz);
00154 fChainData->SetBranchAddress("evt.shwreco",&evtshwreco);
00155 fChainData->SetBranchAddress("evt.shwdiru",&evtshwdiru);
00156 fChainData->SetBranchAddress("evt.shwdirv",&evtshwdirv);
00157 fChainData->SetBranchAddress("evt.shwdirx",&evtshwdirx);
00158 fChainData->SetBranchAddress("evt.shwdiry",&evtshwdiry);
00159 fChainData->SetBranchAddress("evt.shwdirz",&evtshwdirz);
00160 fChainData->SetBranchAddress("evt.shwgevlin",&evtshwgevlin);
00161 fChainData->SetBranchAddress("evt.shwgevdwgt",&evtshwgevdwgt);
00162 fChainData->SetBranchAddress("evt.vtxu",&evtvtxu);
00163 fChainData->SetBranchAddress("evt.vtxv",&evtvtxv);
00164 fChainData->SetBranchAddress("evt.vtxx",&evtvtxx);
00165 fChainData->SetBranchAddress("evt.vtxy",&evtvtxy);
00166 fChainData->SetBranchAddress("evt.vtxz",&evtvtxz);
00167 fChainData->SetBranchAddress("evt.vtxplane",&evtvtxplane);
00168 fChainData->SetBranchAddress("evt.trkplanes",&evttrkplanes);
00169 fChainData->SetBranchAddress("evt.shwplanes",&evtshwplanes);
00170 fChainData->SetBranchAddress("evt.truelength",&truelength);
00171 fChainData->SetBranchAddress("evt.trueemu",&trueemu);
00172 fChainData->SetBranchAddress("evt.trueeshw",&trueeshw);
00173 fChainData->SetBranchAddress("evt.recolength",&recolength);
00174 fChainData->SetBranchAddress("evt.recoemu",&recoemu);
00175 fChainData->SetBranchAddress("evt.recoeshwlin",&recoeshwlin);
00176 fChainData->SetBranchAddress("evt.recoeshwdwgt",&recoeshwdwgt);
00177 fChainData->SetBranchAddress("evt.recoeshwnue",&recoeshwnue);
00178 fChainData->SetBranchAddress("evt.lores",&lores);
00179 fChainData->SetBranchAddress("evt.hires1",&hires1);
00180 fChainData->SetBranchAddress("evt.hires2",&hires2);
00181 fChainData->SetBranchAddress("evt.hires3",&hires3);
00182 fChainData->SetBranchAddress("evt.hires4",&hires4);
00183
00184 if( fDoSolarCycleWeights
00185 && !fReDoSolarCycleWeights ){
00186 fChainData->SetBranchAddress("solar.weight",&solarweight);
00187 }
00188
00189 if( fDoFluxWeights
00190 && !fReDoFluxWeights ){
00191 fChainData->SetBranchAddress("flux.cvupratio",&fluxcvupratio);
00192 }
00193
00194 std::cout << " entries=" << fChainData->GetEntries() << " ktyrs=" << fKtYrsData << std::endl;
00195
00196 return;
00197 }
00198
00199 void AtNuOscFitTemplateMaker::SetKtYrs( Double_t ktyrs )
00200 {
00201 std::cout << " *** AtNuOscFitTemplateMaker::SetExptKtYrs(...) *** " << std::endl;
00202
00203 fKtYrsExpt = ktyrs;
00204
00205 std::cout << " generating templates with normalization: " << fKtYrsExpt << std::endl;
00206 }
00207
00208 void AtNuOscFitTemplateMaker::SetDataType(AtNuDataType::AtNuDataType_t datatype)
00209 {
00210 std::cout << " *** AtNuOscFitTemplateMaker::SetDataType(...) *** " << std::endl;
00211
00212 if( fDataType==AtNuDataType::kUnknown ){
00213 fDataType = datatype;
00214 std::cout << " loading data with type: " << AtNuDataType::AsString(fDataType) << std::endl;
00215 }
00216 else{
00217 std::cout << " ERROR: need to specify data type " << std::endl;
00218 }
00219
00220 return;
00221 }
00222
00223 void AtNuOscFitTemplateMaker::SetTemplateType(AtNuTemplateType::AtNuTemplateType_t datatype)
00224 {
00225 std::cout << " *** AtNuOscFitTemplateMaker::SetTemplateType(...) *** " << std::endl;
00226
00227 if( fTemplateType==AtNuTemplateType::kNone ){
00228 fTemplateType = datatype;
00229 std::cout << " generating templates of type: " << AtNuTemplateType::AsString(fTemplateType) << std::endl;
00230 }
00231 else{
00232 std::cout << " ERROR: need to specify template type " << std::endl;
00233 }
00234
00235 return;
00236 }
00237
00238 void AtNuOscFitTemplateMaker::SetFluxType( AtNuFlux::FluxType_t fluxtype )
00239 {
00240 if( fluxtype==AtNuFlux::kBartol1D ){
00241
00242 }
00243
00244 if( fluxtype==AtNuFlux::kBartol3D ){
00245 fDoFluxWeights = 1;
00246 this->ReweightToBartol();
00247 }
00248
00249 if( fluxtype==AtNuFlux::kFluka3D ){
00250 fDoFluxWeights = 1;
00251 this->ReweightToFluka();
00252 }
00253
00254 if( fluxtype==AtNuFlux::kHonda3D ){
00255 fDoFluxWeights = 1;
00256 this->ReweightToHonda();
00257 }
00258
00259 return;
00260 }
00261
00262 void AtNuOscFitTemplateMaker::SetTimePeriod( Int_t mintime, Int_t maxtime )
00263 {
00264 SolarCycleRW::Instance()->TimePeriod(mintime,maxtime);
00265 }
00266
00267 void AtNuOscFitTemplateMaker::SetEnergyBinning(Int_t nE, Double_t E1, Double_t E2)
00268 {
00269 std::cout << " *** AtNuOscFitTemplateMaker::SetEnergyBinning(...) *** " << std::endl;
00270
00271 if( nE<=0 ){
00272 std::cout << " ERROR: need to specify energy binning... [return] " << std::endl;
00273 return;
00274 }
00275
00276 fNE = nE;
00277 fE1 = E1;
00278 fE2 = E2;
00279
00280 return;
00281 }
00282
00283 void AtNuOscFitTemplateMaker::SetOscillationBinning(Int_t nX, Double_t X1, Double_t X2)
00284 {
00285 std::cout << " *** AtNuOscFitTemplateMaker::SetOscillationBinning(...) *** " << std::endl;
00286
00287 if( nX<=0 ){
00288 std::cout << " ERROR: need to specify oscillation binning... [return] " << std::endl;
00289 return;
00290 }
00291
00292 fNX = nX;
00293 fX1 = X1;
00294 fX2 = X2;
00295
00296 if( fX ) delete [] fX;
00297 if( fXedge ) delete [] fXedge;
00298
00299 fX = new Double_t[fNX];
00300 fXedge = new Double_t[fNX+1];
00301
00302 if( fX1==0 ) fX1 = 1.0e-5;
00303 Double_t tempX1 = log10(fX1);
00304 Double_t tempX2 = log10(fX2);
00305
00306 if( fNX>1 ){
00307 for( Int_t i=0; i<fNX; i++ ){
00308 fX[i] = pow(10.0,tempX1+(tempX2-tempX1)*(i)/(fNX-1.0));
00309 }
00310 for( Int_t i=0; i<fNX+1; i++ ){
00311 fXedge[i] = pow(10.0,tempX1+(tempX2-tempX1)*(i-0.5)/(fNX-1.0));
00312 }
00313 }
00314 else{
00315 fX[0] = fX2;
00316 fXedge[0] = fX1; fXedge[1] = fX2;
00317 }
00318
00319 std::cout << " Dmsq(DmsqBar) Binning: " << std::endl;
00320 std::cout << " " << fNX << " " << fX1 << " " << fX2 << std::endl;
00321
00322 std::cout << " Dmsq(DmsqBar) Values: " << std::endl;
00323 for( Int_t i=0; i<fNX; i++ ){
00324 std::cout << " " << i << ": " << fX[i] << " [" << fXedge[i]<< "->" << fXedge[i+1] << "]" << std::endl;
00325 }
00326
00327 return;
00328 }
00329
00330 void AtNuOscFitTemplateMaker::AddData(AtNuDataType::AtNuDataType_t datatype, const char* file, Double_t ktyrs)
00331 {
00332 std::cout << " *** AtNuOscFitTemplateMaker::AddData(...) *** " << std::endl;
00333 std::cout << " adding file [" << AtNuDataType::AsString(datatype) << "]: "
00334 << file << " (ktyrs: " << ktyrs << ") " << std::endl;
00335
00336 if( fDataType==AtNuDataType::kUnknown ){
00337 this->SetDataType(datatype);
00338 }
00339
00340 if( fDataType!=AtNuDataType::kUnknown
00341 && fDataType==datatype ){
00342
00343 fChainData->Add(file);
00344 fKtYrsData += ktyrs;
00345
00346 std::cout << " added file: entries=" << fChainData->GetEntries() << " ktyrs=" << fKtYrsData << std::endl;
00347 }
00348 else{
00349 std::cout << " failed to add this file... " << std::endl;
00350 }
00351
00352 return;
00353 }
00354
00355 void AtNuOscFitTemplateMaker::Run()
00356 {
00357 std::cout << " *** AtNuOscFitTemplateMaker::Run() *** " << std::endl;
00358
00359 std::cout << " OscFit Template Settings:" << std::endl;
00360 std::cout << " Template Type = " << AtNuTemplateType::AsString(fTemplateType) << std::endl;
00361 std::cout << " KtYrs, Data = " << fKtYrsData << std::endl;
00362 std::cout << " KtYrs, Expt = " << fKtYrsExpt << std::endl;
00363 std::cout << " L/E Binning = [" << fNE << "," << fE1 << "," << fE2 << "]" << std::endl;
00364 std::cout << " Dmsq Binning = [" << fNX << "," << fX1 << "," << fX2 << "]" << std::endl;
00365 std::cout << " Nsystematics = " << fNsystematics << std::endl;
00366 std::cout << " DoProductionHeights = " << fDoProductionHeights << std::endl;
00367 std::cout << " ReDoProductionHeights = " << fReDoProductionHeights << std::endl;
00368 std::cout << " DoSolarCycleWeights = " << fDoSolarCycleWeights << std::endl;
00369 std::cout << " ReDoSolarCycleWeights = " << fReDoSolarCycleWeights << std::endl;
00370 std::cout << " DoFluxWeights = " << fDoFluxWeights << std::endl;
00371 std::cout << " ReDoFluxWeights = " << fReDoFluxWeights << std::endl;
00372 std::cout << " ShieldSignalRate = " << fShieldSigRate << std::endl;
00373 std::cout << " ShieldBackgroundRate = " << fShieldBkgRate << std::endl;
00374 std::cout << " RecombineResolutionBins = " << fRecombineResolutionBins << std::endl;
00375
00376
00377
00378 if( fTemplateType==AtNuTemplateType::kNone ){
00379 std::cout << " ERROR: need to select type of templates... [return]" << std::endl;
00380 return;
00381 }
00382
00383
00384
00385 if( fKtYrsData<=0
00386 && !(fTemplateType==AtNuTemplateType::kAtmosData) ){
00387 std::cout << " ERROR: need to set normalization for data... [return] " << std::endl;
00388 return;
00389 }
00390
00391
00392
00393 Int_t tempNX = 0;
00394
00395 if( fTemplateType==AtNuTemplateType::kAtmosMC
00396 || fTemplateType==AtNuTemplateType::kAtmosNumuCC ){
00397 tempNX = fNX;
00398 }
00399
00400 tempNX += 1;
00401
00402
00403
00404 std::cout << " creating templates... " << std::endl;
00405 AtNuOscFitHistogram* myHistogram =0;
00406
00407 AtNuOscFitTemplate* myTemplate = 0;
00408
00409 fTemplateArray = new TObjArray[tempNX];
00410
00411 for( Int_t nx=0; nx<tempNX; nx++ ){
00412 myTemplate = new AtNuOscFitTemplate(fTemplateType,fNE,fE1,fE2,fNsystematics);
00413 fTemplateArray[nx].Add(myTemplate);
00414 }
00415
00416
00417
00418 Double_t trkenergy,shwenergy;
00419 Double_t recoL,recoE,recoY;
00420 Double_t trueL,trueE,trueY;
00421 Double_t Dmsq,DmsqTemp,DmsqMin,DmsqMax,Sinsq;
00422 Double_t oscweight,normweight,solweight,fluxweight,vetoweight;
00423 Double_t specweight,sigmaspec;
00424 Double_t enustop,enufix;
00425 Double_t p_osc,p_osc_total;
00426 Double_t trkscale,shwscale;
00427 Double_t power,overlap;
00428 Int_t solidnu,soltime;
00429 Int_t trkcharge,trkupdn;
00430 Int_t trueUpDn,trueChg;
00431 Bool_t rangecurve;
00432
00433
00434
00435 Double_t* integratedweightNUE = new Double_t[tempNX];
00436 Double_t* integratedweightCV = new Double_t[tempNX];
00437 Double_t* integratedweightUPMU = new Double_t[tempNX];
00438
00439 for( Int_t nx=0; nx<tempNX; nx++ ){
00440 integratedweightNUE[nx] = 0.0;
00441 integratedweightCV[nx] = 0.0;
00442 integratedweightUPMU[nx] = 0.0;
00443 }
00444
00445
00446
00447 normweight = 0.0;
00448
00449 if( fKtYrsData>0.0 ){
00450 normweight = fKtYrsExpt/fKtYrsData;
00451 }
00452
00453 if( fTemplateType==AtNuTemplateType::kAtmosData ){
00454 normweight = 1.0;
00455 }
00456
00457
00458
00459 std::cout << " filling templates... " << std::endl;
00460
00461 for( Int_t n=0; n<fChainData->GetEntries(); n++ ){
00462 fChainData->GetEntry(n);
00463
00464
00465
00466 if( index==0
00467
00468
00469 && ( goodevent && !spill )
00470 && ( atmosnumu || atmosnue )
00471
00472
00473 && ( ( fDataType==AtNuDataType::kAtmos && simflag==4 && idnu!=0 )
00474 || ( fDataType==AtNuDataType::kUpMu && simflag==4 && idnu!=0 )
00475 || ( fDataType==AtNuDataType::kCosmic && simflag==4 && idnu==0 )
00476 || ( fDataType==AtNuDataType::kData && simflag==1 && idnu==0 ) )
00477
00478
00479 && ( ( fTemplateType==AtNuTemplateType::kAtmosNumuCC
00480 && simflag==4 && (fabs(idnu)==14||fabs(idnu)==16) && idact==1 )
00481 || ( fTemplateType==AtNuTemplateType::kAtmosNueCC
00482 && simflag==4 && fabs(idnu)==12 && idact==1 )
00483 || ( fTemplateType==AtNuTemplateType::kAtmosNumuNC
00484 && simflag==4 && (fabs(idnu)==14||fabs(idnu)==16) && idact==0 )
00485 || ( fTemplateType==AtNuTemplateType::kAtmosNueNC
00486 && simflag==4 && fabs(idnu)==12 && idact==0 )
00487 || ( fTemplateType==AtNuTemplateType::kAtmosMC
00488 && simflag==4 && idnu!=0 )
00489 || ( fTemplateType==AtNuTemplateType::kCosmicMC
00490 && simflag==4 && idnu==0 )
00491 || ( fTemplateType==AtNuTemplateType::kAtmosData
00492 && simflag==1 && (veto==0 || atmosnumuup) )
00493 || ( fTemplateType==AtNuTemplateType::kCosmicData
00494 && simflag==1 && (veto==1 && !atmosnumuup) ) ) ){
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 trueY = 0;
00506 trueL = 0.0;
00507 trueE = 0.0;
00508 trueUpDn = 0;
00509 trueChg = 0;
00510
00511 solweight = 1.0;
00512 fluxweight = 1.0;
00513
00514
00515 if( idnu!=0 ){
00516
00517 trueY = pnuy/enu;
00518 trueE = enu;
00519 trueL = AtNuOscillate::PropagationLength(0,idnu,-trueY,enu);
00520
00521
00522 if( fDoProductionHeights
00523 || fReDoProductionHeights ){
00524 if( fReDoProductionHeights ){
00525 trueL = AtNuOscillate::PropagationLength(2,idnu,-trueY,enu);
00526 }
00527 else{
00528 trueL = truelength;
00529 }
00530 }
00531
00532
00533 if( fDoSolarCycleWeights
00534 || fReDoSolarCycleWeights ){
00535 if( fReDoSolarCycleWeights ){
00536 solidnu = idnu;
00537 if( solidnu==+16 ) solidnu = +14;
00538 if( solidnu==-16 ) solidnu = -14;
00539
00540
00541 soltime = SolarCycleRW::Instance()->GenTime();
00542 solweight = SolarCycleRW::Instance()->Weight(enu,-trueY,
00543 solidnu,soltime);
00544 }
00545 else{
00546 solweight = solarweight;
00547 }
00548 }
00549
00550
00551 if( fDoFluxWeights
00552 || fReDoFluxWeights ){
00553
00554
00555 if( fReDoFluxWeights ){
00556 if( fDataType==AtNuDataType::kUpMu ){
00557 fluxweight = AtNuFluxRW::Instance()->GetBartolFluxRatio(idnu,enu,-trueY);
00558 }
00559 }
00560 else{
00561 if( fDataType==AtNuDataType::kUpMu ){
00562 fluxweight = fluxcvupratio;
00563 }
00564 }
00565
00566
00567 if( fReweightToFluka ){
00568 fluxweight *= AtNuFluxRW::Instance()->GetFlukaFluxRatio(idnu,enu,-trueY);
00569 }
00570 if( fReweightToHonda ){
00571 fluxweight *= AtNuFluxRW::Instance()->GetHondaFluxRatio(idnu,enu,-trueY);
00572 }
00573 }
00574
00575
00576 if( fDataType==AtNuDataType::kUpMu ){
00577 trueUpDn = +1;
00578 }
00579 else{
00580 if( pnuy>=0 ) trueUpDn = +1;
00581 else trueUpDn = -1;
00582 }
00583
00584
00585 if( idnu<0 ) trueChg = +1;
00586 else if( idnu>0 ) trueChg = -1;
00587 else trueChg = 0;
00588 }
00589
00590
00591
00592 recoY = 0.0;
00593 recoL = -999.9;
00594 recoE = 0.0;
00595 trkenergy = 0.0;
00596 shwenergy = 0.0;
00597 trkupdn = 0;
00598 trkcharge = 0;
00599 rangecurve = 1;
00600
00601 if( atmosnue ){
00602 recoY = 0.0;
00603 recoL = recolength;
00604
00605 trkenergy = recoemu;
00606 if( fUseDeweightedEnergy ) shwenergy = recoeshwdwgt;
00607 else shwenergy = recoeshwlin;
00608 recoE = trkenergy+shwenergy;
00609 }
00610
00611 if( atmosnumu ){
00612 recoY = evttrkdiry;
00613 recoL = recolength;
00614
00615 trkenergy = recoemu;
00616 if( !ce ) rangecurve = 0;
00617 if( fUseDeweightedEnergy ) shwenergy = recoeshwdwgt;
00618 else shwenergy = recoeshwlin;
00619 recoE = trkenergy+shwenergy;
00620
00621 trkupdn = 0;
00622 if( gooddirection ){
00623 if( fc || pc ){
00624 if( recoY>=0 ) trkupdn = +1; else trkupdn = -1;
00625 }
00626 else if( upmu ) trkupdn = +1;
00627 }
00628
00629 trkcharge = 0;
00630 if( goodcharge ){
00631 if( negativecharge ) trkcharge = -1;
00632 if( positivecharge ) trkcharge = +1;
00633 }
00634 }
00635
00636
00637
00638 vetoweight = 1.0;
00639
00640
00641 if( !atmosnumuup ){
00642
00643
00644 if( fTemplateType==AtNuTemplateType::kAtmosNumuCC
00645 || fTemplateType==AtNuTemplateType::kAtmosNueCC
00646 || fTemplateType==AtNuTemplateType::kAtmosNumuNC
00647 || fTemplateType==AtNuTemplateType::kAtmosNueNC
00648 || fTemplateType==AtNuTemplateType::kAtmosMC ){
00649 vetoweight = fShieldSigRate;
00650 }
00651
00652
00653 if( fTemplateType==AtNuTemplateType::kCosmicMC
00654 || fTemplateType==AtNuTemplateType::kCosmicData ){
00655 vetoweight = fShieldBkgRate;
00656 }
00657 }
00658
00659
00660
00661
00662 AtNuResolutionType::AtNuResolutionType_t resolutiontype = AtNuResolutionType::kNone;
00663
00664 if( lores ) resolutiontype = AtNuResolutionType::kLoRes;
00665 else if( hires4 ) resolutiontype = AtNuResolutionType::kHiRes4;
00666 else if( hires3 ) resolutiontype = AtNuResolutionType::kHiRes3;
00667 else if( hires2 ) resolutiontype = AtNuResolutionType::kHiRes2;
00668 else if( hires1 ) resolutiontype = AtNuResolutionType::kHiRes1;
00669
00670
00671 if( resolutiontype<=0 ){ trkupdn = 0; trkcharge = 0; }
00672
00673
00674 if( fRecombineResolutionBins ){
00675 if( !lores ) resolutiontype = AtNuResolutionType::kHiRes4;
00676 }
00677
00678
00679
00680 for( Int_t nx=0; nx<tempNX; nx++ ){
00681
00682
00683
00684
00685
00686 Dmsq = 0.0;
00687 DmsqMin = 0.0;
00688 DmsqMax = 0.0;
00689 Sinsq = 0.0;
00690
00691 if( fTemplateType==AtNuTemplateType::kAtmosMC
00692 || fTemplateType==AtNuTemplateType::kAtmosNumuCC ){
00693 if( fNX>0 && nx>0 ){
00694 Dmsq = fX[nx-1];
00695 DmsqMin = fXedge[nx-1];
00696 DmsqMax = fXedge[nx];
00697 Sinsq = 1.0;
00698 }
00699 }
00700
00701
00702
00703 oscweight = 1.0;
00704 if( fabs(idnu)==16 ) oscweight = 0.0;
00705
00706 if( ( idact==1 )
00707 && ( fabs(idnu)==14 || fabs(idnu)==16 ) ){
00708
00709 p_osc = 0.0;
00710
00711 if( Dmsq>0.0 && DmsqMax>=DmsqMin ){
00712 if( fDmsqSubBins<=1 ){
00713
00714 DmsqTemp = Dmsq;
00715 p_osc = Sinsq*pow(sin(1.267*DmsqTemp*trueL/trueE),2.0);
00716
00717 }
00718 else{
00719
00720 p_osc = 0.0;
00721 p_osc_total = 0.0;
00722
00723 for( Int_t idmsq=0; idmsq<fDmsqSubBins; idmsq++ ){
00724
00725
00726
00727
00728 overlap = 1.0+fDmsqOverLap;
00729 power = 0.5*(log10(DmsqMin)+log10(DmsqMax))
00730 + 0.5*overlap*(1.0-2.0*(idmsq)/(fDmsqSubBins-1.0))*(log10(DmsqMin)-log10(DmsqMax));
00731 DmsqTemp = pow(10.0,power);
00732
00733 p_osc += Sinsq*pow(sin(1.267*DmsqTemp*trueL/trueE),2.0);
00734 p_osc_total += 1.0;
00735 }
00736
00737 p_osc = p_osc/p_osc_total;
00738 }
00739 }
00740
00741 if( fabs(idnu)==14 && idact==1 ){
00742 oscweight = 1.0 - p_osc;
00743 }
00744
00745 if( fabs(idnu)==16 && idact==1 ){
00746 oscweight = p_osc;
00747 }
00748 }
00749
00750
00751
00752 myTemplate = (AtNuOscFitTemplate*)(fTemplateArray[nx].At(0));
00753
00754
00755 for( Int_t isys=0; isys<=fNsystematics; isys++ ){
00756 for( Int_t isigma=-2; isigma<=2; isigma++ ){
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766 if( ( isys==0 && isigma==0 )
00767 || ( isys>0 && isigma!=0 ) ){
00768
00769
00770 specweight = 1.0;
00771 trkscale = 0.0;
00772 shwscale = 0.0;
00773
00774 if( idnu!=0 && isigma!=0 ){
00775
00776 sigmaspec = 0.0;
00777 if( isys==1 && idnu>0 ){
00778 sigmaspec = fSigmaSpecNu;
00779 }
00780
00781 if( isys==2 && idnu<0 ){
00782 sigmaspec = fSigmaSpecNuBar;
00783 }
00784
00785
00786 if( fEnuSpec>0.0 && sigmaspec>0.0 ){
00787 if( enu>0 && enu<=fEnuSpec ){
00788 specweight = 1.0 + isigma*(sigmaspec)*(enu-fEnuSpec);
00789
00790
00791 }
00792 else if( enu>fEnuSpec ){
00793 enustop = fEnuSpec*exp(1.0/(2.0*sigmaspec));
00794 enufix = enustop*(1.0-exp(-enu/enustop));
00795
00796 specweight = 1.0 + isigma*(sigmaspec)*(log(enufix/fEnuSpec));
00797 }
00798 }
00799
00800 if( isys==3 ){
00801 if( ce ) trkscale = isigma*fSigmaEmuRange;
00802 }
00803
00804 if( isys==4 ){
00805 if( !ce ) trkscale = isigma*fSigmaEmuCurve;
00806 }
00807
00808 if( isys==5 ){
00809 shwscale = isigma*fSigmaEshw;
00810 }
00811
00812 if( specweight<0.0 ) specweight = 0.0;
00813 if( trkscale<-1.0 ) trkscale = -1.0;
00814 if( shwscale<-1.0 ) shwscale = -1.0;
00815
00816 }
00817
00818
00819 Double_t weight = normweight*oscweight*specweight*solweight*fluxweight*vetoweight;
00820 Double_t fillE = recoE + trkscale*trkenergy + shwscale*shwenergy;
00821 Double_t fillL = recoL;
00822 Double_t fillY = recoY;
00823 Double_t fillLogLE = 0.0;
00824
00825
00826 Double_t fillEmu = trkenergy + trkscale*trkenergy;
00827
00828
00829 if( atmosnue ){
00830 myHistogram = (AtNuOscFitHistogram*)(myTemplate->GetHistogram(AtNuResolutionEventType::kNUE,resolutiontype,trueUpDn,trueChg,isys,isigma));
00831 if( myHistogram ){
00832 myHistogram->Fill( fillLogLE, weight, 0, 0 );
00833 if( isys==0 ) integratedweightNUE[nx] += weight;
00834 }
00835 }
00836
00837
00838 if( atmosnumu && (fc||pc) ){
00839 if( resolutiontype>0 ){
00840 if( fZenithBinningCV ){
00841 fillLogLE = 0.5*(fillY*(fE1-fE2)+(fE1+fE2));
00842 }
00843 else{
00844 fillLogLE = log10(fillL/fillE);
00845 }
00846 }
00847
00848 myHistogram = (AtNuOscFitHistogram*)(myTemplate->GetHistogram(AtNuResolutionEventType::kCV,resolutiontype,trueUpDn,trueChg,isys,isigma));
00849 if( myHistogram ){
00850 myHistogram->Fill( fillLogLE, weight, trkupdn, trkcharge );
00851 if( isys==0 ) integratedweightCV[nx] += weight;
00852 }
00853 }
00854
00855
00856 if( atmosnumu && upmu ){
00857 if( resolutiontype>0 ){
00858 if( fContainmentBinningUPMU ){
00859 AtNuResolutionType::AtNuResolutionType_t resolutionfix = AtNuResolutionType::kNone;
00860 if( ce ) resolutionfix = AtNuResolutionType::kHiRes4;
00861 else resolutionfix = AtNuResolutionType::kHiRes3;
00862 resolutiontype = resolutionfix;
00863 }
00864
00865 if( fZenithBinningUPMU ){
00866 fillLogLE = 0.5*(fillY*(fE1-fE2)+(fE1+fE2));
00867 }
00868 else{
00869 if( fUseUpwardMuonMomentum ) fillLogLE = log10(fillL/fillEmu);
00870 else fillLogLE = log10(fillL/fillE);
00871 }
00872 }
00873
00874 myHistogram = (AtNuOscFitHistogram*)(myTemplate->GetHistogram(AtNuResolutionEventType::kUPMU,resolutiontype,trueUpDn,trueChg,isys,isigma));
00875 if( myHistogram ){
00876 myHistogram->Fill( fillLogLE, weight, trkupdn, trkcharge );
00877 if( isys==0 ) integratedweightUPMU[nx] += weight;
00878 }
00879 }
00880
00881 }
00882
00883 }
00884 }
00885
00886 }
00887
00888 }
00889 }
00890
00891
00892
00893
00894 for( Int_t nx=0; nx<tempNX; nx++ ){
00895
00896 myTemplate = (AtNuOscFitTemplate*)(fTemplateArray[nx].At(0));
00897
00898
00899
00900 if( fNsystematics>1 ){
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955 }
00956
00957 }
00958
00959
00960
00961 std::cout << " total normalizations... " << std::endl;
00962 for( Int_t nx=0; nx<tempNX; nx++ ){
00963 std::cout << " " << nx << ": nue=" << integratedweightNUE[nx] << " cv=" << integratedweightCV[nx] << " upmu=" << integratedweightUPMU[nx] << std::endl;
00964 }
00965
00966
00967
00968 for( Int_t nx=0; nx<tempNX; nx++ ){
00969
00970 TString filename(fTemplateFileName.Data());
00971
00972 if( tempNX>1 ){
00973 TString tempfilename("");
00974 tempfilename.Append(".");
00975 if( nx<100 ) tempfilename.Append("0");
00976 if( nx<10 ) tempfilename.Append("0");
00977 tempfilename+=nx;
00978 if( filename.EndsWith(".root") ){
00979 filename.Insert(filename.Length()-5,tempfilename);
00980 }
00981 else{
00982 filename.Append(tempfilename);
00983 }
00984 }
00985
00986 myTemplate = (AtNuOscFitTemplate*)(fTemplateArray[nx].At(0));
00987 myTemplate->WriteHistograms(filename.Data());
00988
00989 }
00990
00991 delete [] integratedweightNUE;
00992 delete [] integratedweightCV;
00993 delete [] integratedweightUPMU;
00994
00995 return;
00996 }
00997
00998
00999