Changeset 6486 for trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc

Timestamp:

02/15/05 15:36:45 (20 years ago)

Author:

rico

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc (modified) (16 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc

@@ -47 +47 @@
 using namespace std;
 
-const Int_t fNTbins = 2;
-const Double_t fTbin[fNTbins+1] = {0,10,20};
-
+const Int_t ntbins=1;                      // default number of theta bins
+const Double_t tbin[ntbins+1] = {0,90};    // default theta bins bounds
+const Int_t nebins = 10;                   // default number of energy bins
+const Float_t emin = 10.;                  // default minimum energy value
+const Float_t emax = 10000.;               // default maximum energy value
+const Int_t nsubbins = 20;                 // default number of subbins per bin
+const Char_t* deff = "4.e9*pow(x,-2.6+1)"; // default spectrum function
 
 // -------------------------------------------------------------------------
@@ -56 +60 @@
 //
 MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
-  : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
-    fCoeff(NULL), fEffA(NULL), fFile(NULL)
+  : fSpec(NULL),  fEmin(emin), fEmax(emax), fEbins(nebins), fEsubbins(nsubbins), 
+    fWeight(NULL), fCoeff(NULL), fEffA(NULL)//, fFile(NULL)
 {
   // set the default function
-  SetFunction("4.e9*pow(x,-2.6+1)");  
+  SetFunction(deff);  
 
   // create the TChains 
-  // OJO!! make sure they read the appropriate branch
   fCini = new TChain("OriginalMC");
   fCcut = new TChain("Events");
@@ -76 +79 @@
   fCcut->SetBranchAddress("MHillas.",&fHillas);
   fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
+
+  // initial value of the zenith angle binning
+  SetThetaBinning(ntbins,tbin);
   
   // borra
-  fFile = new TFile("test.root","RECREATE");
+  //  fFile = new TFile("coeftest.root","RECREATE");
 }
 
@@ -90 +96 @@
     delete fSpec;
 
+  if(fTbin)
+    delete [] fTbin;
+
   if(fWeight)
     delete [] fWeight;
@@ -99 +108 @@
     delete fEffA;
 
-  if(fHorig) 
-    delete fHorig;
 
   delete fCini;
   delete fCcut;
+}
+
+// -------------------------------------------------------------------------
+//
+// Set the binning in zenith angle
+//
+void MEffAreaAndCoeffCalc::SetThetaBinning(Int_t n, const Double_t* binlist)
+{
+  fNTbins=n;
+  if(fTbin) delete [] fTbin;
+  fTbin = new Double_t[n+1];
+  for(Int_t i=0;i<n+1;i++)
+    fTbin[i] = binlist[i];
 }
 
@@ -144 +164 @@
 void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
 {  
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
-  
-  if(fHorig) delete fHorig;
-  fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
-  fCini->Draw("logenergy>>fhorig","","goff");
-  // borra
-  fHorig->Write();
 }
 
@@ -161 +172 @@
 void MEffAreaAndCoeffCalc::ComputeWeights()
 {  
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
-  const Double_t desub = de/fEsubbins; // subbin size (in log)
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
-
+  // OJO!! maybe this should be hard-coded somewhere else
+  const Float_t abslogmin=0.;
+  const Float_t abslogmax=5.;
+
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_t logemax = TMath::Log10(fEmax);
+  const Float_t de    = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
+  const Int_t esbins   = fEbins*fEsubbins; // total number of subbins
+
+  // compute the binning for weights histogram
+  const Int_t nmindist = (logemin>abslogmin)? Int_t((logemin-abslogmin)/desub) : 0;
+  const Int_t nmaxdist = (logemax<abslogmax)? Int_t((abslogmax-logemax)/desub) : 0;
+
+  fLogEWmin = logemin-desub*nmindist;
+  fLogEWmax = logemax+desub*nmaxdist;
+  fEWbins   = nmindist+esbins+nmaxdist;
+  
   // reset the weights array
   if(fWeight)
     delete [] fWeight;
-  fWeight = new Double_t[esbins];
+  fWeight = new Double_t[fEWbins];
     
-  if(!fHorig)
-    FillOriginalSpectrum();
-  
-  // borra
-  TH1F hw("hw","hw",esbins,fEmin,fEmax);
-  for(Int_t i=0;i<esbins;i++)
-    {
-      const Double_t denom = fHorig->GetBinContent(i+1);              // number of events in initial spectrum
-      const Double_t ew    = TMath::Power(10,logemin+(i+0.5)*desub); // real energy
-      const Double_t numer = fSpec->Eval(ew);                        // number of events for the required spectrum
-      if(denom)
+
+  TH1F* horigs = new TH1F("horigs","Original energy spectrum",fEWbins,fLogEWmin,fLogEWmax);
+  fCini->Draw("logenergy>>horigs","","goff");
+  // borra
+  //  horigs->Write();
+
+  // borra
+  //  TH1F hw("hw","hw",fEWbins,fLogEWmin,fLogEWmax);
+  for(Int_t i=0;i<fEWbins;i++)
+    {
+      const Float_t denom = horigs->GetBinContent(i+1);               // number of events in initial spectrum
+      const Float_t ew    = TMath::Power(10,fLogEWmin+(i+0.5)*desub); // real energy
+      const Float_t numer = fSpec->Eval(ew);                          // number of events for the required spectrum
+      if(denom>10)
         fWeight[i]=numer/denom;
       else
         {
-          cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
+          //      cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
           fWeight[i]=-1;
         }
       // borra
-      hw.SetBinContent(i+1,fWeight[i]);
-    }
-  // borra
-  hw.Write();
+      //      hw.SetBinContent(i+1,fWeight[i]);
+    }
+  // borra
+  //  hw.Write();
+
+  delete horigs;
 }
 
@@ -208 +236 @@
     }
 
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
-  const Double_t desub = de/fEsubbins; // subbin size (in log)
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_t logemax = TMath::Log10(fEmax);
+  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
   const Int_t nentries = Int_t(fCcut->GetEntries());
  
@@ -220 +247 @@
 
   // borra
-  TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
-  TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
-  TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
+  //  TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
+  //  TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
+  //  TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
 
   // reset the coefficients
@@ -238 +265 @@
       Float_t emc   = fMcEvt->GetEnergy();
       Float_t estim = fHillas->GetSize()/0.18/15.;
+
+      if((emc<fEmin && estim<fEmin) || (emc>fEmax && estim>fEmax) ||
+         (emc<fEmin && estim>fEmax) || (emc>fEmax && estim<fEmin))
+        continue;
+
       Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
 
       // compute the bin where the weight is taken from
-      Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub);
+      Int_t wbin = Int_t((TMath::Log10(emc)-fLogEWmin)/desub);
 
       // fill the histograms
-      if(wbin<esbins)
+      if(wbin<fEWbins)
         {
-          hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
-          hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
-          // borra
-          if(theta<fTbin[1])
+          if(fWeight[wbin]>0)
             {
-              hest1->Fill(TMath::Log10(estim),fWeight[wbin]); 
-              hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
+              hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
+              hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
+              // borra
+//            if(theta<fTbin[1])
+//              {
+//                hest1->Fill(TMath::Log10(estim),fWeight[wbin]); 
+//                hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
+//              }
             }
+          else
+            cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has no computed weight (lack of MC statistics), skipping" << endl;
         }
       else
-        cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
+        cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
     }
 
@@ -266 +303 @@
         const Float_t den = hest->GetBinContent(i+1,j+1);
         //borra
-        const Float_t num1 = hmc1->GetBinContent(i+1);
-        const Float_t den1 = hest1->GetBinContent(i+1);
+//      const Float_t num1 = hmc1->GetBinContent(i+1);
+//      const Float_t den1 = hest1->GetBinContent(i+1);
         if(den)
           {
             fCoeff->SetBinContent(i+1,j+1,num/den);
             //borra
-            if(j==0 && den1)
-              coef1->SetBinContent(i+1,num1/den1);
+//          if(j==0 && den1)
+//            coef1->SetBinContent(i+1,num1/den1);
           }
         else
@@ -284 +321 @@
 
   //borra
-  hmc1->Write();
-  hest1->Write();
-  coef1->Write();
+//   hmc1->Write();
+//   hest1->Write();
+//   coef1->Write();
 
   delete hmc;
@@ -305 +342 @@
 
   //OJO!! radius should be read from somewhere!
-  const Double_t radius = 30000.; // generation radius (in cm)
-  //  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
-  const Double_t Atot = 1.; //total area (in cm^2)
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
+  const Float_t radius = 30000.; // generation radius (in cm)
+  const Float_t Atot = 3.14159*radius*radius; //total area (in cm^2)
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_t logemax = TMath::Log10(fEmax);
   const Int_t esbins = fEbins*fEsubbins; // total number of subbins
-  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
-  const Double_t desub = de/fEsubbins; // subbin size (in log)
+  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
 
   // reset the effective areas
@@ -319 +355 @@
   fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
   //borra
-  TH1F eff("eff","Effective area",fEbins,logemin,logemax);
+//   TH1F eff("eff","Effective area",fEbins,logemin,logemax);
   
   // fill the spectrum of the survivors
@@ -343 +379 @@
             else
               {
-                const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
-                const Double_t ww = fSpec->Eval(ew);
+                const Float_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
+                const Float_t ww = fSpec->Eval(ew);
                 effarea+=Atot*numerator/denominator*ww;
                 wtot   += ww;
@@ -357 +393 @@
           {
             fEffA->SetBinContent(i+1,j+1,effarea/wtot); 
-            if(j==0)
-              {
-                cout << "*****" << i << ": " << effarea/wtot << endl;
-                eff.SetBinContent(i+1,effarea/wtot);
-              }
+            // borra
+//          if(j==0)
+//            {
+//              //              cout << "*****" << i << ": " << effarea/wtot << endl;
+//              eff.SetBinContent(i+1,effarea/wtot);
+//            }
           }
       }
 
   // borra
-  eff.Write();
+//   eff.Write();
   
   delete hpass;