Ignore:
Timestamp:
07/28/04 11:01:15 (20 years ago)
Author:
jlopez
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r4405 r4429  
    3939  Double_t GetWidthRealTime() { return fWidthRealTimeBin[fNumberTimeBins-1]; }
    4040
    41   Double_t* GetLiveTimeArray()   { return fLiveTimeBin.GetArray(); }
    42   Double_t* GetMeanRealTimeArray()   { return fMeanRealTimeBin.GetArray(); }
    43   Double_t* GetWidthRealTimeArray()   { return fWidthRealTimeBin.GetArray(); }
     41  TArrayD& GetLiveTimeTArray()      { return fLiveTimeBin; }
     42  TArrayD& GetMeanRealTimeTArray()  { return fMeanRealTimeBin; }
     43  TArrayD& GetWidthRealTimeTArray() { return fWidthRealTimeBin; }
    4444
    4545  void Print(const Option_t*) const;

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

    r4408 r4429  
    4343
    4444Bool_t Debug = kFALSE;
    45 Bool_t PrintNewRun = kFALSE;
     45Bool_t PrintNewRun = kTRUE;
    4646
    4747MLiveTimeCalc::MLiveTimeCalc(const char *name, const char *title) : kSecTomSec(1e3), kDayToSec(24.*60.*60.), fRunHeader(NULL), fEvtHeader(NULL), fPresentEventTime(NULL), fLiveTime(NULL)

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

    r3947 r4429  
    2525}
    2626
    27 //void alpha_plot(TString f_on_name  = "../HillasFiles/20040319_Mrk421_30_15.KDummy*.root",
    28 void alpha_plot(TString f_on_name  = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root",
    29                 TString f_off_name = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root",
     27void alpha_plot(TString f_on_name  = "../HillasFiles/Mrk421/*_H.root",
     28                TString f_off_name = "../HillasFiles/OffMrk421/*_H.root",
    3029                TString f_src_name = "../HillasFiles/20040319_OffMrk421.fake.pos")           
    31 /*void alpha_plot(TString f_on_name  = "../HillasFiles/mrk421OnMisp3015*.root",
    32                 TString f_off_name = "../HillasFiles/mrk421OffMisp3015*.root",
    33                 TString f_src_name = "../HillasFiles/20040215_Mrk421.pos")
    34 */
    3530{
    3631
    … …  
    3833   
    3934    //cuts
    40     Float_t sizemin   = 500.; //[ADC]
     35    Float_t sizemin   = 2000.; //[ADC]
    4136    Float_t sizemax   = 10000000000.; //[ADC]
    42     Float_t widthmin  = 0.;
    43     Float_t widthmax  = 0.8;
    44     Float_t lengthmin = 0.;
    45     Float_t lengthmax = 0.8;
    46     Float_t distmin   = 0.;
    47     Float_t distmax   = 2.;
     37    Float_t widthmin  = 0.06;
     38    Float_t widthmax  = 0.12;
     39    Float_t lengthmin = 0.10;
     40    Float_t lengthmax = 0.26;
     41    Float_t distmin   = 0.3;
     42    Float_t distmax   = 1.2;
    4843    Float_t alphamin   = 0.;
    4944    Float_t alphamax   = 90.;
    … …  
    5247    Float_t sigexccmin = 0.;
    5348    Float_t sigexccmax = 30.;
    54     Float_t bkgnormmin = 0.;
     49    Float_t bkgnormmin = 30.;
    5550    Float_t bkgnormmax = 90.;
    5651   
    … …  
    173168
    174169    //cuts
    175     TString sizestr = "MHillas.fSize < ";
     170    TString sizestr = "(MHillas.fSize < ";
    176171    sizestr += sizemin;
    177     sizestr += " || ";
     172    sizestr += ") || (";
    178173    sizestr += "MHillas.fSize > ";
    179174    sizestr += sizemax;
     175    sizestr += ")";
    180176    MF sizefilter(sizestr);
    181177   
    182     TString widthstr = "{MHillas.fWidth/315.} < ";
     178    TString widthstr = "({MHillas.fWidth/315.} < ";
    183179    widthstr += widthmin;
    184     widthstr += " || ";
     180    widthstr += ") || (";
    185181    widthstr += "{MHillas.fWidth/315.} > ";
    186182    widthstr += widthmax;
     183    widthstr += ")";
    187184    MF widthfilter(widthstr);
    188185   
    189     TString lengthstr = "{MHillas.fLength/315.} < ";
     186    TString lengthstr = "({MHillas.fLength/315.} < ";
    190187    lengthstr += lengthmin;
    191     lengthstr += " || ";
     188    lengthstr += ") || (";
    192189    lengthstr += "{MHillas.fLength/315.} > ";
    193190    lengthstr += lengthmax;
     191    lengthstr += ")";
    194192    MF lengthfilter(lengthstr);
    195193   
    196     TString diststr = "{MHillasSrc.fDist/315.} < ";
     194    TString diststr = "({MHillasSrc.fDist/315.} < ";
    197195    diststr += distmin;
    198     diststr += " || ";
     196    diststr += ") || (";
    199197    diststr += "{MHillasSrc.fDist/315.} > ";
    200198    diststr += distmax;
     199    diststr += ")";
    201200    MF distfilter(diststr);
    202201   
    203     TString alphastr = "{abs(MHillasSrc.fAlpha)} < ";
     202    TString alphastr = "({abs(MHillasSrc.fAlpha)} < ";
    204203    alphastr += alphamin;
    205     alphastr += " || ";
     204    alphastr += ") || (";
    206205    alphastr += "{abs(MHillasSrc.fAlpha)} > ";
    207206    alphastr += alphamax;
     207    alphastr += ")";
    208208    MF alphafilter(alphastr);
    209209   
    … …  
    219219    MContinue cont_odd(&oddfilter);
    220220   
    221     MSrcPosFromFile srccalc(f_src_name);
     221    //    MSrcPosFromFile srccalc(f_src_name);
     222    MSrcPlace srccalc;
    222223   
    223224    MHillasSrcCalc csrc_on;
    … …  
    245246    tlist_on.AddToList(&csrc_on);
    246247    tlist_on.AddToList(&fsrcpos_on);
    247     tlist_on.AddToList(&cont_odd);
     248//    tlist_on.AddToList(&cont_odd);
    248249    tlist_on.AddToList(&cont_size);
    249250    tlist_on.AddToList(&cont_width);
    … …  
    357358    read_off.DisableAutoScheme();
    358359   
    359     srccalc.SetMode(MSrcPosFromFile::kOff);
     360    srccalc.SetMode(MSrcPlace::kOff);
    360361
    361362    MHillasSrcCalc csrc_off;
    … …  
    376377    tlist_off.AddToList(&csrc_off);
    377378    tlist_off.AddToList(&fsrcpos_off);
    378     tlist_off.AddToList(&cont_even);
     379//    tlist_off.AddToList(&cont_even);
    379380    tlist_off.AddToList(&cont_size);
    380381    tlist_off.AddToList(&cont_width);

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

    r4294 r4429  
    4747
    4848
    49 void findstars(const TString filename="dc_2004_03_17_01_16_51_20440_Mrk421.root", const TString directory="/nfs/magic/CaCodata/rootdata/Mrk421/Period015/2004_03_17/", const UInt_t numEvents = 0)
     49void findstars(const TString filename="20040422_*_D_Mrk421_E.root", const TString directory="/local_disk/Data/rootdata/Mrk421/Period016/2004_04_22/", const UInt_t numEvents = 0)
    5050{
    5151
    … …  
    8181  MGeomApply geomapl;
    8282  TString continuoslightfile =
    83     //    "/home/Javi/mnt_magic_data/CaCo/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
    84      "/nfs/magic/CaCodata/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
     83    "/local_disk/CaCoData/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root";
    8584
    8685  Float_t mindc = 0.9; //[uA]
    … …  
    9392  const TArrayS blindpixels(numblind,(Short_t*)x);
    9493  Float_t ringinterest = 100; //[mm]
    95   Float_t tailcut = 2.5;
    96   UInt_t integratedevents = 1;
     94  Float_t tailcut = 3.5;
     95  UInt_t integratedevents = 4;
    9796
    9897  MFindStars findstars;

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

    r4035 r4429  
    2727       
    2828    // define file, tree, canvas
    29     TFile *fdata = new TFile(onname);
     29    TChain *fdata = new TChain;
     30    fdata->AddFile(onname);
    3031    TTree *t = Parameters;                                                     
    3132   
    3233    // define file, tree, canvas
    33     TFile *fdata2 = new TFile(offname);
     34    TChain *fdata2 = new TChain;
     35    fdata2->AddFile(offname);
    3436    TTree *t2 = Parameters; 
    3537   

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

    r3977 r4429  
    4747
    4848
    49 void hvnotnominal(const TString filename="20040319_20821_D_Mrk421_S.root", const TString directory="/nfs/magic/CaCodata/2004_03_19/")
     49void hvnotnominal(const TString filename="20040319_20821_D_Mrk421_S.root", const TString directory="/nfs/magic/CaCodata/2004_03_19/", Float_t percent = 0.01)
    5050{
    5151
    … …  
    8686  MGeomApply geomapl;
    8787  MFHVNotNominal fHVNominal;
    88   TString hvnominal = "/nfs/ifae.es/user/j/jlopez/Mars/ConfigFiles/HVSettings_FF35q.conf";
     88  TString hvnominal = "/mnt/users/jlopez/Mars/Files4Mars/Config/HVSettings_FF35q.conf";
    8989  fHVNominal.SetHVNominalValues(hvnominal);
    9090  //  fHVNominal.SetMaxNumPixelsDeviated(10);
    … …  
    105105  //
    106106
    107   if (!evtloop.Eventloop())
     107  MHCamera hvnotnominal(geomcam);
     108  MHCamera lowhvnotnominal(geomcam);
     109  MHCamera uphvnotnominal(geomcam);
     110
     111  if (!evtloop.PreProcess())
    108112    return;
    109113
     114  TArrayD nominal = fHVNominal.GetHVNominal();
     115  TArrayD lownominal = nominal;
     116  TArrayD upnominal = nominal;
     117 
     118  for (UInt_t pix=0; pix<nominal.GetSize(); pix++)
     119    {
     120      lownominal[pix] *= (1-percent);
     121      upnominal[pix]  *= (1+percent);
     122    }
     123 
     124 
     125  while(tlist.Process())
     126    {
     127     
     128      hvnotnominal.CntCamContent(hvcam,lownominal,0,kFALSE);
     129      hvnotnominal.CntCamContent(hvcam,upnominal,0,kTRUE);
     130     
     131      lowhvnotnominal.CntCamContent(hvcam,lownominal,0,kFALSE);
     132      uphvnotnominal.CntCamContent(hvcam,upnominal,0,kTRUE);
     133    }
     134 
     135
     136  evtloop.PostProcess();
     137 
    110138  tlist.PrintStatistics();
     139 
     140  TCanvas c1;
     141  c1.Divide(2,1);
     142  c1.cd(1);
     143  lowhvnotnominal.SetPrettyPalette();
     144  lowhvnotnominal.Draw();
     145  gPad->cd(1);
     146  gPad->Modified();
     147  gPad->Update();
     148  c1.cd(2);
     149  uphvnotnominal.SetPrettyPalette();
     150  uphvnotnominal.Draw();
     151  gPad->cd(2);
     152  gPad->Modified();
     153  gPad->Update();
     154
     155  TCanvas c2;
     156  hvnotnominal.SetPrettyPalette();
     157  hvnotnominal.Draw();
     158  gPad->cd(1);
     159  gPad->Modified();
     160  gPad->Update();
     161
     162  if (!HandleInput())
     163    break;
    111164
    112165}

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

    r4408 r4429  
    1 Double_t ChiSquareNDof(TH1D *h1, TH1D *h2)
    2 {
    3     Double_t chiq = 0.;
    4     Double_t chi;
    5     Double_t error;
    6     Int_t nbinsnozero = 0;
    7 
    8     Int_t nbins = h1->GetNbinsX();
    9     if (nbins != h2->GetNbinsX() || nbins == 0)
    10         return -1;
    11 
    12     for (UInt_t bin=1; bin<=nbins; bin++)
    13     {
    14         error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) +
    15                            h2->GetBinError(bin)*h2->GetBinError(bin));
    16         if (error != 0)
    17         {
    18             chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error;
    19             chiq += chi*chi;
    20             nbinsnozero++;
    21         }
    22     }
    23 
    24     return (nbinsnozero>0?chiq/nbinsnozero:0);
    25 }
    26 
    27 Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[])
    28 {
    29 
    30   Int_t result = -1;
    31 
    32   Int_t lowerbin = 0;
    33   Int_t upperbin = numberSizeBins;
    34   Int_t bin;
    35 
    36   Int_t count = 0;
    37 
    38   if (size >= sizeBins[0])
    39     {
    40       while (upperbin - lowerbin > 1 && count++<=numberSizeBins)
    41         {
    42           bin = (upperbin+lowerbin)/2;
    43           if (size >= sizeBins[bin])
    44             lowerbin = bin;
    45           else
    46             upperbin = bin;
    47         }
    48       result = count<=numberSizeBins?lowerbin:-1;
    49     }
    50 
    51   return result;
    52 
    53 }
    54 
     1
     2Double_t ChiSquareNDof(TH1D *h1, TH1D *h2);
     3Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[]);
    554
    565void lightcurve(TString f_on_name = "../HillasFiles/Mrk421/*_H.root",
    … …  
    6716  // Configuration
    6817  const Bool_t debug = kFALSE;
    69   const Double_t timebin = 600.; //[sec]
    70   TString psname = "./20040420_Mrk421.600s.ps";
     18  const Double_t timebin = 1380.; //[sec]
     19  TString psname = "./20040421_Mrk421.1380s.ps";
    7120
    7221  //Constanst
    … …  
    9443  TObjArray timediffHistoOn;
    9544 
    96   const Int_t numberSizeBins = 4;
    97   Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]
     45  const Int_t numberSizeBins = 3;
     46  Double_t sizeBins[numberSizeBins] = {1897., 2785., 4087.}; //[Photons]
    9847
    9948  //cuts
    10049
    101   Double_t widthmin[numberSizeBins]  = { 0.06, 0.06, 0.06, 0.06 };
    102   Double_t widthmax[numberSizeBins]  = { 0.10, 0.10, 0.12, 0.12 };
    103   Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 };
    104   Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
    105   Double_t distmin[numberSizeBins]   = { 0.30, 0.30, 0.30, 0.30 };
    106   Double_t distmax[numberSizeBins]   = { 1.20, 1.20, 1.20, 1.20 };
     50  Double_t widthmin[numberSizeBins]  = { 0.06, 0.06, 0.06 };
     51  Double_t widthmax[numberSizeBins]  = { 0.10, 0.12, 0.12 };
     52  Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10 };
     53  Double_t lengthmax[numberSizeBins] = { 0.25, 0.26, 0.36 };
     54  Double_t distmin[numberSizeBins]   = { 0.30, 0.30, 0.30 };
     55  Double_t distmax[numberSizeBins]   = { 1.20, 1.20, 1.20 };
     56
     57//   const Int_t numberSizeBins = 4;
     58//   Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]
     59
     60//   //cuts
     61
     62//   Double_t widthmin[numberSizeBins]  = { 0.06, 0.06, 0.06, 0.06 };
     63//   Double_t widthmax[numberSizeBins]  = { 0.10, 0.10, 0.12, 0.12 };
     64//   Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 };
     65//   Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
     66//   Double_t distmin[numberSizeBins]   = { 0.30, 0.30, 0.30, 0.30 };
     67//   Double_t distmax[numberSizeBins]   = { 1.20, 1.20, 1.20, 1.20 };
    10768
    10869  //alpha plot integration for gammas
    10970  Double_t sigexccmin = 0.;
    110   Double_t sigexccmax = 20.;
     71  Double_t sigexccmax = 15.;
    11172  Double_t bkgnormmin = 40.;
    11273  Double_t bkgnormmax = 80.;
    … …  
    12889  MRawEvtHeader evtheader_on;
    12990  MTime              time_on;
     91  MLiveTime      livetime_on;
    13092  MHillas          hillas_on;
    13193  MHillasSrc    hillassrc_on;
    … …  
    13799  plist_on.AddToList(&evtheader_on);
    138100  plist_on.AddToList(&time_on);
     101  plist_on.AddToList(&livetime_on);
    139102  plist_on.AddToList(&hillas_on);
    140103  plist_on.AddToList(&hillassrc_on);
    … …  
    186149  MSrcPlace srcplace;
    187150  MHillasSrcCalc csrc_on;
    188  
     151
     152  MLiveTimeCalc livetimecalc_on;
     153  livetimecalc_on.SetRealTimeBinSize(timebin);
     154
    189155//   // prints
    190156//   MPrint pevent("MRawEvtHeader");
    … …  
    195161  //tasklist
    196162  tlist_on.AddToList(&read_on);
     163  tlist_on.AddToList(&livetimecalc_on);
    197164  tlist_on.AddToList(&srcplace);
    198165  tlist_on.AddToList(&csrc_on);
    … …  
    212179    return;
    213180 
    214   Double_t mjdFirstEventofBin=0;
    215   Double_t mjdFirstEvent=0;
    216 
    217   MTime mtimeFirstEventofBin;
    218   MTime mtimeFirstEvent;
    219 
    220   Double_t mjdLastEvent=0;
    221   MTime  mtimeLastEvent=0;
    222   UInt_t   evtLastEvent;
    223   UInt_t   runLastEvent=0;
    224   MTime    mtimeLastEvent;
    225 
    226   Bool_t flag = kFALSE;
    227 
    228181  numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
    229182  zenithMinimumOn[numberTimeBins-1] = 100.;
    … …  
    252205  TF1 *f1 = new TF1("exp","expo",mintimediff,maxtimediff);           
    253206
     207  Double_t lastSrcPosX = 0;
     208  Double_t lastSrcPosY = 0;
     209
    254210  while(tlist_on.Process())
    255211    {
    256       // Compute live time
    257       UInt_t   run = runheader_on.GetRunNumber();
    258       UInt_t   evt = evtheader_on.GetDAQEvtNumber();
    259       Double_t mjd = time_on.GetMjd();
     212      numberOnEventsAfterCleaning[numberTimeBins-1]++;
    260213     
    261       numberOnEventsAfterCleaning[numberTimeBins-1]++;
    262                                  
    263       if (mjd == 0)
    264       {
    265        
    266           if (debug)
    267             {
    268               cout << "MTime::GetTime() == 0" << endl;
    269               cout.precision(15);
    270               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    271               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    272               mtimeLastEvent.Print();
    273               time_on.Print();
    274             }
    275          
    276           flag = kTRUE;
    277        }
    278       else if (mjd-mjdLastEvent < 0 && flag)
    279         {
    280          
    281           if (debug)
    282             {
    283               cout << "mjd-mjdLastEvent < 0" << endl;
    284               cout.precision(15);
    285               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    286               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    287               mtimeLastEvent.Print();
    288               time_on.Print(); 
    289             }
    290          
    291           flag = kTRUE;
    292        }
    293       else if (mjd-mjdLastEvent == 0 && flag)
    294         {
    295          
    296           if (debug)
    297             {
    298               cout << "mjd-mjdLastEvent == 0" << endl;
    299               cout.precision(15);
    300               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    301               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    302               mtimeLastEvent.Print();
    303               time_on.Print(); 
    304             }
    305 
    306           flag = kTRUE;
    307        }
    308       else if (evt-evtLastEvent <= 0)
    309         {
    310          
    311           if (debug)
    312             {
    313               cout << "evt-evtLastEvent <= 0" << endl;
    314               cout.precision(15);
    315               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    316               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    317               mtimeLastEvent.Print();
    318               time_on.Print(); 
    319             }
    320 
    321           flag = kTRUE;
    322         }
    323 
    324       else if ((Int_t)mjd == mjd)
    325         {
    326          
    327           if (debug)
    328             {
    329               cout << "(Int_t)mjd == mjd" << endl;
    330               cout.precision(15);
    331               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    332               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    333               mtimeLastEvent.Print();
    334               time_on.Print(); 
    335             }
    336 
    337           flag = kTRUE;
    338         }
     214      if (srcpos_on.GetX() == 0. && srcpos_on.GetY() == 0.)
     215        srcpos_on.SetXY(lastSrcPosX,lastSrcPosY);
    339216      else
    340217        {
    341 
    342           if (flag && debug)
     218          lastSrcPosX = srcpos_on.GetX();
     219          lastSrcPosY = srcpos_on.GetY();
     220        }
     221      srcplace.CallProcess();
     222      csrc_on.CallProcess();
     223
     224     
     225      Double_t size = hillas_on.GetSize();
     226      Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg();
     227      Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg();
     228      Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg();
     229      Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg();
     230      Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany);
     231      Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg();
     232      Double_t alpha = hillassrc_on.GetAlpha();
     233      Double_t srcposx = srcpos_on.GetX();
     234      Double_t srcposy = srcpos_on.GetY();
     235      Double_t zenith = drive_on.GetNominalZd();
     236         
     237
     238      Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
     239         
     240      if (sizebin >= 0)
     241        {
     242          if (width > widthmin[sizebin] && width < widthmax[sizebin])
    343243            {
    344            
    345               cout << "flag = TRUE" << endl;
    346               cout.precision(15);
    347               cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    348               cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    349               mtimeLastEvent.Print();
    350               time_on.Print(); 
    351              
    352               flag = kFALSE;
    353 
    354             }
    355 
    356           if ( run !=  runLastEvent )
    357             {
    358               if ( mjdLastEvent != 0 )
     244              if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
    359245                {
    360                   timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
    361                  
    362                   cout.precision(10);
    363                   cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
    364                   if (flag && debug)
    365                     {
    366                       cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
    367                       cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
    368                       mtimeLastEvent.Print();
    369                       time_on.Print(); 
    370                     }
    371                  
    372                 }
    373              
    374               mjdFirstEvent = mjd;
    375             }
    376           else if (mjdFirstEventofBin != 0)
    377             hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec);
    378          
    379           if (mjdFirstEventofBin == 0)
    380             {
    381               mjdFirstEvent = mjd;
    382               mjdFirstEventofBin = mjd;
    383             }
    384 
    385           evtLastEvent = evt;
    386           runLastEvent = run;
    387           mjdLastEvent = mjd;
    388           mtimeLastEvent = time_on;
    389      
    390           Double_t size = hillas_on.GetSize();
    391           Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg();
    392           Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg();
    393           Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg();
    394           Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg();
    395           Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany);
    396           Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg();
    397           Double_t alpha = hillassrc_on.GetAlpha();
    398           Double_t srcposx = srcpos_on.GetX();
    399           Double_t srcposy = srcpos_on.GetY();
    400           Double_t zenith = drive_on.GetNominalZd();
    401          
    402 
    403           Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
    404          
    405           if (sizebin >= 0)
    406             {
    407               if (width > widthmin[sizebin] && width < widthmax[sizebin])
    408                 {
    409                   if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
    410                     {
    411                       if (dist > distmin[sizebin] && centerdist < distmax[sizebin])
    412                         {     
    413 
    414                           //General histos
    415                           hSize_on->Fill(log10(size));
    416                           hWidth_on->Fill(width);
    417                           hLength_on->Fill(length);
    418                           hDist_on->Fill(dist);
    419                           hAlpha_on_abs->Fill(TMath::Abs(alpha));
    420                           hAlpha_on->Fill(alpha);
    421                           hSrcPos_on->Fill(srcposx,srcposy);
    422                          
    423                           // Time bin histos
    424                           hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
    425                           hSrcPos_on_timebin->Fill(srcposx,srcposy);
    426                           hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad));
    427                          
    428                           if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1])
    429                             zenithMinimumOn[numberTimeBins-1] = zenith;
    430                           if (zenith>zenithMaximumOn[numberTimeBins-1])
    431                             zenithMaximumOn[numberTimeBins-1] = zenith;
    432                          
    433                         }
     246                  if (dist > distmin[sizebin] && centerdist < distmax[sizebin])
     247                    {     
     248                     
     249                      //General histos
     250                      hSize_on->Fill(log10(size));
     251                      hWidth_on->Fill(width);
     252                      hLength_on->Fill(length);
     253                      hDist_on->Fill(dist);
     254                      hAlpha_on_abs->Fill(TMath::Abs(alpha));
     255                      hAlpha_on->Fill(alpha);
     256                      hSrcPos_on->Fill(srcposx,srcposy);
     257                     
     258                      // Time bin histos
     259                      hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
     260                      hSrcPos_on_timebin->Fill(srcposx,srcposy);
     261                      hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad));
     262                     
     263                      if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1])
     264                        zenithMinimumOn[numberTimeBins-1] = zenith;
     265                      if (zenith>zenithMaximumOn[numberTimeBins-1])
     266                        zenithMaximumOn[numberTimeBins-1] = zenith;
     267                     
    434268                    }
    435269                }
    436270            }
     271        }
     272     
     273     
     274      // Check if you are overload the maxim time bin
     275      if (numberTimeBins != livetime_on.GetNumberTimeBins())
     276        {
     277          timeOn[numberTimeBins-1] = livetime_on.GetLiveTimeTArray().At(numberTimeBins-1);
     278          meanTimeBinOn[numberTimeBins-1] = livetime_on.GetMeanRealTimeTArray().At(numberTimeBins-1);
     279          widthTimeBinOn[numberTimeBins-1] = livetime_on.GetWidthRealTimeTArray().At(numberTimeBins-1);
     280
     281          //Compute the number of on events
     282          numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90);
     283          numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90);         
    437284         
    438  
    439       // Check if you are overload the maxim time bin
    440         if ((mjd-mjdFirstEventofBin)*kDay >= timebin)
    441           {
    442               //Compute the time on
    443               timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
    444               widthTimeBinOn[numberTimeBins-1] = (mjdLastEvent-mjdFirstEventofBin)/2;
    445               meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
    446              
    447               cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
    448               cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " widthTimeBinOn " <<  widthTimeBinOn[numberTimeBins-1] << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << ' ' <<  mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1] << endl;
    449              
    450               //Compute the number of on events
    451               numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
    452               numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);           
    453               //Initialize variables
    454               mjdFirstEvent = mjd;
    455               mjdFirstEventofBin = mjd;
    456 
    457               hTimeDiff_on_timebin->Fit("exp","RQ0");
    458               deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
    459               cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" <<  f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl;
    460               cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1)))  << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
    461 
    462               alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
    463               srcposHistoOn.AddLast(hSrcPos_on_timebin);
    464               coszenithHistoOn.AddLast(hCosZenith_on_timebin);
    465               timediffHistoOn.AddLast(hTimeDiff_on_timebin);
    466              
    467 //            cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
    468 
    469               //Increase the number of time bins and all needed arrays
    470               numberTimeBins++;
    471              
    472               timeOn.Set(numberTimeBins);
    473               numberOnEvents.Set(numberTimeBins);
    474               numberBkgEventsToNormOn.Set(numberTimeBins);
    475               widthTimeBinOn.Set(numberTimeBins);
    476               meanTimeBinOn.Set(numberTimeBins);
    477               zenithMinimumOn.Set(numberTimeBins);
    478               zenithMaximumOn.Set(numberTimeBins);
    479               deadFractionOn.Set(numberTimeBins);
    480               numberOnEventsAfterCleaning.Set(numberTimeBins);
    481 
    482               timeOn[numberTimeBins-1] = 0.;
    483               zenithMinimumOn[numberTimeBins-1] = 100.;
    484               zenithMaximumOn[numberTimeBins-1] = 0.;
    485               numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
    486              
    487               alphaTitle =  Form("%s%02i","hAlphaOn",numberTimeBins-1);         
    488               hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
    489              
    490               srcposTitle =  Form("%s%02i","hSrcPosOn",numberTimeBins-1);
    491               hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
    492 
    493               coszenithTitle =  Form("%s%02i","hCosZenithOn",numberTimeBins-1);
    494               hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
    495 
    496               timediffTitle =  Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
    497               hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
    498 
    499               flag = kTRUE;
    500           }
     285//        hTimeDiff_on_timebin->Fit("exp","RQ0");
     286//        deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
     287//        cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" <<  f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl;
     288//        cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1)))  << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
     289          deadFractionOn[numberTimeBins-1] = 1.;
     290         
     291          alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
     292          srcposHistoOn.AddLast(hSrcPos_on_timebin);
     293          coszenithHistoOn.AddLast(hCosZenith_on_timebin);
     294          timediffHistoOn.AddLast(hTimeDiff_on_timebin);
     295         
     296          numberTimeBins = livetime_on.GetNumberTimeBins();
     297         
     298          timeOn.Set(numberTimeBins);
     299          numberOnEvents.Set(numberTimeBins);
     300          numberBkgEventsToNormOn.Set(numberTimeBins);
     301          widthTimeBinOn.Set(numberTimeBins);
     302          meanTimeBinOn.Set(numberTimeBins);
     303          zenithMinimumOn.Set(numberTimeBins);
     304          zenithMaximumOn.Set(numberTimeBins);
     305          deadFractionOn.Set(numberTimeBins);
     306          numberOnEventsAfterCleaning.Set(numberTimeBins);
     307         
     308          timeOn[numberTimeBins-1] = 0.;
     309          zenithMinimumOn[numberTimeBins-1] = 100.;
     310          zenithMaximumOn[numberTimeBins-1] = 0.;
     311          numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
     312         
     313          alphaTitle =  Form("%s%02i","hAlphaOn",numberTimeBins-1);         
     314          hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
     315         
     316          srcposTitle =  Form("%s%02i","hSrcPosOn",numberTimeBins-1);
     317          hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
     318         
     319          coszenithTitle =  Form("%s%02i","hCosZenithOn",numberTimeBins-1);
     320          hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
     321         
     322          timediffTitle =  Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
     323          hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
     324         
     325        }
     326     
     327    }
    501328       
    502         }
    503        
    504     }
    505 
    506  
    507   //Compute the time on for last time bin
    508   timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
    509   widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2;
    510   meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
     329  loop_on.PostProcess();
     330 
     331  tlist_on.PrintStatistics();
     332 
     333  timeOn[numberTimeBins-1] = livetime_on.GetLiveTimeTArray().At(numberTimeBins-1);
     334  meanTimeBinOn[numberTimeBins-1] = livetime_on.GetMeanRealTimeTArray().At(numberTimeBins-1);
     335  widthTimeBinOn[numberTimeBins-1] = livetime_on.GetWidthRealTimeTArray().At(numberTimeBins-1);
     336 
    511337  //Compute the number of on events for the last time bin
    512   numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
    513   numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);       
    514 
    515   hTimeDiff_on_timebin->Fit("exp","RQ0");
    516   deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
    517 
    518   cout.precision(5);
    519   cout << "2-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" <<  f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl;
    520   cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/TMath::Abs(f1->GetParameter(1)) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1)))  << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
    521 
    522   cout.precision(10);
    523   cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
    524   cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << " widthTimeBinOn " <<  widthTimeBinOn[numberTimeBins-1] << endl;
     338  numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90);
     339  numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90);         
     340
     341//   hTimeDiff_on_timebin->Fit("exp","RQ0");
     342//   deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries());
     343
     344//   cout.precision(5);
     345//   cout << "2-> Exp(" << f1->GetParameter(0) << " +
     346//   cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/TMath::Abs(f1->GetParameter(1)) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1)))  << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl;
     347  deadFractionOn[numberTimeBins-1] = 1.;
    525348
    526349  alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
    … …  
    529352  timediffHistoOn.AddLast(hTimeDiff_on_timebin);
    530353
    531   //  cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
    532 
    533   loop_on.PostProcess();
    534  
    535   tlist_on.PrintStatistics();
    536  
    537   for (UInt_t bin=0; bin<numberTimeBins; bin++)
    538     cout << bin << " timeOn " << timeOn[bin] << " min-max zenithOn " << zenithMinimumOn[bin] << "-" <<  zenithMaximumOn[bin] << " min-max cos(zenithOn) " << TMath::Cos(zenithMinimumOn[bin]*kConvDegToRad) << "-" <<  TMath::Cos(zenithMaximumOn[bin]*kConvDegToRad) << " numberOnEvents " << numberOnEvents[bin] << " numberOnEventsAfterCleaning " << numberOnEventsAfterCleaning[bin] << " deadFractionOn " << deadFractionOn[bin] << endl;
    539 
    540 
    541      
    542 
    543354  //-----------------------OFF------------------------
    544355
    … …  
    552363  TArrayD numberBkgEventsToNormOff(numberTimeBins);
    553364  TArrayD onOffNormFactor(numberTimeBins);
     365  TArrayD onOffNormFactorWithEvents(numberTimeBins);
     366  TArrayD onOffNormFactorWithLiveTime(numberTimeBins);
    554367
    555368  TArrayD numberExcessEvents(numberTimeBins);
    … …  
    590403  MRawEvtHeader evtheader_off;
    591404  MTime              time_off;
     405  MLiveTime      livetime_off;
    592406  MHillas          hillas_off;
    593407  MHillasSrc    hillassrc_off;
    … …  
    599413  plist_off.AddToList(&evtheader_off);
    600414  plist_off.AddToList(&time_off);
     415  plist_off.AddToList(&livetime_off);
    601416  plist_off.AddToList(&hillas_off);
    602417  plist_off.AddToList(&hillassrc_off);
    … …  
    627442  MHillasSrcCalc csrc_off;
    628443 
    629   //tasklist
     444  MLiveTimeCalc livetimecalc_off;
     445
     446 //tasklist
    630447  tlist_off.AddToList(&read_off);
     448  tlist_off.AddToList(&livetimecalc_off);
    631449  tlist_off.AddToList(&srcplace_timebin); // This is just to run the preprocess correctely
    632450  tlist_off.AddToList(&csrc_off); // This is just to run the preprocess correctely
    … …  
    714532            }
    715533        }
    716      
    717       //       if (zenith >  zenithMaximumOn[numberTimeBins-1])
    718       //        {
    719       //          cout << "Exit loop because we arrive to zenith = " << zenith << endl;
    720       //          break;
    721       //        }
    722534    }
    723535 
    … …  
    734546  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    735547    {
    736      
     548      cout.precision(5);
     549      cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl;
     550    }
     551  livetime_off.Print("last");
     552  livetime_off.Print("all");
     553  cout << "livetime_off.GetLiveTime() " << livetime_off.GetLiveTime() << endl;
     554
     555  for (UInt_t bin=0; bin<numberTimeBins; bin++)
     556    {
     557      timeOff[bin] = livetime_off.GetLiveTime();
     558
     559      //
    737560      meanTriggerRateOn[bin] = numberOnEventsAfterCleaning[bin]/timeOn[bin];
    738561      errorMeanTriggerRateOn[bin] = TMath::Sqrt(numberOnEventsAfterCleaning[bin])/timeOn[bin];
    … …  
    745568      widthTimeBinOnInSec[bin] = widthTimeBinOn[bin]*kDay;
    746569
    747       numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
    748       numberBkgEventsToNormOff[bin] =  (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
     570      numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90);
     571      numberBkgEventsToNormOff[bin] =  (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90);
    749572      if (numberOffEvents[bin] != 0 && numberBkgEventsToNormOff[bin] != 0)
    750573        {
    751           onOffNormFactor[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin];
     574          onOffNormFactorWithEvents[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin];
     575          onOffNormFactorWithLiveTime[bin] = timeOn[bin]/timeOff[bin];
     576          onOffNormFactor[bin] = onOffNormFactorWithEvents[bin];
    752577          numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin];
    753578          errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]);
    … …  
    756581          numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.);
    757582          errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin];
     583
    758584        }
    759585    }
    … …  
    762588    {
    763589      cout.precision(5);
    764       cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl;
     590      cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOnInSec[bin] << "-" << widthTimeBinOnInSec[bin] << endl;
    765591      cout << " numberOnEvents\t " << numberOnEvents[bin] << endl;
    766592      cout << " numberOffEvents\t " << numberOffEvents[bin] << endl;
    767593      cout << " numberBkgEventsToNormOn\t " << numberBkgEventsToNormOn[bin] << endl;
    768594      cout << " numberBkgEventsToNormOff\t " << numberBkgEventsToNormOff[bin] << endl;
    769       cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl;
     595      cout << " onOffNormFactorWithEvents\t " << onOffNormFactorWithEvents[bin] << endl;
     596      cout << " onOffNormFactorWithLiveTime\t " << onOffNormFactorWithLiveTime[bin] << endl;
    770597      cout << " numberExcessEvents\t " <<  numberExcessEvents[bin] <<  " +- " << errorNumberExcessEvents[bin] << endl;
    771598      cout << " deadFraction\t" << deadFractionOn[bin] << endl;
    … …  
    796623  TGraphErrors* cosmicrategraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),meanTriggerRateOn.GetArray(),widthTimeBinOnInSec.GetArray(),errorMeanTriggerRateOn.GetArray());
    797624  cosmicrategraph->SetTitle("Cosmic Rate");
    798   cosmicrategraph->SetMinimum(0.);
    799625  cosmicrategraph->SetMarkerStyle(21);
    800626  cosmicrategraph->SetMarkerSize(0.03);
    … …  
    810636  Float_t maxAlphaHistoHeight = 0;
    811637 
    812   for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
     638  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    813639    {
    814640      for (UInt_t i=1; i<=nbins_abs; i++)
    … …  
    816642          maxAlphaHistoHeight = ((TH1F*)alphaHistoOn[bin])->GetBinContent(i);
    817643    }     
    818  
     644 
     645  cout << "maxAlphaHistoHeight " << maxAlphaHistoHeight << endl;
    819646
    820647  for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
    … …  
    897724  // ############################################################################
    898725
    899   Double_t norm_on_abs  = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
    900   Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
    901   Double_t norm_off_abs  = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
    902   Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
     726  Double_t norm_on_abs  = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90);
     727  Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90);
     728  Double_t norm_off_abs  = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90);
     729  Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90);
    903730  Double_t norm = norm_on_abs/norm_off_abs;
    904731
    … …  
    12961123
    12971124
    1298 
     1125Double_t ChiSquareNDof(TH1D *h1, TH1D *h2)
     1126{
     1127    Double_t chiq = 0.;
     1128    Double_t chi;
     1129    Double_t error;
     1130    Int_t nbinsnozero = 0;
     1131
     1132    Int_t nbins = h1->GetNbinsX();
     1133    if (nbins != h2->GetNbinsX() || nbins == 0)
     1134        return -1;
     1135
     1136    for (UInt_t bin=1; bin<=nbins; bin++)
     1137    {
     1138        error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) +
     1139                           h2->GetBinError(bin)*h2->GetBinError(bin));
     1140        if (error != 0)
     1141        {
     1142            chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error;
     1143            chiq += chi*chi;
     1144            nbinsnozero++;
     1145        }
     1146    }
     1147
     1148    return (nbinsnozero>0?chiq/nbinsnozero:0);
     1149}
     1150
     1151Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[])
     1152{
     1153
     1154  Int_t result = -1;
     1155
     1156  Int_t lowerbin = 0;
     1157  Int_t upperbin = numberSizeBins;
     1158  Int_t bin;
     1159
     1160  Int_t count = 0;
     1161
     1162  if (size >= sizeBins[0])
     1163    {
     1164      while (upperbin - lowerbin > 1 && count++<=numberSizeBins)
     1165        {
     1166          bin = (upperbin+lowerbin)/2;
     1167          if (size >= sizeBins[bin])
     1168            lowerbin = bin;
     1169          else
     1170            upperbin = bin;
     1171        }
     1172      result = count<=numberSizeBins?lowerbin:-1;
     1173    }
     1174
     1175  return result;
     1176
     1177}
     1178

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

    r4310 r4429  
    7171
    7272  MJPedestal pedloop0;
     73  pedloop0.SetOutputPath("./");
    7374  pedloop0.SetInput(&pedcaliter);
    7475  pedloop0.SetExtractor(&extractor);
    … …  
    8687  MJCalibration calloop;
    8788
     89  calloop.SetOutputPath("./");
    8890  calloop.SetExtractor(&extractor);
    8991  calloop.SetInput(&caliter);
    … …  
    101103  // First Compute the pedestals
    102104  MJPedestal pedloop;
     105  pedloop.SetOutputPath("./");
    103106  pedloop.SetInput(&pediter);
    104107
    … …  
    225228  read4.AddTree("Drive");
    226229  static_cast<MRead&>(read4).AddFiles(datiter);
     230  read4.AddToBranchList("MReportCurrents.*");
    227231
    228232  // Task that use Slow control information
    … …  
    252256  MSrcPosFromStars srcposfromstar;
    253257   
    254   // set bad pixels
    255   MBlindPixelCalc   blind;
    256   MBlindPixelCalc   blind2;
    257   const Short_t x[16] = {330,395,329,396,389,
    258                          323,388,322,384,385,
    259                          386,387,321,320,319,
    260                          394};
    261   const TArrayS bp(16,(Short_t*)x);
    262   blind.SetPixelIndices(bp);
    263   blind2.SetPixelIndices(bp);
     258  MBadPixelsTreat   interpolatebadpixels;
     259  interpolatebadpixels.SetUseInterpolation();
    264260 
    265261  MImgCleanStd      clean(lcore,ltail);
    … …  
    315311  tlist4.AddToList(&extractor,"Events");
    316312  tlist4.AddToList(&photcalc,"Events");
    317   //tlist4.AddToList(&blind,"Events");
     313  if(calflag==11)
     314    tlist4.AddToList(&interpolatebadpixels);
    318315  tlist4.AddToList(&clean,"Events");
    319316
    … …  
    332329  //tlist4.AddToList(&blind2,"Events");
    333330  tlist4.AddToList(&hcalc,"Events");
    334   //  tlist4.AddToList(&srcposcalc,"Events");
    335331  tlist4.AddToList(&csrc1,"Events");
    336332  tlist4.AddToList(&write,"Events");
    … …  
    350346
    351347  tlist4.PrintStatistics();   
     348
    352349}
    353350//-------------------------------------------------------------------------------
    … …  
    476473    }
    477474
    478    pedcaliter.Reset();
    479    pediter.Reset();
    480    caliter.Reset();
    481    datiter.Reset();
     475  pedcaliter.Reset();
     476  pediter.Reset();
     477  caliter.Reset();
     478  datiter.Reset();
    482479  TString pfile;
    483480

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

    r4259 r4429  
    66                     Int_t iniRun=0,
    77                     Int_t finRun=99999,
    8                      TString filename="/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/AllCrabNebula_3015_kDefault_MispE.root",
     8                     TString filename="/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/AllCrabNebula_3015_kDefault_MispE.root"
    99)
    1010{   
Note: See TracChangeset for help on using the changeset viewer.