loon: AtNuOscFitTemplateMaker.cxx Source File

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   // initializations
00017   // ===============
00018 
00019   // number of systematics
00020   fNsystematics = 0;      // initialize as zero
00021 
00022   // systematic parameters
00023   fEnuSpec = 3.0;          // spectrum switching point 
00024   fSigmaSpecNu = 0.06;     // error in spectral index (neutrinos)
00025   fSigmaSpecNuBar = 0.06;  // error in spectral index (anti-neutrinos)
00026   fSigmaEmuRange = 0.03;   // error in muon momentum (from range)
00027   fSigmaEmuCurve = 0.05;   // error in muon momentum (from curvature)
00028   fSigmaEshw = 0.15;       // error in shower energy
00029 
00030   // energy settings
00031   fUseDeweightedEnergy = 1;    // use deweighted (not linear) shower energy
00032   fUseUpwardMuonMomentum = 1;  // ignore shower energy for upward muons
00033 
00034   // bin by zenith angle
00035   fZenithBinningCV = 0;
00036   fZenithBinningUPMU = 0; 
00037 
00038   // bin by containment
00039   fContainmentBinningCV = 0;
00040   fContainmentBinningUPMU = 0;
00041 
00042   // recombine resolution bins
00043   fRecombineResolutionBins = 0;
00044 
00045   // production heights
00046   fDoProductionHeights = 1;
00047   fReDoProductionHeights = 0;
00048 
00049   // solar cycle weighting
00050   fDoSolarCycleWeights = 1;
00051   fReDoSolarCycleWeights = 0;
00052 
00053   // flux reweighting
00054   fDoFluxWeights = 1;
00055   fReDoFluxWeights = 0;
00056 
00057   // switch flux model
00058   fReweightToFluka = 0;
00059   fReweightToHonda = 0;
00060 
00061   // veto shield
00062   fShieldSigRate = 1.0;
00063   fShieldBkgRate = 0.0;
00064 
00065   // default energy binning: 25 bins of log(L/E)
00066   fNE = 25;
00067   fE1 = -0.5;
00068   fE2 = +4.5;
00069 
00070   // oscillation bins - start with no binning
00071   fNX = 0;                // number of bins
00072   fX1 = 0.0;              // first bin
00073   fX2 = 0.0;              // last bin
00074   fX = 0;                 // bin centres [array]
00075   fXedge = 0;             // bin edges   [array]
00076   fDmsqSubBins = 0;       // sub-bins
00077   fDmsqOverLap = 0.0;     // overlap
00078 
00079   // inputted data
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   // re-initializations
00099   // ==================
00100 
00101   fKtYrsData = 0.0;       // normalization (inputted data)
00102   fTemplateType = AtNuTemplateType::kNone;  // template type
00103   fDataType = AtNuDataType::kUnknown;       // data type
00104 
00105   fNsystematics = 0;      // number of systematics
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     // not implemented
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   // need a data type
00377   // ================
00378   if( fTemplateType==AtNuTemplateType::kNone ){
00379     std::cout << "   ERROR: need to select type of templates... [return]" << std::endl;
00380     return;
00381   }
00382     
00383   // need normalization for inputted data
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   // oscillation grid
00392   // ================
00393   Int_t tempNX = 0;
00394 
00395   if( fTemplateType==AtNuTemplateType::kAtmosMC
00396    || fTemplateType==AtNuTemplateType::kAtmosNumuCC ){
00397     tempNX = fNX;
00398   }
00399 
00400   tempNX += 1; // real data or null oscillations monte carlo
00401 
00402   // create array of histograms
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   // analysis variables
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   // total normalizations
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   // normalizations
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   // Loop over Data
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     // event selection
00465     // ===============
00466     if( index==0
00467 
00468      // selected event
00469      && ( goodevent && !spill )         
00470      && ( atmosnumu || atmosnue )
00471 
00472      // correct data type
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      // correct template type
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       // Production heights
00497       // ==================
00498       // (calculation moved to ntuple maker)
00499       // Double_t fR = 6370.0;  // radius of earth
00500       // Double_t fH = 15.0;    // production height
00501       // Double_t fD = 1.0;     // detector depth
00502       
00503       // Monte Carlo truth
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       // neutrino weights
00515       if( idnu!=0 ){
00516         
00517         trueY = pnuy/enu;
00518         trueE = enu;
00519         trueL = AtNuOscillate::PropagationLength(0,idnu,-trueY,enu);
00520 
00521         // production heights
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         // solar cycle reweighting
00533         if( fDoSolarCycleWeights 
00534          || fReDoSolarCycleWeights ){
00535           if( fReDoSolarCycleWeights ){
00536             solidnu = idnu;
00537             if( solidnu==+16 ) solidnu = +14;
00538             if( solidnu==-16 ) solidnu = -14;  
00539             // soltime = unixtime;    // set to random time during run period
00540             // soltime = 1136073600;  // set to mid-point (1st January 2006 for now)
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         // flux weighting
00551         if( fDoFluxWeights
00552          || fReDoFluxWeights ){
00553 
00554           // Bartol96 to Bartol3D
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           // Bartol3D to Fluka3D or Honda3D
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         // true direction
00576         if( fDataType==AtNuDataType::kUpMu ){
00577           trueUpDn = +1; 
00578         }
00579         else{
00580           if( pnuy>=0 ) trueUpDn = +1; 
00581           else trueUpDn = -1;  
00582         }
00583 
00584         // true charge sign
00585         if( idnu<0 ) trueChg = +1; 
00586           else if( idnu>0 ) trueChg = -1; 
00587             else trueChg = 0;
00588       }
00589 
00590       // event kinematics
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 ){  // electron neutrino events
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 ){  // muon neutrino events
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       // shield efficiency
00637       // =================
00638       vetoweight = 1.0;
00639 
00640       // apply to containment selection 
00641       if( !atmosnumuup ){
00642 
00643         // atmospheric neutrinos
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         // cosmic muons
00653         if( fTemplateType==AtNuTemplateType::kCosmicMC 
00654          || fTemplateType==AtNuTemplateType::kCosmicData ){
00655           vetoweight = fShieldBkgRate;
00656         }
00657       }
00658 
00659       // event resolution
00660       // ================
00661       // (use external resolution class)
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       // (use low resolution for normalization only)
00671       if( resolutiontype<=0 ){ trkupdn = 0; trkcharge = 0; }
00672 
00673       // (combine resolution bins into a single bin) 
00674       if( fRecombineResolutionBins ){
00675         if( !lores ) resolutiontype = AtNuResolutionType::kHiRes4;
00676       }
00677 
00678       // loop over oscillation bins
00679       // ==========================
00680       for( Int_t nx=0; nx<tempNX; nx++ ){    
00681             
00682         // oscillation parameters
00683         // ======================
00684         // either use centre of bin, or average fDmsqSubBins sub-bins
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         // oscillation weight
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                 // DmsqTemp = pow(10.0, log10(DmsqMin)
00726                 //               + (log10(DmsqMax)-log10(DmsqMin))*(idmsq)/(fDmsqSubBins-1.0));
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 ){ // numu CC
00742             oscweight = 1.0 - p_osc;
00743           }
00744 
00745           if( fabs(idnu)==16 && idact==1 ){ // nutau CC
00746             oscweight = p_osc;
00747           }
00748         }
00749 
00750         // fill histograms
00751         // ===============
00752         myTemplate = (AtNuOscFitTemplate*)(fTemplateArray[nx].At(0));
00753 
00754         // loop over systematics
00755         for( Int_t isys=0; isys<=fNsystematics; isys++ ){
00756           for( Int_t isigma=-2; isigma<=2; isigma++ ){ 
00757 
00758             // systematic shifts
00759             // 0: nominal - no shift
00760             // 1: SpecNu    (spectral index, neutrinos)
00761             // 2: SpecNuBar (spectral index, anti-neutrinos)
00762             // 3: EmuRange  (muon momentum from range)
00763             // 4: EmuCurve  (muon momentum from curvature)
00764             // 5: Eshw      (shower energy)
00765 
00766             if( ( isys==0 && isigma==0 )      // nominal spectrum
00767              || ( isys>0  && isigma!=0 ) ){   // shifted spectra
00768 
00769               // systematic shifts
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 ){  // isys=1: spectral index (neutrinos)
00778                   sigmaspec = fSigmaSpecNu;
00779                 }
00780 
00781                 if( isys==2 && idnu<0 ){  // isys=2: spectral index (anti-neutrinos)
00782                   sigmaspec = fSigmaSpecNuBar;
00783                 }
00784 
00785                 // spectral index systematic
00786                 if( fEnuSpec>0.0 && sigmaspec>0.0 ){
00787                   if( enu>0 && enu<=fEnuSpec ){ 
00788                     specweight = 1.0 + isigma*(sigmaspec)*(enu-fEnuSpec);
00789                     // old function (which preserved first derivative):
00790                     //  specweight = 1.0 + isigma*(sigmaspec/fEnuSpec)*(enu-fEnuSpec);
00791                   }
00792                   else if( enu>fEnuSpec ){
00793                     enustop = fEnuSpec*exp(1.0/(2.0*sigmaspec)); // a fix to keep
00794                     enufix = enustop*(1.0-exp(-enu/enustop));    // -1<weight<+1
00795 
00796                     specweight = 1.0 + isigma*(sigmaspec)*(log(enufix/fEnuSpec));
00797                   }
00798                 }
00799 
00800                 if( isys==3 ){  // isys=3: track energy scale
00801                   if( ce ) trkscale = isigma*fSigmaEmuRange;  
00802                 }
00803 
00804                 if( isys==4 ){  // isys=4: track energy scale
00805                   if( !ce ) trkscale = isigma*fSigmaEmuCurve;
00806                 }
00807 
00808                 if( isys==5 ){  // isys=5: shower energy scale
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               // fill histograms
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               // separate muon energy
00826               Double_t fillEmu = trkenergy + trkscale*trkenergy; 
00827 
00828               // electron neutrinos
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               // contained muon neutrinos
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               // rock muon neutrinos
00856               if( atmosnumu && upmu ){
00857                 if( resolutiontype>0 ){
00858                   if( fContainmentBinningUPMU ){ // change binning to containment rather than momentum
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             } // if( ( isys==0 && isigma==0 ) || ( isys>0 && isigma!=0 ) ...
00882 
00883           } // for( Int_t isys=0; isys<=fNsystematics; isys++ ){ ...
00884         } // for( Int_t isigma=-2; isigma<=2; isigma++ ){ ...
00885 
00886       } // for( Int_t nx=0; nx<tempNX; nx++ ){  // Dmsq
00887 
00888     } // if( index==0 && ... // event selection
00889   } // for( Int_t n=0... // loop over events
00890 
00891 
00892   // re-normalize templates
00893   // ======================
00894   for( Int_t nx=0; nx<tempNX; nx++ ){  
00895 
00896     myTemplate = (AtNuOscFitTemplate*)(fTemplateArray[nx].At(0));
00897 
00898     // systematic shifts
00899     // =================
00900     if( fNsystematics>1 ){
00901   
00902       // 0: nominal - no shift
00903       // 1: SpecNu    (spectral index, neutrinos)
00904       // 2: SpecNuBar (spectral index, anti-neutrinos)
00905 
00906       // Double_t nominal_weight = 0.0;
00907       // Double_t systematic_weight = 0.0;
00908 
00909       //
00910       // spectral index, neutrinos
00911       // =========================  
00912       // nominal_weight = 0.0;
00913       // for( Int_t idir=-1; idir<=+1; idir++ ){
00914       //   nominal_weight += myTemplate->GetTotalContent(idir,-1,0,0);
00915       // }
00916       //
00917       // if( nominal_weight>0.0 ){
00918       //   for( Int_t isigma=-2; isigma<=2; isigma++ ){
00919       //     systematic_weight = 0.0;
00920       //     for( Int_t idir=-1; idir<=+1; idir++ ){
00921       //       systematic_weight += myTemplate->GetTotalContent(idir,-1,1,isigma);
00922       //     }
00923       //     if( systematic_weight>0.0 ){
00924       //       for( Int_t idir=-1; idir<=+1; idir++ ){
00925       //         myTemplate->Scale(idir,-1,1,isigma,nominal_weight/systematic_weight);
00926       //       }
00927       //     }
00928       //   }
00929       // }
00930       //
00931       
00932       //
00933       // spectral index, anti-neutrinos
00934       // ==============================  
00935       // nominal_weight = 0.0;
00936       // for( Int_t idir=-1; idir<=+1; idir++ ){
00937       //   nominal_weight += myTemplate->GetTotalContent(idir,+1,0,0);
00938       // }
00939       //
00940       // if( nominal_weight>0.0 ){
00941       //   for( Int_t isigma=-2; isigma<=2; isigma++ ){
00942       //     systematic_weight = 0.0;
00943       //     for( Int_t idir=-1; idir<=+1; idir++ ){
00944       //       systematic_weight += myTemplate->GetTotalContent(idir,+1,2,isigma);
00945       //     }
00946       //     if( systematic_weight>0.0 ){
00947       //       for( Int_t idir=-1; idir<=+1; idir++ ){
00948       //         myTemplate->Scale(idir,+1,2,isigma,nominal_weight/systematic_weight);
00949       //       }
00950       //     }
00951       //   }
00952       // }
00953       //
00954 
00955     } // if( fNsystematics>1 ){ ...
00956 
00957   } // for( Int_t nx=0; nx<tempNX; nx++ ){ // Dmsq
00958 
00959   // print out normalizations
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   // write out templates
00967   // ===================
00968   for( Int_t nx=0; nx<tempNX; nx++ ){    // Dmsq
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