Ignore:
Timestamp:
07/20/04 11:57:23 (20 years ago)
Author:
jlopez
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mtemp/mifae/macros
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

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

    r4397 r4408  
    6767  // Configuration
    6868  const Bool_t debug = kFALSE;
    69   const Double_t timebin = 480.; //[sec]
    70   TString psname = "./20040420_Mrk421.480s.ps";
     69  const Double_t timebin = 600.; //[sec]
     70  TString psname = "./20040420_Mrk421.600s.ps";
    7171
    7272  //Constanst
    7373  const Double_t kConvDegToRad = TMath::Pi()/180.;
    74   const Double_t kSec = 1e6;  //[sec/microsec]
     74  const Double_t kSec = 1e3;  //[sec/milisec]
    7575  const Double_t kDay = 24.*60.*60.;  //[Day/sec]
    7676
    7777  UInt_t numberTimeBins = 1;
    7878  TArrayI numberOnEvents(1);
     79  TArrayD numberOnEventsAfterCleaning(1);
    7980  TArrayD numberBkgEventsToNormOn(1);
    8081  TArrayD timeOn(1);
    81 
     82 
    8283  TArrayD meanTimeBinOn(1);
    8384  TArrayD widthTimeBinOn(1);
    … …  
    8586  TArrayD zenithMinimumOn(1);
    8687  TArrayD zenithMaximumOn(1);
     88
     89  TArrayD deadFractionOn(1);
    8790 
    8891  TObjArray coszenithHistoOn;
    8992  TObjArray alphaHistoOn;
    9093  TObjArray srcposHistoOn;
    91  
    92 
     94  TObjArray timediffHistoOn;
    9395 
    9496  const Int_t numberSizeBins = 4;
    … …  
    101103  Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 };
    102104  Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
    103   Double_t distmin[numberSizeBins]   = { 0.30, 0.30, 0.30, 0.30,};
    104   Double_t distmax[numberSizeBins]   = { 1.20, 1.20, 1.20, 1.20,};
     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 };
    105107
    106108  //alpha plot integration for gammas
    … …  
    170172  TH1F *hAlpha_on = new TH1F ("hAlphaOn","",nbins,minalpha,maxalpha);
    171173
    172   Int_t nbins_srcpos = 400;
    173   Double_t minsrcpos = -600.;
    174   Double_t maxsrcpos =  600.;  //[mm]  //!!! position precition 3mm ~ 0.01 deg
     174  Int_t nbins_srcpos = 200;
     175  Double_t minsrcpos = -200.;
     176  Double_t maxsrcpos =  200.;  //[mm]  //!!! position precition 3mm ~ 0.01 deg
    175177  TH2F *hSrcPos_on = new TH2F ("hSrcPosOn","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
    176178
    … …  
    224226  Bool_t flag = kFALSE;
    225227
     228  numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
    226229  zenithMinimumOn[numberTimeBins-1] = 100.;
    227230  zenithMaximumOn[numberTimeBins-1] = 0.;
    … …  
    241244  TString coszenithTitle =  Form("%s%02i","hCosZenithOn",numberTimeBins-1);
    242245  TH1F *hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
     246
     247  Int_t nbins_timediff = 2000;
     248  Double_t mintimediff = 0.; 
     249  Double_t maxtimediff = 40.; 
     250  TString timediffTitle =  Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
     251  TH1F *hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
     252  TF1 *f1 = new TF1("exp","expo",mintimediff,maxtimediff);           
    243253
    244254  while(tlist_on.Process())
    … …  
    248258      UInt_t   evt = evtheader_on.GetDAQEvtNumber();
    249259      Double_t mjd = time_on.GetMjd();
    250 
     260     
     261      numberOnEventsAfterCleaning[numberTimeBins-1]++;
     262                                 
    251263      if (mjd == 0)
    252         {
    253          
     264      {
     265       
    254266          if (debug)
    255267            {
    … …  
    362374              mjdFirstEvent = mjd;
    363375            }
     376          else if (mjdFirstEventofBin != 0)
     377            hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec);
    364378         
    365379          if (mjdFirstEventofBin == 0)
    … …  
    368382              mjdFirstEventofBin = mjd;
    369383            }
    370          
     384
    371385          evtLastEvent = evt;
    372386          runLastEvent = run;
    373387          mjdLastEvent = mjd;
    374388          mtimeLastEvent = time_on;
    375 
    376      
    377389     
    378390          Double_t size = hillas_on.GetSize();
    379391          Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg();
    380392          Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg();
    381           Double_t dist = hillassrc_on.GetDist()*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();
    382397          Double_t alpha = hillassrc_on.GetAlpha();
    383398          Double_t srcposx = srcpos_on.GetX();
    … …  
    394409                  if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
    395410                    {
    396                       if (dist > distmin[sizebin] && dist < distmax[sizebin])
     411                      if (dist > distmin[sizebin] && centerdist < distmax[sizebin])
    397412                        {     
    398413
    … …  
    439454              mjdFirstEvent = mjd;
    440455              mjdFirstEventofBin = mjd;
    441              
     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
    442462              alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
    443463              srcposHistoOn.AddLast(hSrcPos_on_timebin);
    444464              coszenithHistoOn.AddLast(hCosZenith_on_timebin);
     465              timediffHistoOn.AddLast(hTimeDiff_on_timebin);
    445466             
    446467//            cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
    … …  
    456477              zenithMinimumOn.Set(numberTimeBins);
    457478              zenithMaximumOn.Set(numberTimeBins);
    458              
     479              deadFractionOn.Set(numberTimeBins);
     480              numberOnEventsAfterCleaning.Set(numberTimeBins);
     481
    459482              timeOn[numberTimeBins-1] = 0.;
    460483              zenithMinimumOn[numberTimeBins-1] = 100.;
    461484              zenithMaximumOn[numberTimeBins-1] = 0.;
     485              numberOnEventsAfterCleaning[numberTimeBins-1] = 0;
    462486             
    463487              alphaTitle =  Form("%s%02i","hAlphaOn",numberTimeBins-1);         
    … …  
    470494              hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);
    471495
     496              timediffTitle =  Form("%s%02i","hTimeDiffOn",numberTimeBins-1);
     497              hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff);
     498
    472499              flag = kTRUE;
    473500          }
    474         }       
     501       
     502        }
     503       
    475504    }
    476505
    … …  
    484513  numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);       
    485514
     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
    486522  cout.precision(10);
    487523  cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
    … …  
    491527  srcposHistoOn.AddLast(hSrcPos_on_timebin);
    492528  coszenithHistoOn.AddLast(hCosZenith_on_timebin);
     529  timediffHistoOn.AddLast(hTimeDiff_on_timebin);
    493530
    494531  //  cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;
    … …  
    499536 
    500537  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    501     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] << endl;
    502 
     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     
    503542
    504543  //-----------------------OFF------------------------
    … …  
    616655      Double_t width = hillas_off.GetWidth()*geomcam.GetConvMm2Deg();
    617656      Double_t length = hillas_off.GetLength()*geomcam.GetConvMm2Deg();
     657      Double_t meanx = hillas_off.GetMeanX()*geomcam.GetConvMm2Deg();
     658      Double_t meany = hillas_off.GetMeanY()*geomcam.GetConvMm2Deg();
     659      Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany);
    618660     
    619661      Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
    … …  
    625667              if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
    626668                {
    627                   //General histos
    628 
    629                   if (dist > distmin[sizebin] && dist < distmax[sizebin])
     669                  if (centerdist < distmax[sizebin])
    630670                    {
    631                       hSize_off->Fill(log10(size));
    632                       hWidth_off->Fill(width);
    633                       hLength_off->Fill(length);
     671                      //General histos
    634672
    635673                      srcplace_timebin.SetPositionHisto(hSrcPos_on);
    636                       srcplace_timebin.CallProcess();
    637                       csrc_off.CallProcess();
    638                  
    639                       Double_t dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
    640                       Double_t alpha = hillassrc_off.GetAlpha();
    641                       Double_t srcposx = srcpos_off.GetX();
    642                       Double_t srcposy = srcpos_off.GetY();
    643 
    644                       hDist_off->Fill(dist);
    645                       hAlpha_off_abs->Fill(TMath::Abs(alpha));
    646                       hAlpha_off->Fill(alpha);
    647                       hSrcPos_off->Fill(srcposx,srcposy);
    648                     }
    649 
    650                   // Time bin histos
    651                   for (UInt_t bin=0; bin<numberTimeBins; bin++)
    652                     {
    653                       srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
    654674                      srcplace_timebin.CallProcess();
    655675                      csrc_off.CallProcess();
    … …  
    660680                      Double_t srcposy = srcpos_off.GetY();
    661681                     
    662                       if (dist > distmin[sizebin] && dist < distmax[sizebin])
     682                      if (dist > distmin[sizebin] )
    663683                        {
    664                           ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
    665                           ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
     684                          hSize_off->Fill(log10(size));
     685                          hWidth_off->Fill(width);
     686                          hLength_off->Fill(length);
     687                         
     688                          hDist_off->Fill(dist);
     689                          hAlpha_off_abs->Fill(TMath::Abs(alpha));
     690                          hAlpha_off->Fill(alpha);
     691                          hSrcPos_off->Fill(srcposx,srcposy);
     692                        }
     693                     
     694                      // Time bin histos
     695                      for (UInt_t bin=0; bin<numberTimeBins; bin++)
     696                        {
     697                          srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
     698                          srcplace_timebin.CallProcess();
     699                          csrc_off.CallProcess();
     700                         
     701                          dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
     702                          alpha = hillassrc_off.GetAlpha();
     703                          srcposx = srcpos_off.GetX();
     704                          srcposy = srcpos_off.GetY();
     705                         
     706                          if (dist > distmin[sizebin])
     707                            {
     708                              ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
     709                              ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
     710                            }
    666711                        }
    667712                    }
    … …  
    684729  TArrayD widthTimeBinOnInSec(numberTimeBins);
    685730 
     731  TArrayD meanTriggerRateOn(numberTimeBins);
     732  TArrayD errorMeanTriggerRateOn(numberTimeBins);
     733
    686734  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    687735    {
    688736     
     737      meanTriggerRateOn[bin] = numberOnEventsAfterCleaning[bin]/timeOn[bin];
     738      errorMeanTriggerRateOn[bin] = TMath::Sqrt(numberOnEventsAfterCleaning[bin])/timeOn[bin];
     739
    689740      meanTimeBinOnInSec[bin]  = (meanTimeBinOn[bin]-(Int_t)meanTimeBinOn[bin]);
    690741      if (meanTimeBinOnInSec[bin] > 0.5)
    … …  
    701752          numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin];
    702753          errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]);
    703           numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin];
     754          numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.);
    704755          errorNumberExcessEventsPerSec[bin] = errorNumberExcessEvents[bin]/timeOn[bin];
    705           numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin];
     756          numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.);
    706757          errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin];
    707758        }
    … …  
    718769      cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl;
    719770      cout << " numberExcessEvents\t " <<  numberExcessEvents[bin] <<  " +- " << errorNumberExcessEvents[bin] << endl;
     771      cout << " deadFraction\t" << deadFractionOn[bin] << endl;
    720772      cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl;
     773      cout << " Trigger Rate\t " << meanTriggerRateOn[bin] << " +- " << errorMeanTriggerRateOn[bin] << endl;
    721774    }
    722775
    … …  
    739792  c0->Print(openpsname);
    740793     
     794  TCanvas *c00 = new TCanvas;
     795  c00->cd(1);
     796  TGraphErrors* cosmicrategraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),meanTriggerRateOn.GetArray(),widthTimeBinOnInSec.GetArray(),errorMeanTriggerRateOn.GetArray());
     797  cosmicrategraph->SetTitle("Cosmic Rate");
     798  cosmicrategraph->SetMinimum(0.);
     799  cosmicrategraph->SetMarkerStyle(21);
     800  cosmicrategraph->SetMarkerSize(0.03);
     801  cosmicrategraph->Draw("AP");
     802  cosmicrategraph->GetYaxis()->SetTitle("[Hz]");
     803  cosmicrategraph->GetXaxis()->SetTitle("UTC Time");
     804  cosmicrategraph->GetXaxis()->SetTimeDisplay(1);
     805  c00->Print(psname);
     806
    741807  TCanvas** c = new TCanvas*[numberTimeBins];
    742808 
     809  //Compute the maximum of all hAlphaOn histograms
     810  Float_t maxAlphaHistoHeight = 0;
     811 
     812  for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
     813    {
     814      for (UInt_t i=1; i<=nbins_abs; i++)
     815        if (((TH1F*)alphaHistoOn[bin])->GetBinContent(i) > maxAlphaHistoHeight)
     816          maxAlphaHistoHeight = ((TH1F*)alphaHistoOn[bin])->GetBinContent(i);
     817    }     
     818 
     819
    743820  for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
    744821    {
    … …  
    758835      ((TH1F*)alphaHistoOff[bin])->SetLineColor(46);
    759836      ((TH1F*)alphaHistoOff[bin])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
     837      ((TH1F*)alphaHistoOn[bin])->SetMaximum(maxAlphaHistoHeight*1.2);
    760838      ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.);
    761839      alphaTitle =  Form("%s%02i","hAlphaOnOff",bin);
    … …  
    786864  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46);
    787865  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
     866  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMaximum(maxAlphaHistoHeight*1.2);
    788867  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.);
    789868  alphaTitle =  Form("%s%02i","hAlphaOnOff",numberTimeBins-1);
    … …  
    798877  c[numberTimeBins-1]->Print(psname);
    799878
    800 //   TFile input("./prueba.root", "RECREATE");
    801 
    802 //    for (UInt_t bin=0; bin<numberTimeBins; bin++)
    803 //      {
    804 //        ((TH1F*)alphaHistoOn[bin])->Write();
    805 //        ((TH1F*)alphaHistoOff[bin])->Write();
    806 //        ((TH2F*)srcposHistoOn[bin])->Write();
    807 //        ((TH2F*)srcposHistoOff[bin])->Write();
    808 //        ((TH1F*)coszenithHistoOn[bin])->Write();
    809 //      }
    810 
    811 //    input.Close();
     879  //  TString rootname = psname.ReplaceAll(".ps",".root");
     880  TString rootname = "./prueba.root";
     881  TFile input(rootname, "RECREATE");
     882
     883   for (UInt_t bin=0; bin<numberTimeBins; bin++)
     884     {
     885       ((TH1F*)alphaHistoOn[bin])->Write();
     886       ((TH2F*)srcposHistoOn[bin])->Write();
     887       ((TH1F*)coszenithHistoOn[bin])->Write();
     888       ((TH1F*)timediffHistoOn[bin])->Write();
     889       ((TH1F*)alphaHistoOff[bin])->Write();
     890       ((TH2F*)srcposHistoOff[bin])->Write();
     891    }
     892 
     893  input.Close();
    812894
    813895  // ############################################################################
Note: See TracChangeset for help on using the changeset viewer.