Ignore:
Timestamp:
02/15/05 15:36:45 (20 years ago)
Author:
rico
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mtemp/mifae/Changelog

    r6465 r6486  
    1919                                                 -*-*- END OF LINE -*-*-
    2020
    21  2005/02/10 Javier Rico
     21 2005/02/15 Javier Rico
     22    * library/MEffAreaAndCoeffCalc.[h,cc]
     23    * macros/computeCoeff.C
     24      - extend the energy bounds to compute weights
     25      - add possibility to set the zenith angle binning from outside
     26
     27 2005/02/14 Javier Rico
    2228    * library/MEffAreaAndCoeffCalc.[h,cc]
    2329    * macros/computeCoeff.C

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

    r6465 r6486  
    4747using namespace std;
    4848
    49 const Int_t fNTbins = 2;
    50 const Double_t fTbin[fNTbins+1] = {0,10,20};
    51 
     49const Int_t ntbins=1;                      // default number of theta bins
     50const Double_t tbin[ntbins+1] = {0,90};    // default theta bins bounds
     51const Int_t nebins = 10;                   // default number of energy bins
     52const Float_t emin = 10.;                  // default minimum energy value
     53const Float_t emax = 10000.;               // default maximum energy value
     54const Int_t nsubbins = 20;                 // default number of subbins per bin
     55const Char_t* deff = "4.e9*pow(x,-2.6+1)"; // default spectrum function
    5256
    5357// -------------------------------------------------------------------------
    … …  
    5660//
    5761MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
    58   : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
    59     fCoeff(NULL), fEffA(NULL), fFile(NULL)
     62  : fSpec(NULL),  fEmin(emin), fEmax(emax), fEbins(nebins), fEsubbins(nsubbins),
     63    fWeight(NULL), fCoeff(NULL), fEffA(NULL)//, fFile(NULL)
    6064{
    6165  // set the default function
    62   SetFunction("4.e9*pow(x,-2.6+1)"); 
     66  SetFunction(deff); 
    6367
    6468  // create the TChains
    65   // OJO!! make sure they read the appropriate branch
    6669  fCini = new TChain("OriginalMC");
    6770  fCcut = new TChain("Events");
    … …  
    7679  fCcut->SetBranchAddress("MHillas.",&fHillas);
    7780  fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
     81
     82  // initial value of the zenith angle binning
     83  SetThetaBinning(ntbins,tbin);
    7884 
    7985  // borra
    80   fFile = new TFile("test.root","RECREATE");
     86  //  fFile = new TFile("coeftest.root","RECREATE");
    8187}
    8288
    … …  
    9096    delete fSpec;
    9197
     98  if(fTbin)
     99    delete [] fTbin;
     100
    92101  if(fWeight)
    93102    delete [] fWeight;
    … …  
    99108    delete fEffA;
    100109
    101   if(fHorig)
    102     delete fHorig;
    103110
    104111  delete fCini;
    105112  delete fCcut;
     113}
     114
     115// -------------------------------------------------------------------------
     116//
     117// Set the binning in zenith angle
     118//
     119void MEffAreaAndCoeffCalc::SetThetaBinning(Int_t n, const Double_t* binlist)
     120{
     121  fNTbins=n;
     122  if(fTbin) delete [] fTbin;
     123  fTbin = new Double_t[n+1];
     124  for(Int_t i=0;i<n+1;i++)
     125    fTbin[i] = binlist[i];
    106126}
    107127
    … …  
    144164void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
    145165{ 
    146   const Double_t logemin = TMath::Log10(fEmin);
    147   const Double_t logemax = TMath::Log10(fEmax);
    148   const Int_t esbins = fEbins*fEsubbins; // total number of subbins
    149  
    150   if(fHorig) delete fHorig;
    151   fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
    152   fCini->Draw("logenergy>>fhorig","","goff");
    153   // borra
    154   fHorig->Write();
    155166}
    156167
    … …  
    161172void MEffAreaAndCoeffCalc::ComputeWeights()
    162173{ 
    163   const Double_t logemin = TMath::Log10(fEmin);
    164   const Double_t logemax = TMath::Log10(fEmax);
    165   const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
    166   const Double_t desub = de/fEsubbins; // subbin size (in log)
    167   const Int_t esbins = fEbins*fEsubbins; // total number of subbins
    168 
     174  // OJO!! maybe this should be hard-coded somewhere else
     175  const Float_t abslogmin=0.;
     176  const Float_t abslogmax=5.;
     177
     178  const Float_t logemin = TMath::Log10(fEmin);
     179  const Float_t logemax = TMath::Log10(fEmax);
     180  const Float_t de    = (logemax-logemin)/fEbins; // bin size (in log)
     181  const Float_t desub = de/fEsubbins; // subbin size (in log)
     182  const Int_t esbins   = fEbins*fEsubbins; // total number of subbins
     183
     184  // compute the binning for weights histogram
     185  const Int_t nmindist = (logemin>abslogmin)? Int_t((logemin-abslogmin)/desub) : 0;
     186  const Int_t nmaxdist = (logemax<abslogmax)? Int_t((abslogmax-logemax)/desub) : 0;
     187
     188  fLogEWmin = logemin-desub*nmindist;
     189  fLogEWmax = logemax+desub*nmaxdist;
     190  fEWbins   = nmindist+esbins+nmaxdist;
     191 
    169192  // reset the weights array
    170193  if(fWeight)
    171194    delete [] fWeight;
    172   fWeight = new Double_t[esbins];
     195  fWeight = new Double_t[fEWbins];
    173196   
    174   if(!fHorig)
    175     FillOriginalSpectrum();
    176  
    177   // borra
    178   TH1F hw("hw","hw",esbins,fEmin,fEmax);
    179   for(Int_t i=0;i<esbins;i++)
    180     {
    181       const Double_t denom = fHorig->GetBinContent(i+1);              // number of events in initial spectrum
    182       const Double_t ew    = TMath::Power(10,logemin+(i+0.5)*desub); // real energy
    183       const Double_t numer = fSpec->Eval(ew);                        // number of events for the required spectrum
    184       if(denom)
     197
     198  TH1F* horigs = new TH1F("horigs","Original energy spectrum",fEWbins,fLogEWmin,fLogEWmax);
     199  fCini->Draw("logenergy>>horigs","","goff");
     200  // borra
     201  //  horigs->Write();
     202
     203  // borra
     204  //  TH1F hw("hw","hw",fEWbins,fLogEWmin,fLogEWmax);
     205  for(Int_t i=0;i<fEWbins;i++)
     206    {
     207      const Float_t denom = horigs->GetBinContent(i+1);               // number of events in initial spectrum
     208      const Float_t ew    = TMath::Power(10,fLogEWmin+(i+0.5)*desub); // real energy
     209      const Float_t numer = fSpec->Eval(ew);                          // number of events for the required spectrum
     210      if(denom>10)
    185211        fWeight[i]=numer/denom;
    186212      else
    187213        {
    188           cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
     214          //      cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
    189215          fWeight[i]=-1;
    190216        }
    191217      // borra
    192       hw.SetBinContent(i+1,fWeight[i]);
    193     }
    194   // borra
    195   hw.Write();
     218      //      hw.SetBinContent(i+1,fWeight[i]);
     219    }
     220  // borra
     221  //  hw.Write();
     222
     223  delete horigs;
    196224}
    197225
    … …  
    208236    }
    209237
    210   const Double_t logemin = TMath::Log10(fEmin);
    211   const Double_t logemax = TMath::Log10(fEmax);
    212   const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
    213   const Double_t desub = de/fEsubbins; // subbin size (in log)
    214   const Int_t esbins = fEbins*fEsubbins; // total number of subbins
     238  const Float_t logemin = TMath::Log10(fEmin);
     239  const Float_t logemax = TMath::Log10(fEmax);
     240  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
     241  const Float_t desub = de/fEsubbins; // subbin size (in log)
    215242  const Int_t nentries = Int_t(fCcut->GetEntries());
    216243 
    … …  
    220247
    221248  // borra
    222   TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
    223   TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
    224   TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
     249  //  TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
     250  //  TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
     251  //  TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
    225252
    226253  // reset the coefficients
    … …  
    238265      Float_t emc   = fMcEvt->GetEnergy();
    239266      Float_t estim = fHillas->GetSize()/0.18/15.;
     267
     268      if((emc<fEmin && estim<fEmin) || (emc>fEmax && estim>fEmax) ||
     269         (emc<fEmin && estim>fEmax) || (emc>fEmax && estim<fEmin))
     270        continue;
     271
    240272      Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
    241273
    242274      // compute the bin where the weight is taken from
    243       Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub);
     275      Int_t wbin = Int_t((TMath::Log10(emc)-fLogEWmin)/desub);
    244276
    245277      // fill the histograms
    246       if(wbin<esbins)
     278      if(wbin<fEWbins)
    247279        {
    248           hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]);
    249           hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
    250           // borra
    251           if(theta<fTbin[1])
     280          if(fWeight[wbin]>0)
    252281            {
    253               hest1->Fill(TMath::Log10(estim),fWeight[wbin]);
    254               hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
     282              hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]);
     283              hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
     284              // borra
     285//            if(theta<fTbin[1])
     286//              {
     287//                hest1->Fill(TMath::Log10(estim),fWeight[wbin]);
     288//                hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
     289//              }
    255290            }
     291          else
     292            cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has no computed weight (lack of MC statistics), skipping" << endl;
    256293        }
    257294      else
    258         cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
     295        cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
    259296    }
    260297
    … …  
    266303        const Float_t den = hest->GetBinContent(i+1,j+1);
    267304        //borra
    268         const Float_t num1 = hmc1->GetBinContent(i+1);
    269         const Float_t den1 = hest1->GetBinContent(i+1);
     305//      const Float_t num1 = hmc1->GetBinContent(i+1);
     306//      const Float_t den1 = hest1->GetBinContent(i+1);
    270307        if(den)
    271308          {
    272309            fCoeff->SetBinContent(i+1,j+1,num/den);
    273310            //borra
    274             if(j==0 && den1)
    275               coef1->SetBinContent(i+1,num1/den1);
     311//          if(j==0 && den1)
     312//            coef1->SetBinContent(i+1,num1/den1);
    276313          }
    277314        else
    … …  
    284321
    285322  //borra
    286   hmc1->Write();
    287   hest1->Write();
    288   coef1->Write();
     323//   hmc1->Write();
     324//   hest1->Write();
     325//   coef1->Write();
    289326
    290327  delete hmc;
    … …  
    305342
    306343  //OJO!! radius should be read from somewhere!
    307   const Double_t radius = 30000.; // generation radius (in cm)
    308   //  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
    309   const Double_t Atot = 1.; //total area (in cm^2)
    310   const Double_t logemin = TMath::Log10(fEmin);
    311   const Double_t logemax = TMath::Log10(fEmax);
     344  const Float_t radius = 30000.; // generation radius (in cm)
     345  const Float_t Atot = 3.14159*radius*radius; //total area (in cm^2)
     346  const Float_t logemin = TMath::Log10(fEmin);
     347  const Float_t logemax = TMath::Log10(fEmax);
    312348  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
    313   const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
    314   const Double_t desub = de/fEsubbins; // subbin size (in log)
     349  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
     350  const Float_t desub = de/fEsubbins; // subbin size (in log)
    315351
    316352  // reset the effective areas
    … …  
    319355  fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
    320356  //borra
    321   TH1F eff("eff","Effective area",fEbins,logemin,logemax);
     357//   TH1F eff("eff","Effective area",fEbins,logemin,logemax);
    322358 
    323359  // fill the spectrum of the survivors
    … …  
    343379            else
    344380              {
    345                 const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
    346                 const Double_t ww = fSpec->Eval(ew);
     381                const Float_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
     382                const Float_t ww = fSpec->Eval(ew);
    347383                effarea+=Atot*numerator/denominator*ww;
    348384                wtot   += ww;
    … …  
    357393          {
    358394            fEffA->SetBinContent(i+1,j+1,effarea/wtot);
    359             if(j==0)
    360               {
    361                 cout << "*****" << i << ": " << effarea/wtot << endl;
    362                 eff.SetBinContent(i+1,effarea/wtot);
    363               }
     395            // borra
     396//          if(j==0)
     397//            {
     398//              //              cout << "*****" << i << ": " << effarea/wtot << endl;
     399//              eff.SetBinContent(i+1,effarea/wtot);
     400//            }
    364401          }
    365402      }
    366403
    367404  // borra
    368   eff.Write();
     405//   eff.Write();
    369406 
    370407  delete hpass;

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

    r6465 r6486  
    1919
    2020  TF1* fSpec;        // function used to parametrize the spectrum
    21   TH1F* fHorig;      // histogram with the original sample energy spectrum
    2221
    23   Double_t fEmin;    // Minimum energy in GeV
    24   Double_t fEmax;    // Maximum energy in GeV
     22  Float_t fEmin;    // Minimum energy in GeV
     23  Float_t fEmax;    // Maximum energy in GeV
    2524  Int_t fEbins;      // number of bins to build spectrum
    2625  Int_t fEsubbins;   // number of subbins per big bin (to compute weights, eff areas...)
     26
     27  Float_t fLogEWmin; // Log Minimum energy for weights (in GeV)
     28  Float_t fLogEWmax; // Log Maximum energy for weights (in GeV)
     29  Int_t    fEWbins;   // Number of bins for weights
     30
     31  Int_t fNTbins;      // Number of bins in zenith angle
     32  Double_t* fTbin;     // array containing bin boundaries (size must be fNTbins+)
    2733
    2834  Double_t* fWeight; // array containing weights
    … …  
    5864  void SetEmax(Float_t x)   {fEmax=x;}
    5965
     66  void SetThetaBinning(Int_t n, const Double_t* binlist);
     67
    6068  void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);}
    6169

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

    r6465 r6486  
    1010  /* now it's just a cross check with the same MC */
    1111  /************************************************/
    12   const Int_t ebins = 10;              // number of bins to build spectrum
     12  const Int_t ebins = 20;              // number of bins to build spectrum
    1313  const Double_t emin=10;
    14   const Double_t emax=600;
    15   const Double_t logemin = TMath::Log10(emin);
    16   const Double_t logemax = TMath::Log10(emax);
     14  const Double_t emax=200;
     15  const Int_t tbins = 2;
     16  const Double_t thetab[tbins+1] = {0,10,20};
    1717  calc.SetEbins(ebins);
    18   //  calc.SetEsubbins(1);
    19 
     18  calc.SetEmin(emin);
     19  calc.SetEmax(emax);
     20  calc.SetThetaBinning(tbins,thetab);
     21 
    2022  // define the funtion of the desired spectrum
    21   calc.SetFunction("4.e9*pow(x,-2.6+1)",emin,emax);
     23  calc.SetFunction("4.e9*pow(x,-2.6+1)");
    2224  calc.ComputeAllFactors();
    2325 
    … …  
    3133  ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
    3234
     35
     36  const Double_t logemin = TMath::Log10(emin);
     37  const Double_t logemax = TMath::Log10(emax);
    3338  TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
    3439  hspec->Sumw2();
    … …  
    4247      Float_t corrval;
    4348      if(effa)
    44         corrval = uncorrval*unfold/effa;
     49        corrval = uncorrval*unfold/effa*1e9;
    4550      else
    4651        corrval = 0;
Note: See TracChangeset for help on using the changeset viewer.