Changeset 7666 for trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

Timestamp:

04/30/06 15:41:17 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

@@ -141 +141 @@
 
     c.cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
     gPad->SetLogx();
     h.DrawCopy();
@@ -151 +153 @@
 
     c.cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
     gPad->SetLogx();
     MH::SetPalette("pretty");
@@ -156 +160 @@
 
     c.cd(4);
+    gPad->SetBorderMode(0);
+    gPad->SetFrameBorderMode(0);
     gPad->SetLogx();
     MH::SetPalette("pretty");
@@ -265 +271 @@
 }
 
-Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
+Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const
 {
     TChain chain1("Events");
     set.AddFilesOff(chain1);
 
+    num = chain1.GetEntries();
+    if (num<100)
+    {
+        *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
+        return -1;
+    }
+
+    TChain chain("EffectiveOnTime");
+    set.AddFilesOff(chain);
+
+    chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
+
+    TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
+    if (!h)
+    {
+        *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
+        return -1;
+    }
+
+    const Double_t ontime = h->Integral();
+    delete h;
+
+    if (ontime<1)
+    {
+        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
+        return -1;
+    }
+
+    *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl;
+
+    return num/ontime;
+}
+
+Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const
+{
+    TChain chain1("Events");
+    set.AddFilesOn(chain1);
+
     const Double_t num = chain1.GetEntries();
     if (num<100)
@@ -277 +321 @@
     }
 
-    TChain chain("EffectiveOnTime");
-    set.AddFilesOff(chain);
-
-    chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
-
-    TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
-    if (!h)
-    {
-        *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
-        return -1;
-    }
-
-    const Double_t ontime = h->Integral();
-    delete h;
-
-    if (ontime<1)
-    {
-        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
-        return -1;
-    }
-
-    *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl;
-
-    return num/ontime;
+    return num;
 }
 
@@ -312 +333 @@
 
     // Target spectrum
-    TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
-    flx.SetParameter(0, 1e-7);
+    TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max);
+    flx.SetParameter(0, 1e-6);
 
     // Number n0 of events this spectrum would produce per s and m^2
@@ -324 +345 @@
     const Double_t R = n0*A;                       //[Hz]
 
+    *fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl;
+
     // Number N of events produced (in trainings sample)
     const Double_t N = num;                        //[#]
 
+    *fLog << "Events produced by MC inside the production area:         " << num << endl;
+
     // This correponds to an observation time T [s]
     const Double_t T = N/R;                        //[s]
 
+    *fLog << "Total time produced by the Monte Carlo:                   " << T << "s" << endl;
+
     // With an average data rate after star of
-    const Double_t r = GetDataRate(fDataSetTrain); //[Hz]
+    Double_t data=0;
+    const Double_t r = GetDataRate(fDataSetTrain, data); //[Hz]
+
+    *fLog << "Events measured per second effective on time:             " << r << "Hz" << endl;
 
     // this yields a number of n events to be read for training
     const Double_t n = r*T;                        //[#]
 
+    *fLog << "Events to be read from the data sample:                   " << n    << endl;
+    *fLog << "Events available in data sample:                          " << data << endl;
+
     if (r<0)
         return kFALSE;
 
-    *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
-    *fLog << "For a data rate of " << r << "Hz this corresponds to " << TMath::Nint(n) << " data events." << endl;
-
-    fNumTrainOn  = (UInt_t)-1;
-    fNumTrainOff = TMath::Nint(n);
+    Double_t nummc = GetNumMC(fDataSetTrain);
+
+    *fLog << "Events available in MC sample:                            " << nummc << endl;
+
+    *fLog << "MC read probability:                                      " << data/n << endl;
+
+    // more data requested than available => Scale down num MC events
+    Double_t on, off;
+    if (data<n)
+    {
+        on  = TMath::Nint(nummc*data/n); //(UInt_t)-1;
+        off = TMath::Nint(data);
+        *fLog << "Not enough data events available... scaling by " << data/n << endl;
+    }
+    else
+    {
+        on  = TMath::Nint(nummc);
+        off = TMath::Nint(n);
+    }
+
+    if (fNumTrainOn>0 && fNumTrainOn<on)
+    {
+        fNumTrainOff = TMath::Nint(off*fNumTrainOn/on);
+        *fLog << "Less MC events requested... scaling by " << fNumTrainOn/on << endl;
+    }
+    else
+    {
+        fNumTrainOn  = TMath::Nint(on);
+        fNumTrainOff = TMath::Nint(off);
+    }
+
+    *fLog << "Target number of MC events:   " << fNumTrainOn << endl;
+    *fLog << "Target number of data events: " << fNumTrainOff << endl;
 
     /*
@@ -359 +420 @@
      MF f("f * MEventRate.fRate < rand");
      */
-
 
     return kTRUE;