Changeset 7666 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
04/30/06 15:41:17 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mjtrain
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

    r7665 r7666  
    141141
    142142    c.cd(1);
     143    gPad->SetBorderMode(0);
     144    gPad->SetFrameBorderMode(0);
    143145    gPad->SetLogx();
    144146    h.DrawCopy();
     
    151153
    152154    c.cd(2);
     155    gPad->SetBorderMode(0);
     156    gPad->SetFrameBorderMode(0);
    153157    gPad->SetLogx();
    154158    MH::SetPalette("pretty");
     
    156160
    157161    c.cd(4);
     162    gPad->SetBorderMode(0);
     163    gPad->SetFrameBorderMode(0);
    158164    gPad->SetLogx();
    159165    MH::SetPalette("pretty");
     
    265271}
    266272
    267 Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
     273Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const
    268274{
    269275    TChain chain1("Events");
    270276    set.AddFilesOff(chain1);
    271277
     278    num = chain1.GetEntries();
     279    if (num<100)
     280    {
     281        *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
     282        return -1;
     283    }
     284
     285    TChain chain("EffectiveOnTime");
     286    set.AddFilesOff(chain);
     287
     288    chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
     289
     290    TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
     291    if (!h)
     292    {
     293        *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
     294        return -1;
     295    }
     296
     297    const Double_t ontime = h->Integral();
     298    delete h;
     299
     300    if (ontime<1)
     301    {
     302        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
     303        return -1;
     304    }
     305
     306    *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl;
     307
     308    return num/ontime;
     309}
     310
     311Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const
     312{
     313    TChain chain1("Events");
     314    set.AddFilesOn(chain1);
     315
    272316    const Double_t num = chain1.GetEntries();
    273317    if (num<100)
     
    277321    }
    278322
    279     TChain chain("EffectiveOnTime");
    280     set.AddFilesOff(chain);
    281 
    282     chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
    283 
    284     TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
    285     if (!h)
    286     {
    287         *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
    288         return -1;
    289     }
    290 
    291     const Double_t ontime = h->Integral();
    292     delete h;
    293 
    294     if (ontime<1)
    295     {
    296         *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
    297         return -1;
    298     }
    299 
    300     *fLog << inf << "Found " << num << " background events in " << ontime << "s" << endl;
    301 
    302     return num/ontime;
     323    return num;
    303324}
    304325
     
    312333
    313334    // Target spectrum
    314     TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
    315     flx.SetParameter(0, 1e-7);
     335    TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max);
     336    flx.SetParameter(0, 1e-6);
    316337
    317338    // Number n0 of events this spectrum would produce per s and m^2
     
    324345    const Double_t R = n0*A;                       //[Hz]
    325346
     347    *fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl;
     348
    326349    // Number N of events produced (in trainings sample)
    327350    const Double_t N = num;                        //[#]
    328351
     352    *fLog << "Events produced by MC inside the production area:         " << num << endl;
     353
    329354    // This correponds to an observation time T [s]
    330355    const Double_t T = N/R;                        //[s]
    331356
     357    *fLog << "Total time produced by the Monte Carlo:                   " << T << "s" << endl;
     358
    332359    // With an average data rate after star of
    333     const Double_t r = GetDataRate(fDataSetTrain); //[Hz]
     360    Double_t data=0;
     361    const Double_t r = GetDataRate(fDataSetTrain, data); //[Hz]
     362
     363    *fLog << "Events measured per second effective on time:             " << r << "Hz" << endl;
    334364
    335365    // this yields a number of n events to be read for training
    336366    const Double_t n = r*T;                        //[#]
    337367
     368    *fLog << "Events to be read from the data sample:                   " << n    << endl;
     369    *fLog << "Events available in data sample:                          " << data << endl;
     370
    338371    if (r<0)
    339372        return kFALSE;
    340373
    341     *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
    342     *fLog << "For a data rate of " << r << "Hz this corresponds to " << TMath::Nint(n) << " data events." << endl;
    343 
    344     fNumTrainOn  = (UInt_t)-1;
    345     fNumTrainOff = TMath::Nint(n);
     374    Double_t nummc = GetNumMC(fDataSetTrain);
     375
     376    *fLog << "Events available in MC sample:                            " << nummc << endl;
     377
     378    *fLog << "MC read probability:                                      " << data/n << endl;
     379
     380    // more data requested than available => Scale down num MC events
     381    Double_t on, off;
     382    if (data<n)
     383    {
     384        on  = TMath::Nint(nummc*data/n); //(UInt_t)-1;
     385        off = TMath::Nint(data);
     386        *fLog << "Not enough data events available... scaling by " << data/n << endl;
     387    }
     388    else
     389    {
     390        on  = TMath::Nint(nummc);
     391        off = TMath::Nint(n);
     392    }
     393
     394    if (fNumTrainOn>0 && fNumTrainOn<on)
     395    {
     396        fNumTrainOff = TMath::Nint(off*fNumTrainOn/on);
     397        *fLog << "Less MC events requested... scaling by " << fNumTrainOn/on << endl;
     398    }
     399    else
     400    {
     401        fNumTrainOn  = TMath::Nint(on);
     402        fNumTrainOff = TMath::Nint(off);
     403    }
     404
     405    *fLog << "Target number of MC events:   " << fNumTrainOn << endl;
     406    *fLog << "Target number of data events: " << fNumTrainOff << endl;
    346407
    347408    /*
     
    359420     MF f("f * MEventRate.fRate < rand");
    360421     */
    361 
    362422
    363423    return kTRUE;
  • trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h

    r7665 r7666  
    3030
    3131    Bool_t   GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const;
    32     Double_t GetDataRate(MDataSet &set) const;
     32    Double_t GetDataRate(MDataSet &set, Double_t &num) const;
     33    Double_t GetNumMC(MDataSet &set) const;
    3334    Bool_t   AutoTrain();
    3435
Note: See TracChangeset for help on using the changeset viewer.