Changeset 6465 for trunk/MagicSoft

Timestamp:

02/14/05 19:04:58 (20 years ago)

Author:

rico

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae

Files:

• Changelog (modified) (1 diff)
• library/MEffAreaAndCoeffCalc.cc (modified) (13 diffs)
• library/MEffAreaAndCoeffCalc.h (modified) (4 diffs)
• macros/computeCoeff.C (modified) (2 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/Changelog

@@ -19 +19 @@
                                                  -*-*- END OF LINE -*-*-
 
+ 2005/02/10 Javier Rico
+    * library/MEffAreaAndCoeffCalc.[h,cc]
+    * macros/computeCoeff.C
+      - compute the coefficients in the big bins
+      - migrate to new data format
+      
+        
  2005/02/10 Javier Rico
     * library/MEffAreaAndCoeffCalc.[h,cc]

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

@@ -38 +38 @@
 #include "MMcEvt.hxx"
 #include "TH2F.h"
+#include "TFile.h"
 
 #include "MLog.h"
@@ -55 +56 @@
 //
 MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
-  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.),
-    fCoeff(NULL), fEffA(NULL)
+  : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
+    fCoeff(NULL), fEffA(NULL), fFile(NULL)
 {
   // set the default function
-  // OJO!! this should be changed
   SetFunction("4.e9*pow(x,-2.6+1)");  
 
   // create the TChains 
   // OJO!! make sure they read the appropriate branch
-  fCini = new TChain("Events");
+  fCini = new TChain("OriginalMC");
   fCcut = new TChain("Events");
 
-  // define some usefull aliases
-  fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
-  fCini->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
-
-  // OJO!! the estimated energy must be changed for a real one
+  // define some useful aliases
+  fCini->SetAlias("logenergy","log10(MMcEvtBasic.fEnergy)");
+  fCini->SetAlias("theta","MMcEvtBasic.fTelescopeTheta*180./3.14159");
+
   fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
-  fCcut->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
-
-  fCcut->SetBranchStatus("*",0);
-  fCcut->SetBranchStatus("MHillas.*",1);
-  fCcut->SetBranchStatus("MMcEvt.*",1);
+  fCcut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
+
   fCcut->SetBranchAddress("MHillas.",&fHillas);
   fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
+  
+  // borra
+  fFile = new TFile("test.root","RECREATE");
 }
 
@@ -119 +118 @@
       emin = fEmin;
       emax = fEmax;
+    }
+  else
+    {
+      fEmin = emin;
+      fEmax = emax;
     }
   fSpec = new TF1("fspec",chfunc,emin,emax);
@@ -147 +151 @@
   fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
   fCini->Draw("logenergy>>fhorig","","goff");
+  // borra
+  fHorig->Write();
 }
 
@@ -169 +175 @@
     FillOriginalSpectrum();
   
+  // borra
+  TH1F hw("hw","hw",esbins,fEmin,fEmax);
   for(Int_t i=0;i<esbins;i++)
     {
@@ -181 +189 @@
           fWeight[i]=-1;
         }
-    }
+      // borra
+      hw.SetBinContent(i+1,fWeight[i]);
+    }
+  // borra
+  hw.Write();
 }
 
@@ -207 +219 @@
   TH2F* hmc  = new TH2F("hmc","MC energy",fEbins,logemin,logemax,fNTbins,fTbin);
 
+  // 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);
+
   // reset the coefficients
   if(fCoeff)
@@ -220 +237 @@
       // OJO!! Estimated energy will be taken directly from some input container
       Float_t emc   = fMcEvt->GetEnergy();
-      Float_t estim = fHillas->GetSize()/15.;
+      Float_t estim = fHillas->GetSize()/0.18/15.;
       Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
 
@@ -231 +248 @@
           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 energy out of bounds, skpping" << endl;
+        cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
     }
 
@@ -242 +265 @@
         const Float_t num = hmc->GetBinContent(i+1,j+1);
         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);
         if(den)
-          fCoeff->SetBinContent(i+1,j+1,num/den);
+          {
+            fCoeff->SetBinContent(i+1,j+1,num/den);
+            //borra
+            if(j==0 && den1)
+              coef1->SetBinContent(i+1,num1/den1);
+          }
         else
           {
@@ -252 +283 @@
 
 
+  //borra
+  hmc1->Write();
+  hest1->Write();
+  coef1->Write();
+
   delete hmc;
   delete hest;
@@ -270 +306 @@
   //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 = 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 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)
 
   // reset the effective areas
   if(fEffA)
     delete fEffA;
-  fEffA = new TH2F("feffa","Effective area",esbins,logemin,logemax,fNTbins,fTbin);
-
+  fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
+  //borra
+  TH1F eff("eff","Effective area",fEbins,logemin,logemax);
+  
   // fill the spectrum of the survivors
   TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin);
@@ -289 +330 @@
   // compute the effective areas
   for(Int_t j=0;j<fNTbins;j++)
-    for(Int_t i=0;i<esbins;i++)
+    for(Int_t i=0;i<fEbins;i++)
       {
-        if(horig->GetBinContent(i+1,j+1)<=0)
+        Float_t effarea =0;
+        Float_t wtot = 0;
+        for(Int_t k=0;k<fEsubbins;k++)
+          {
+            Float_t numerator = hpass->GetBinContent(i*fEsubbins+k+1,j+1);
+            Float_t denominator = horig->GetBinContent(i*fEsubbins+k+1,j+1);
+            
+            if(denominator<=0)
+              cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy subbin " << i*fEsubbins+k <<", theta bin " << j << " contains no events" << endl;
+            else
+              {
+                const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
+                const Double_t ww = fSpec->Eval(ew);
+                effarea+=Atot*numerator/denominator*ww;
+                wtot   += ww;
+              }
+          }
+        if(!wtot)
           {
             cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl;
             fEffA->SetBinContent(i+1,j+1,0);
+          }         
+        else
+          {
+            fEffA->SetBinContent(i+1,j+1,effarea/wtot); 
+            if(j==0)
+              {
+                cout << "*****" << i << ": " << effarea/wtot << endl;
+                eff.SetBinContent(i+1,effarea/wtot);
+              }
           }
-        else
-          fEffA->SetBinContent(i+1,j+1,Atot*hpass->GetBinContent(i+1,j+1)/horig->GetBinContent(i+1,j+1));
       }
-  
-
+
+  // borra
+  eff.Write();
   
   delete hpass;

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

@@ -18 +18 @@
  private:
 
-  TF1* fSpec;        // function used to parametrized the spectrum
-  TH1F* fHorig;      // histogram with the original energy spectrum
+  TF1* fSpec;        // function used to parametrize the spectrum
+  TH1F* fHorig;      // histogram with the original sample energy spectrum
 
-  Int_t fEbins;      // number of bins to build spectrum
-  Int_t fEsubbins;   // number of subbins per big bin
   Double_t fEmin;    // Minimum energy in GeV
   Double_t fEmax;    // Maximum energy in GeV
+  Int_t fEbins;      // number of bins to build spectrum
+  Int_t fEsubbins;   // number of subbins per big bin (to compute weights, eff areas...)
 
   Double_t* fWeight; // array containing weights
@@ -30 +30 @@
   TH2F* fEffA;       // histogram containing effective areas
 
-  TChain* fCini;     // chain for initial MC events (even those not triggering)
+  TChain* fCini;     // chain for initial MC files (before trigger)
   TChain* fCcut;     // chain for surviving MC events (after cuts)
 
   MHillas* fHillas;  // pointer to the MHillas Branch
   MMcEvt*  fMcEvt;   // pointer to the MMcEvt Branch
+
+  TFile* fFile; // output file (for debugging only)
 
  protected:
@@ -44 +46 @@
 
  public:
+
   MEffAreaAndCoeffCalc();
 
@@ -55 +58 @@
   void SetEmax(Float_t x)   {fEmax=x;}
 
-  void AddFileToInitialMC(const Char_t* name) {fCini->Add(name);}
-  void AddFileToFinalMC(const Char_t* name) {fCcut->Add(name);}
+  void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);}
 
   TH2F* GetEffectiveAreaHisto() {return fEffA;}

trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C

@@ -4 +4 @@
 {    
   MEffAreaAndCoeffCalc calc;
+  calc.AddFile("/data/star_gamma_test.root");
 
-  // initial MC files
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1000to1009_w0.root");
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1130to1139_w0.root");
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1520to1529_w0.root");
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1260to1269_w0.root");
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1390to1399_w0.root");
-  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1650to1659_w0.root");
-  
-
-  // files with remaining events after cuts (or where to apply cuts) see warning above
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root");
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root");
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root");
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root");
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root");
-  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root");
-
-
-  // define the funtion of the desired spectrum
-  calc.SetFunction("4.e9*pow(x,-2.6+1)",10.,10000.);
-  calc.ComputeAllFactors();
-  
   /************************************************/
   /* Build spectrum                               */
@@ -32 +11 @@
   /************************************************/
   const Int_t ebins = 10;              // number of bins to build spectrum
-  const Int_t esubbins = 20;           // number of subbins per big bin
   const Double_t emin=10;
-  const Double_t emax=1000;
+  const Double_t emax=600;
   const Double_t logemin = TMath::Log10(emin);
   const Double_t logemax = TMath::Log10(emax);
-  const Double_t de = (logemax-logemin)/ebins; // bin size (in log)
-  const Double_t desub = de/esubbins; // subbin size (in log)
-  const Int_t esbins = ebins*esubbins; // total number of subbins
+  calc.SetEbins(ebins);
+  //  calc.SetEsubbins(1);
 
-  // remaining events after cuts
-  const UInt_t ncutfiles=6;
-  Char_t* cutName[ncutfiles]={
-    "/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root",
-    "/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root",
-    "/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root",
-    "/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root",
-    "/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root",
-    "/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root"
-  };
+  // define the funtion of the desired spectrum
+  calc.SetFunction("4.e9*pow(x,-2.6+1)",emin,emax);
+  calc.ComputeAllFactors();
+  
+  const UInt_t ncutfiles=1;
+  Char_t* cutName[ncutfiles]={"/data/star_gamma_test.root"};
 
   TChain* ccut = new TChain("Events");
   for(Int_t i = 0; i < ncutfiles; i++)
     ccut->Add(cutName[i]);
-  ccut->SetAlias("logestenergy","log10(MHillas.fSize/15.)");
-  const Int_t nentries = Int_t(ccut->GetEntries());
-
-  MHillas*  hillas;
-  MMcEvt*   mcevt;
-  ccut->SetBranchStatus("*",0);
-  ccut->SetBranchStatus("MHillas.*",1);
-  ccut->SetBranchStatus("MMcEvt.*",1);
-  ccut->SetBranchAddress("MHillas.",&hillas);
-  ccut->SetBranchAddress("MMcEvt.",&mcevt);
+  ccut->SetAlias("logestenergy","log10(MHillas.fSize/0.18/15.)");
+  ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
 
   TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
   hspec->Sumw2();
+  ccut->Draw("logestenergy>>hspec","theta<10","goff");
 
-  for(Int_t i=0;i<nentries;i++)
+  for(Int_t i=0;i<ebins;i++)
     {
-      ccut->GetEntry(i);  
-      
-      // OJO!! Estimated energy will be taken directly from some input container
-      Float_t estim  = hillas->GetSize()/15.;
-
-      UInt_t effabin = UInt_t((TMath::Log10(estim)-logemin)/desub);
-      UInt_t coefbin = UInt_t((TMath::Log10(estim)-logemin)/de);
-      
-      Float_t effa  =  calc.GetEffectiveAreaHisto()->GetBinContent(effabin+1,1);
-      Float_t unfold = calc.GetCoefficientHisto()->GetBinContent(coefbin+1,1);
-      
+      const Float_t uncorrval = hspec->GetBinContent(i+1);
+      const Float_t effa  =  calc.GetEffectiveAreaHisto()->GetBinContent(i+1,1);
+      const Float_t unfold = calc.GetCoefficientHisto()->GetBinContent(i+1,1);
+      Float_t corrval;
       if(effa)
-        hspec->Fill(TMath::Log10(estim),unfold/effa*1e9);
-    }      
+        corrval = uncorrval*unfold/effa;
+      else
+        corrval = 0;
+      hspec->SetBinContent(i+1,corrval);
+    }
 
   // SAVE RESULTS