Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc (revision 4407) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc (revision 4408) @@ -148,8 +148,12 @@ fLiveTime->SetRealTime(mean,width); - *fLog << inf << GetName() << ": New time bin" << endl; - Print("all"); - *fLog << GetName() << ": First event time " << setprecision(10) << fFirstEventMjd << setprecision(5) << endl; - fLiveTime->Print("last"); + + if (Debug) + { + *fLog << inf << GetName() << ": New time bin" << endl; + Print("all"); + *fLog << GetName() << ": First event time " << setprecision(10) << fFirstEventMjd << setprecision(5) << endl; + fLiveTime->Print("last"); + } fLiveTime->AddBin(); Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C (revision 4407) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C (revision 4408) @@ -67,17 +67,18 @@ // Configuration const Bool_t debug = kFALSE; - const Double_t timebin = 480.; //[sec] - TString psname = "./20040420_Mrk421.480s.ps"; + const Double_t timebin = 600.; //[sec] + TString psname = "./20040420_Mrk421.600s.ps"; //Constanst const Double_t kConvDegToRad = TMath::Pi()/180.; - const Double_t kSec = 1e6; //[sec/microsec] + const Double_t kSec = 1e3; //[sec/milisec] const Double_t kDay = 24.*60.*60.; //[Day/sec] UInt_t numberTimeBins = 1; TArrayI numberOnEvents(1); + TArrayD numberOnEventsAfterCleaning(1); TArrayD numberBkgEventsToNormOn(1); TArrayD timeOn(1); - + TArrayD meanTimeBinOn(1); TArrayD widthTimeBinOn(1); @@ -85,10 +86,11 @@ TArrayD zenithMinimumOn(1); TArrayD zenithMaximumOn(1); + + TArrayD deadFractionOn(1); TObjArray coszenithHistoOn; TObjArray alphaHistoOn; TObjArray srcposHistoOn; - - + TObjArray timediffHistoOn; const Int_t numberSizeBins = 4; @@ -101,6 +103,6 @@ Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 }; Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 }; - Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30,}; - Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20,}; + Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30 }; + Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20 }; //alpha plot integration for gammas @@ -170,7 +172,7 @@ TH1F *hAlpha_on = new TH1F ("hAlphaOn","",nbins,minalpha,maxalpha); - Int_t nbins_srcpos = 400; - Double_t minsrcpos = -600.; - Double_t maxsrcpos = 600.; //[mm] //!!! position precition 3mm ~ 0.01 deg + Int_t nbins_srcpos = 200; + Double_t minsrcpos = -200.; + Double_t maxsrcpos = 200.; //[mm] //!!! position precition 3mm ~ 0.01 deg TH2F *hSrcPos_on = new TH2F ("hSrcPosOn","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos); @@ -224,4 +226,5 @@ Bool_t flag = kFALSE; + numberOnEventsAfterCleaning[numberTimeBins-1] = 0; zenithMinimumOn[numberTimeBins-1] = 100.; zenithMaximumOn[numberTimeBins-1] = 0.; @@ -241,4 +244,11 @@ TString coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1); TH1F *hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith); + + Int_t nbins_timediff = 2000; + Double_t mintimediff = 0.; + Double_t maxtimediff = 40.; + TString timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1); + TH1F *hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff); + TF1 *f1 = new TF1("exp","expo",mintimediff,maxtimediff); while(tlist_on.Process()) @@ -248,8 +258,10 @@ UInt_t evt = evtheader_on.GetDAQEvtNumber(); Double_t mjd = time_on.GetMjd(); - + + numberOnEventsAfterCleaning[numberTimeBins-1]++; + if (mjd == 0) - { - + { + if (debug) { @@ -362,4 +374,6 @@ mjdFirstEvent = mjd; } + else if (mjdFirstEventofBin != 0) + hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec); if (mjdFirstEventofBin == 0) @@ -368,16 +382,17 @@ mjdFirstEventofBin = mjd; } - + evtLastEvent = evt; runLastEvent = run; mjdLastEvent = mjd; mtimeLastEvent = time_on; - - Double_t size = hillas_on.GetSize(); Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg(); Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg(); - Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg(); + Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg(); + Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg(); + Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany); + Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg(); Double_t alpha = hillassrc_on.GetAlpha(); Double_t srcposx = srcpos_on.GetX(); @@ -394,5 +409,5 @@ if (length > lengthmin[sizebin] && length < lengthmax[sizebin]) { - if (dist > distmin[sizebin] && dist < distmax[sizebin]) + if (dist > distmin[sizebin] && centerdist < distmax[sizebin]) { @@ -439,8 +454,14 @@ mjdFirstEvent = mjd; mjdFirstEventofBin = mjd; - + + hTimeDiff_on_timebin->Fit("exp","RQ0"); + deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); + cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl; + 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; + alphaHistoOn.AddLast(hAlpha_on_abs_timebin); srcposHistoOn.AddLast(hSrcPos_on_timebin); coszenithHistoOn.AddLast(hCosZenith_on_timebin); + timediffHistoOn.AddLast(hTimeDiff_on_timebin); // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl; @@ -456,8 +477,11 @@ zenithMinimumOn.Set(numberTimeBins); zenithMaximumOn.Set(numberTimeBins); - + deadFractionOn.Set(numberTimeBins); + numberOnEventsAfterCleaning.Set(numberTimeBins); + timeOn[numberTimeBins-1] = 0.; zenithMinimumOn[numberTimeBins-1] = 100.; zenithMaximumOn[numberTimeBins-1] = 0.; + numberOnEventsAfterCleaning[numberTimeBins-1] = 0; alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1); @@ -470,7 +494,12 @@ hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith); + timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1); + hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff); + flag = kTRUE; } - } + + } + } @@ -484,4 +513,11 @@ numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1); + hTimeDiff_on_timebin->Fit("exp","RQ0"); + deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); + + cout.precision(5); + cout << "2-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl; + 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; + cout.precision(10); cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl; @@ -491,4 +527,5 @@ srcposHistoOn.AddLast(hSrcPos_on_timebin); coszenithHistoOn.AddLast(hCosZenith_on_timebin); + timediffHistoOn.AddLast(hTimeDiff_on_timebin); // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl; @@ -499,6 +536,8 @@ for (UInt_t bin=0; bin lengthmin[sizebin] && length < lengthmax[sizebin]) { - //General histos - - if (dist > distmin[sizebin] && dist < distmax[sizebin]) + if (centerdist < distmax[sizebin]) { - hSize_off->Fill(log10(size)); - hWidth_off->Fill(width); - hLength_off->Fill(length); + //General histos srcplace_timebin.SetPositionHisto(hSrcPos_on); - srcplace_timebin.CallProcess(); - csrc_off.CallProcess(); - - Double_t dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg(); - Double_t alpha = hillassrc_off.GetAlpha(); - Double_t srcposx = srcpos_off.GetX(); - Double_t srcposy = srcpos_off.GetY(); - - hDist_off->Fill(dist); - hAlpha_off_abs->Fill(TMath::Abs(alpha)); - hAlpha_off->Fill(alpha); - hSrcPos_off->Fill(srcposx,srcposy); - } - - // Time bin histos - for (UInt_t bin=0; bin distmin[sizebin] && dist < distmax[sizebin]) + if (dist > distmin[sizebin] ) { - ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha)); - ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy); + hSize_off->Fill(log10(size)); + hWidth_off->Fill(width); + hLength_off->Fill(length); + + hDist_off->Fill(dist); + hAlpha_off_abs->Fill(TMath::Abs(alpha)); + hAlpha_off->Fill(alpha); + hSrcPos_off->Fill(srcposx,srcposy); + } + + // Time bin histos + for (UInt_t bin=0; bin distmin[sizebin]) + { + ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha)); + ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy); + } } } @@ -684,7 +729,13 @@ TArrayD widthTimeBinOnInSec(numberTimeBins); + TArrayD meanTriggerRateOn(numberTimeBins); + TArrayD errorMeanTriggerRateOn(numberTimeBins); + for (UInt_t bin=0; bin 0.5) @@ -701,7 +752,7 @@ numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin]; errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]); - numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin]; + numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.); errorNumberExcessEventsPerSec[bin] = errorNumberExcessEvents[bin]/timeOn[bin]; - numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]; + numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.); errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin]; } @@ -718,5 +769,7 @@ cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl; cout << " numberExcessEvents\t " << numberExcessEvents[bin] << " +- " << errorNumberExcessEvents[bin] << endl; + cout << " deadFraction\t" << deadFractionOn[bin] << endl; cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl; + cout << " Trigger Rate\t " << meanTriggerRateOn[bin] << " +- " << errorMeanTriggerRateOn[bin] << endl; } @@ -739,6 +792,30 @@ c0->Print(openpsname); + TCanvas *c00 = new TCanvas; + c00->cd(1); + TGraphErrors* cosmicrategraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),meanTriggerRateOn.GetArray(),widthTimeBinOnInSec.GetArray(),errorMeanTriggerRateOn.GetArray()); + cosmicrategraph->SetTitle("Cosmic Rate"); + cosmicrategraph->SetMinimum(0.); + cosmicrategraph->SetMarkerStyle(21); + cosmicrategraph->SetMarkerSize(0.03); + cosmicrategraph->Draw("AP"); + cosmicrategraph->GetYaxis()->SetTitle("[Hz]"); + cosmicrategraph->GetXaxis()->SetTitle("UTC Time"); + cosmicrategraph->GetXaxis()->SetTimeDisplay(1); + c00->Print(psname); + TCanvas** c = new TCanvas*[numberTimeBins]; + //Compute the maximum of all hAlphaOn histograms + Float_t maxAlphaHistoHeight = 0; + + for (UInt_t bin=0; binGetBinContent(i) > maxAlphaHistoHeight) + maxAlphaHistoHeight = ((TH1F*)alphaHistoOn[bin])->GetBinContent(i); + } + + for (UInt_t bin=0; binSetLineColor(46); ((TH1F*)alphaHistoOff[bin])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + ((TH1F*)alphaHistoOn[bin])->SetMaximum(maxAlphaHistoHeight*1.2); ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.); alphaTitle = Form("%s%02i","hAlphaOnOff",bin); @@ -786,4 +864,5 @@ ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46); ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMaximum(maxAlphaHistoHeight*1.2); ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.); alphaTitle = Form("%s%02i","hAlphaOnOff",numberTimeBins-1); @@ -798,16 +877,19 @@ c[numberTimeBins-1]->Print(psname); -// TFile input("./prueba.root", "RECREATE"); - -// for (UInt_t bin=0; binWrite(); -// ((TH1F*)alphaHistoOff[bin])->Write(); -// ((TH2F*)srcposHistoOn[bin])->Write(); -// ((TH2F*)srcposHistoOff[bin])->Write(); -// ((TH1F*)coszenithHistoOn[bin])->Write(); -// } - -// input.Close(); + // TString rootname = psname.ReplaceAll(".ps",".root"); + TString rootname = "./prueba.root"; + TFile input(rootname, "RECREATE"); + + for (UInt_t bin=0; binWrite(); + ((TH2F*)srcposHistoOn[bin])->Write(); + ((TH1F*)coszenithHistoOn[bin])->Write(); + ((TH1F*)timediffHistoOn[bin])->Write(); + ((TH1F*)alphaHistoOff[bin])->Write(); + ((TH2F*)srcposHistoOff[bin])->Write(); + } + + input.Close(); // ############################################################################ Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/runlivetime.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/runlivetime.C (revision 4408) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/runlivetime.C (revision 4408) @@ -0,0 +1,96 @@ +void runlivetime(TString filename) +{ + + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + + // + // Make a loop only for the ON data: + // + + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + //Containers + + MLiveTime livetime; + + plist.AddToList(&livetime); + + // + //tasks + // + + MReadTree read("Parameters"); + read.DisableAutoScheme(); + + TString tmpfile = "./runlivetime.ls.tmp"; + TString cmd = "ls " + filename + " > " + tmpfile; + gSystem->Exec(cmd); + TString tmpline; + ifstream in(tmpfile); + while(1) + { + in >> tmpline; + if(in.eof()) + break; + read.AddFile(tmpline); + } + TString cmd = "rm " + tmpfile; + gSystem->Exec(cmd); + + Double_t timebin = 500.; //[sec] + MLiveTimeCalc livetimecalc; + livetimecalc.SetRealTimeBinSize(timebin); + + tlist.AddToList(&read); + tlist.AddToList(&livetimecalc); + + // + // Create and setup the eventloop + // + MEvtLoop loop; + loop.SetParList(&plist); + + + if (!loop.Eventloop()) + return; + +// if (!loop.PreProcess()) +// return; + +// while(loop.Process()) +// {} + +// loop.PostProcess(); + + tlist.PrintStatistics(); + + UInt_t numbertimebins = livetime.GetNumberTimeBins(); + TArrayD y(numbertimebins); + TArrayD erry(numbertimebins); + + for (UInt_t bin=0; bincd(1); + TGraphErrors *graph = new TGraphErrors(numbertimebins,livetime.GetMeanRealTimeArray(),y.GetArray(),livetime.GetWidthRealTimeArray(),erry.GetArray()); + graph->SetTitle(""); + graph->SetMinimum(0.); + graph->SetMarkerStyle(21); + graph->SetMarkerSize(0.03); + graph->Draw("AP"); + graph->GetXaxis()->SetLabelSize(0.03); + graph->GetXaxis()->SetTitle("Time [MJD]"); + // lightcurvegraph->GetXaxis()->SetTimeDisplay(1); + + +}