- Timestamp:
- 07/20/04 11:57:23 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mtemp/mifae
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc
r4407 r4408 148 148 149 149 fLiveTime->SetRealTime(mean,width); 150 *fLog << inf << GetName() << ": New time bin" << endl; 151 Print("all"); 152 *fLog << GetName() << ": First event time " << setprecision(10) << fFirstEventMjd << setprecision(5) << endl; 153 fLiveTime->Print("last"); 150 151 if (Debug) 152 { 153 *fLog << inf << GetName() << ": New time bin" << endl; 154 Print("all"); 155 *fLog << GetName() << ": First event time " << setprecision(10) << fFirstEventMjd << setprecision(5) << endl; 156 fLiveTime->Print("last"); 157 } 154 158 155 159 fLiveTime->AddBin(); -
trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C
r4397 r4408 67 67 // Configuration 68 68 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"; 71 71 72 72 //Constanst 73 73 const Double_t kConvDegToRad = TMath::Pi()/180.; 74 const Double_t kSec = 1e 6; //[sec/microsec]74 const Double_t kSec = 1e3; //[sec/milisec] 75 75 const Double_t kDay = 24.*60.*60.; //[Day/sec] 76 76 77 77 UInt_t numberTimeBins = 1; 78 78 TArrayI numberOnEvents(1); 79 TArrayD numberOnEventsAfterCleaning(1); 79 80 TArrayD numberBkgEventsToNormOn(1); 80 81 TArrayD timeOn(1); 81 82 82 83 TArrayD meanTimeBinOn(1); 83 84 TArrayD widthTimeBinOn(1); … … 85 86 TArrayD zenithMinimumOn(1); 86 87 TArrayD zenithMaximumOn(1); 88 89 TArrayD deadFractionOn(1); 87 90 88 91 TObjArray coszenithHistoOn; 89 92 TObjArray alphaHistoOn; 90 93 TObjArray srcposHistoOn; 91 92 94 TObjArray timediffHistoOn; 93 95 94 96 const Int_t numberSizeBins = 4; … … 101 103 Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 }; 102 104 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 }; 105 107 106 108 //alpha plot integration for gammas … … 170 172 TH1F *hAlpha_on = new TH1F ("hAlphaOn","",nbins,minalpha,maxalpha); 171 173 172 Int_t nbins_srcpos = 400;173 Double_t minsrcpos = - 600.;174 Double_t maxsrcpos = 600.; //[mm] //!!! position precition 3mm ~ 0.01 deg174 Int_t nbins_srcpos = 200; 175 Double_t minsrcpos = -200.; 176 Double_t maxsrcpos = 200.; //[mm] //!!! position precition 3mm ~ 0.01 deg 175 177 TH2F *hSrcPos_on = new TH2F ("hSrcPosOn","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos); 176 178 … … 224 226 Bool_t flag = kFALSE; 225 227 228 numberOnEventsAfterCleaning[numberTimeBins-1] = 0; 226 229 zenithMinimumOn[numberTimeBins-1] = 100.; 227 230 zenithMaximumOn[numberTimeBins-1] = 0.; … … 241 244 TString coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1); 242 245 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); 243 253 244 254 while(tlist_on.Process()) … … 248 258 UInt_t evt = evtheader_on.GetDAQEvtNumber(); 249 259 Double_t mjd = time_on.GetMjd(); 250 260 261 numberOnEventsAfterCleaning[numberTimeBins-1]++; 262 251 263 if (mjd == 0) 252 253 264 { 265 254 266 if (debug) 255 267 { … … 362 374 mjdFirstEvent = mjd; 363 375 } 376 else if (mjdFirstEventofBin != 0) 377 hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec); 364 378 365 379 if (mjdFirstEventofBin == 0) … … 368 382 mjdFirstEventofBin = mjd; 369 383 } 370 384 371 385 evtLastEvent = evt; 372 386 runLastEvent = run; 373 387 mjdLastEvent = mjd; 374 388 mtimeLastEvent = time_on; 375 376 377 389 378 390 Double_t size = hillas_on.GetSize(); 379 391 Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg(); 380 392 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(); 382 397 Double_t alpha = hillassrc_on.GetAlpha(); 383 398 Double_t srcposx = srcpos_on.GetX(); … … 394 409 if (length > lengthmin[sizebin] && length < lengthmax[sizebin]) 395 410 { 396 if (dist > distmin[sizebin] && dist < distmax[sizebin])411 if (dist > distmin[sizebin] && centerdist < distmax[sizebin]) 397 412 { 398 413 … … 439 454 mjdFirstEvent = mjd; 440 455 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 442 462 alphaHistoOn.AddLast(hAlpha_on_abs_timebin); 443 463 srcposHistoOn.AddLast(hSrcPos_on_timebin); 444 464 coszenithHistoOn.AddLast(hCosZenith_on_timebin); 465 timediffHistoOn.AddLast(hTimeDiff_on_timebin); 445 466 446 467 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl; … … 456 477 zenithMinimumOn.Set(numberTimeBins); 457 478 zenithMaximumOn.Set(numberTimeBins); 458 479 deadFractionOn.Set(numberTimeBins); 480 numberOnEventsAfterCleaning.Set(numberTimeBins); 481 459 482 timeOn[numberTimeBins-1] = 0.; 460 483 zenithMinimumOn[numberTimeBins-1] = 100.; 461 484 zenithMaximumOn[numberTimeBins-1] = 0.; 485 numberOnEventsAfterCleaning[numberTimeBins-1] = 0; 462 486 463 487 alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1); … … 470 494 hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith); 471 495 496 timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1); 497 hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff); 498 472 499 flag = kTRUE; 473 500 } 474 } 501 502 } 503 475 504 } 476 505 … … 484 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); 485 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 486 522 cout.precision(10); 487 523 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl; … … 491 527 srcposHistoOn.AddLast(hSrcPos_on_timebin); 492 528 coszenithHistoOn.AddLast(hCosZenith_on_timebin); 529 timediffHistoOn.AddLast(hTimeDiff_on_timebin); 493 530 494 531 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl; … … 499 536 500 537 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 503 542 504 543 //-----------------------OFF------------------------ … … 616 655 Double_t width = hillas_off.GetWidth()*geomcam.GetConvMm2Deg(); 617 656 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); 618 660 619 661 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins); … … 625 667 if (length > lengthmin[sizebin] && length < lengthmax[sizebin]) 626 668 { 627 //General histos 628 629 if (dist > distmin[sizebin] && dist < distmax[sizebin]) 669 if (centerdist < distmax[sizebin]) 630 670 { 631 hSize_off->Fill(log10(size)); 632 hWidth_off->Fill(width); 633 hLength_off->Fill(length); 671 //General histos 634 672 635 673 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 histos651 for (UInt_t bin=0; bin<numberTimeBins; bin++)652 {653 srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);654 674 srcplace_timebin.CallProcess(); 655 675 csrc_off.CallProcess(); … … 660 680 Double_t srcposy = srcpos_off.GetY(); 661 681 662 if (dist > distmin[sizebin] && dist < distmax[sizebin])682 if (dist > distmin[sizebin] ) 663 683 { 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 } 666 711 } 667 712 } … … 684 729 TArrayD widthTimeBinOnInSec(numberTimeBins); 685 730 731 TArrayD meanTriggerRateOn(numberTimeBins); 732 TArrayD errorMeanTriggerRateOn(numberTimeBins); 733 686 734 for (UInt_t bin=0; bin<numberTimeBins; bin++) 687 735 { 688 736 737 meanTriggerRateOn[bin] = numberOnEventsAfterCleaning[bin]/timeOn[bin]; 738 errorMeanTriggerRateOn[bin] = TMath::Sqrt(numberOnEventsAfterCleaning[bin])/timeOn[bin]; 739 689 740 meanTimeBinOnInSec[bin] = (meanTimeBinOn[bin]-(Int_t)meanTimeBinOn[bin]); 690 741 if (meanTimeBinOnInSec[bin] > 0.5) … … 701 752 numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin]; 702 753 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.); 704 755 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.); 706 757 errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin]; 707 758 } … … 718 769 cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl; 719 770 cout << " numberExcessEvents\t " << numberExcessEvents[bin] << " +- " << errorNumberExcessEvents[bin] << endl; 771 cout << " deadFraction\t" << deadFractionOn[bin] << endl; 720 772 cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl; 773 cout << " Trigger Rate\t " << meanTriggerRateOn[bin] << " +- " << errorMeanTriggerRateOn[bin] << endl; 721 774 } 722 775 … … 739 792 c0->Print(openpsname); 740 793 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 741 807 TCanvas** c = new TCanvas*[numberTimeBins]; 742 808 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 743 820 for (UInt_t bin=0; bin<numberTimeBins-1; bin++) 744 821 { … … 758 835 ((TH1F*)alphaHistoOff[bin])->SetLineColor(46); 759 836 ((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); 760 838 ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.); 761 839 alphaTitle = Form("%s%02i","hAlphaOnOff",bin); … … 786 864 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46); 787 865 ((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); 788 867 ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.); 789 868 alphaTitle = Form("%s%02i","hAlphaOnOff",numberTimeBins-1); … … 798 877 c[numberTimeBins-1]->Print(psname); 799 878 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(); 812 894 813 895 // ############################################################################
Note:
See TracChangeset
for help on using the changeset viewer.