Changeset 3760 for trunk/MagicSoft/Mars/macros

Timestamp:

04/15/04 19:54:59 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/macros/pedphotcalc.C (modified) (17 diffs)

trunk/MagicSoft/Mars/macros/pedphotcalc.C

@@ -26 +26 @@
 !
 \* ======================================================================== */
-
+/////////////////////////////////////////////////////////////////////////////
+//
+//  pedphotcalc.C
+//
+//  This macro is an example of the use of MPedPhotCalc. It computes
+//  the pedestal mean and rms from pedestal files undergoing the
+//  signal extraction + calibration chain, in units of photons. It
+//  produces plots of the relevant computed quantities.
+//
+//  Needs as arguments the run number of a pedestal file ("*_P_*.root"), 
+//  one of a calibration file ("*_C_*.root").
+//  performs the pedestal calculation, the calibration 
+/// constants calculation and the calibration of the pedestals. 
+//
+//  The TString inpath has to be set correctly.
+//
+//  The macro searches for the pulser colour which corresponds to the calibration
+//  run number. If the run number is smaller than 20000, pulser colour "CT1" 
+//  is assumed, otherwise, it searches for the strings "green", "blue", "uv" or 
+//  "ct1" in the filenames. If no colour or multiple colours are found, the 
+//  execution is aborted.  
+//
+////////////////////////////////////////////////////////////////////////////////////
 #include "MParList.h"
 #include "MTaskList.h"
+#include "MJPedestal.h"
+#include "MJCalibration.h"
 #include "MPedestalCam.h"
+#include "MPedestalPix.h"
 #include "MReadMarsFile.h"
 #include "MGeomApply.h"
-#include "MPedCalcPedRun.h"
 #include "MGeomCamMagic.h"
 #include "MEvtLoop.h"
+#include "MCalibrationCam.h"
 #include "MCalibrationChargeCam.h"
-#include "MHCalibrationChargeCam.h"
-#include "MCalibrationChargePINDiode.h"
-#include "MHCalibrationChargePINDiode.h"
-#include "MCalibrationChargeBlindPix.h"
-#include "MHCalibrationChargeBlindPix.h"
-#include "MCalibrationChargeCalc.h"
+#include "MCalibrationChargePix.h"
+#include "MCalibrationQECam.h"
+#include "MCalibrationQEPix.h"
 #include "MExtractedSignalCam.h"
 #include "MExtractSignal.h" 
-#include "MExtractedSignalPINDiode.h"
-#include "MExtractPINDiode.h" 
-#include "MExtractedSignalBlindPixel.h"
-#include "MExtractBlindPixel.h" 
 #include "MCerPhotEvt.h"
 #include "MCalibrate.h"
 #include "MPedPhotCalc.h"
 #include "MPedPhotCam.h"
-#include "MBadPixelsCam.h"
+#include "MPedPhotPix.h"
+#include "MHCamera.h"
+#include "MRunIter.h"
+#include "MDirIter.h"
+#include "MStatusDisplay.h"
 #include "MHCamera.h"
 
@@ -58 +80 @@
 #include "TString.h"
 #include "TCanvas.h"
-
-#include <iostream>
-
-
-const TString pedfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
-const TString calfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root"; 
-
-void pedphotcalc(TString pedname=pedfile,TString calname=calfile)
+#include "TStyle.h"
+#include "TF1.h"
+#include "TLegend.h"
+
+#include <iostream.h>
+#include "Getline.h"
+
+const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/";
+const Int_t dpedrun  = 14607;
+const Int_t dcalrun1 = 14608;
+const Int_t dcalrun2 = 0;
+const Bool_t usedisplay = kTRUE;
+
+static MCalibrationCam::PulserColor_t FindColor(MDirIter* run) 
 {
   
-  /*
-    This macro is an example of the use of MPedPhotCalc. It computes
-    the pedestal mean and rms from pedestal files undergoing the
-    signal extraction + calibration chain, in units of photons. It
-    produces plots of the relevant computed quantities.
-  */
-  
+  MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE;
+
+  TString filenames;
+
+  while (!(filenames=run->Next()).IsNull())
+    {
+
+      filenames.ToLower();
+
+      if (filenames.Contains("green"))
+        if (col == MCalibrationCam::kNONE)
+          {
+            cout << "Found colour: Green  in " << filenames << endl;
+            col = MCalibrationCam::kGREEN;
+          }
+        else if (col != MCalibrationCam::kGREEN)
+          {
+            cout << "Different colour found in " << filenames << "... abort" << endl;
+            return MCalibrationCam::kNONE;
+          }
+
+      if (filenames.Contains("blue"))
+        if (col == MCalibrationCam::kNONE)
+          {
+            cout << "Found colour: Blue  in " << filenames << endl;
+            col = MCalibrationCam::kBLUE;
+          }
+        else if (col != MCalibrationCam::kBLUE)
+          {
+            cout << "Different colour found in " << filenames << "... abort" << endl;
+            return MCalibrationCam::kNONE;
+          }
+
+      if (filenames.Contains("uv"))
+        if (col == MCalibrationCam::kNONE)
+          {
+            cout << "Found colour: Uv  in " << filenames << endl;
+            col = MCalibrationCam::kUV;
+          }
+        else if (col != MCalibrationCam::kUV)
+          {
+            cout << "Different colour found in " << filenames << "... abort" << endl;
+            return MCalibrationCam::kNONE;
+          }
+
+      if (filenames.Contains("ct1"))
+        if (col == MCalibrationCam::kNONE)
+          {
+            cout << "Found colour: Ct1  in " << filenames << endl;
+            col = MCalibrationCam::kCT1;
+          }
+        else if (col != MCalibrationCam::kCT1)
+          {
+            cout << "Different colour found in " << filenames << "... abort" << endl;
+            return MCalibrationCam::kNONE;
+          }
+      
+    }
+  
+
+      
+  if (col == MCalibrationCam::kNONE)
+    cout <<  "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString() 
+         << "... abort" << endl;
+  
+  return col;      
+}
+
+
+
+void DrawProjection(MHCamera *obj1, Int_t fit) 
+{
+    TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
+    obj2->SetDirectory(0);
+    obj2->Draw();
+    obj2->SetBit(kCanDelete);
+    gPad->SetLogy();
+
+    if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
+    {
+        TArrayI s0(3);
+        s0[0] = 6;
+        s0[1] = 1;
+        s0[2] = 2;
+
+        TArrayI s1(3);
+        s1[0] = 3;
+        s1[1] = 4;
+        s1[2] = 5;
+
+        TArrayI inner(1);
+        inner[0] = 0;
+
+        TArrayI outer(1);
+        outer[0] = 1;
+
+        // Just to get the right (maximum) binning
+        TH1D *half[4];
+        half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
+        half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
+        half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
+        half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
+
+        for (int i=0; i<4; i++)
+        {
+            half[i]->SetLineColor(kRed+i);
+            half[i]->SetDirectory(0);
+            half[i]->SetBit(kCanDelete);
+            half[i]->Draw("same");
+        }
+    }
+
+    const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
+    const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
+    const Double_t integ = obj2->Integral("width")/2.5066283;
+    const Double_t mean  = obj2->GetMean();
+    const Double_t rms   = obj2->GetRMS();
+    const Double_t width = max-min;
+
+    if (rms==0 || width==0)
+        return;
+
+    const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
+                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])");
+
+    TF1 *f=0;
+    switch (fit)
+    {
+        // FIXME: MAYBE add function to TH1?
+    case 0:
+        f = new TF1("sgaus", "gaus(0)", min, max);
+        f->SetLineColor(kYellow);
+        f->SetBit(kCanDelete);
+        f->SetParNames("Area", "#mu", "#sigma");
+        f->SetParameters(integ/rms, mean, rms);
+        f->SetParLimits(0, 0,   integ);
+        f->SetParLimits(1, min, max);
+        f->SetParLimits(2, 0,   width/1.5);
+        obj2->Fit(f, "QLRM");
+        break;
+
+    case 1:
+        f = new TF1("dgaus", dgausformula, min, max);
+        f->SetLineColor(kYellow);
+        f->SetBit(kCanDelete);
+        f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
+        f->SetParameters(integ,         (min+mean)/2, width/4,
+                         integ/width/2, (max+mean)/2, width/4);
+        // The left-sided Gauss
+        f->SetParLimits(0, integ-1.5,      integ+1.5);
+        f->SetParLimits(1, min+(width/10), mean);
+        f->SetParLimits(2, 0,              width/2);
+        // The right-sided Gauss
+        f->SetParLimits(3, 0,    integ);
+        f->SetParLimits(4, mean, max-(width/10));
+        f->SetParLimits(5, 0,    width/2);
+        obj2->Fit(f, "QLRM");
+        break;
+
+    default:
+        obj2->Fit("gaus", "Q");
+        obj2->GetFunction("gaus")->SetLineColor(kYellow);
+        break;
+    }
+}
+
+void CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
+{
+    c.cd(x);
+    gPad->SetBorderMode(0);
+    MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
+
+    c.cd(x+y);
+    gPad->SetBorderMode(0);
+    obj1->Draw();
+
+    c.cd(x+2*y);
+    gPad->SetBorderMode(0);
+    DrawProjection(obj1, fit);
+}
+
+void pedphotcalc(const Int_t prun=dpedrun, // pedestal file
+                 const Int_t crun1=dcalrun1, const Int_t crun2=dcalrun2 // calibration file(s)
+                 )
+{
+  
+  MRunIter pruns;
+  MRunIter cruns;
+
+  pruns.AddRun(prun,inpath);
+
+  if (crun2==0)
+    cruns.AddRun(crun1,inpath);
+  else
+    cruns.AddRuns(crun1,crun2,inpath);
+
+  MCalibrationCam::PulserColor_t color;
+
+  if (crun1 < 20000)
+    color = MCalibrationCam::kCT1;
+  else
+    color = FindColor((MDirIter*)&cruns);
+
+  if (color == MCalibrationCam::kNONE)
+    {
+      TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
+      
+      while (1)
+        {
+          timer.TurnOn();
+          TString input = Getline("Could not find the correct colour: Type 'q' to exit, "
+                                  "green, blue, uv or ct1 to go on: ");
+          timer.TurnOff();
+          
+          if (input=="q\n")
+            return ;
+          
+          if (input=="green")
+            color = MCalibrationCam::kGREEN;
+          if (input=="blue")
+            color = MCalibrationCam::kBLUE;
+          if (input=="uv")
+            color = MCalibrationCam::kUV;
+          if (input=="ct1")
+            color = MCalibrationCam::kCT1;
+        }
+    }
+
+  //
+  // Now setup the tasks and tasklist for the pedestals:
+  // ---------------------------------------------------
+  //
+  MGeomCamMagic     geomcam;
+  MGeomApply        geomapl;
+  MStatusDisplay   *display = NULL;
+
   /************************************/
   /* FIRST LOOP: PEDESTAL COMPUTATION */
   /************************************/
   
-  MParList plist;
-  MTaskList tlist;
-  plist.AddToList(&tlist);
-  
-  // containers
-  MPedestalCam   pedcam;
-  plist.AddToList(&pedcam);
-    
-  //tasks
-  MReadMarsFile read("Events", pedname);
-  MGeomApply     geomapl;
-  MPedCalcPedRun pedcalc;
-  MGeomCamMagic  geomcam;
-
-  read.DisableAutoScheme();
-    
-  tlist.AddToList(&read);
-  tlist.AddToList(&geomapl);
-  tlist.AddToList(&pedcalc);
-  
-  // Create and execute the event looper
-  MEvtLoop evtloop;
-  evtloop.SetParList(&plist);
+  MJPedestal pedloop;
+  pedloop.SetInput(&pruns);
+  if (usedisplay)
+    {
+      display = new MStatusDisplay;
+      display->SetUpdateTime(3000);
+      display->Resize(850,700);
+      display->SetBit(kCanDelete);
+      pedloop.SetDisplay(display);
+    }
   
   cout << "*************************" << endl;
   cout << "** COMPUTING PEDESTALS **" << endl;
   cout << "*************************" << endl;
-  if (!evtloop.Eventloop())
+
+  if (!pedloop.Process())
     return;
-  
-  tlist.PrintStatistics();
+
+  MPedestalCam pedcam = pedloop.GetPedestalCam();
 
   /****************************/
@@ -115 +360 @@
   /****************************/
   
-  MParList  plist2;
-  MTaskList tlist2;
-  plist2.AddToList(&tlist2);
-  
-  // containers  
-  MCalibrationChargeCam calcam;
-  MExtractedSignalCam   sigcam;
-  MBadPixelsCam         badcam;
-
-  MCalibrationChargeCam        chargecam;
-  MCalibrationChargePINDiode   chargepin;
-  MCalibrationChargeBlindPixel chargeblind;
-  chargeblind.SetColor(MCalibrationChargeBlindPix::kECT1);
-  chargepin.SetColor(  MCalibrationChargePINDiode::kECT1);
-  //
-  plist2.AddToList(&geomcam);
-  plist2.AddToList(&pedcam);
-  plist2.AddToList(&calcam);
-  plist2.AddToList(&sigcam);
-  plist2.AddToList(&badcam);  
-  plist2.AddToList(&chargecam);
-  plist2.AddToList(&chargepin);  
-  plist2.AddToList(&chargeblind);  
-
-  //tasks
-  MReadMarsFile read2("Events", calname);
-  read2.DisableAutoScheme();
-  
-  MExtractPINDiode       pincalc;
-  MExtractBlindPixel     blindcalc;
-  MExtractSignal         sigcalc;
-  MCalibrationChargeCalc calcalc;
-  
-  MFillH fillpin("MHCalibrationChargePINDiode",     "MExtractedSignalPINDiode");
-  MFillH fillblind("MHCalibrationChargeBlindPixel", "MExtractedSignalBlindPixel");
-  MFillH fillcam("MHCalibrationChargeCam"     ,     "MExtractedSignalCam");
-  // 
-  // Apply a filter against cosmics
-  // (was directly in MCalibrationChargeCalc in earlier versions)
-  //
-  MFCosmics            cosmics;
-  MContinue            cont(&cosmics);
-
-  tlist2.AddToList(&read2);
-  tlist2.AddToList(&geomapl);
-  tlist2.AddToList(&blindcalc);
-  tlist2.AddToList(&pincalc);
-  tlist2.AddToList(&sigcalc);
-  //
-  // In case, you want to skip the cosmics rejection, 
-  // uncomment the next line
-  //
-  tlist2.AddToList(&cont);
-  //
-  tlist2.AddToList(&fillpin);
-  tlist2.AddToList(&fillblind);
-  tlist2.AddToList(&fillcam);
-  tlist2.AddToList(&calcalc);
-  
-  // Execute second loop (calibration)
-  MEvtLoop evtloop2;
-  evtloop2.SetParList(&plist2);
+  //
+  // Now setup the new tasks for the calibration:
+  // ---------------------------------------------------
+  //
+  MJCalibration     calloop;
+  calloop.SetColor(color);
+  calloop.SetInput(&cruns);
+  //
+  // Use as signal extractor MExtractSignal:
+  //
+  calloop.SetExtractorLevel(1);
+  //
+  // The next two commands are for the display:
+  //
+  if (usedisplay)
+    {
+      calloop.SetDisplay(display);
+      calloop.SetDataCheck();
+    }
   
   cout << "***********************************" << endl;
   cout << "** COMPUTING CALIBRATION FACTORS **" << endl;
   cout << "***********************************" << endl;
-  if (!evtloop2.Eventloop())
+
+  if (!calloop.Process(pedcam))
     return;
   
-  tlist2.PrintStatistics();
-
-  
+  MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
+  MCalibrationQECam     &qecam  = calloop.GetQECam();
+
   /************************************************************************/
   /* THIRD LOOP: PEDESTAL COMPUTATION USING EXTRACTED SIGNAL (IN PHOTONS) */
@@ -202 +405 @@
   
   plist3.AddToList(&geomcam);
-  plist3.AddToList(&pedcam);
-  plist3.AddToList(&calcam);
+  plist3.AddToList(&pedcam );
+  plist3.AddToList(&calcam );
+  plist3.AddToList(&qecam  );
   plist3.AddToList(&photcam);
   plist3.AddToList(&extedsig);
@@ -209 +413 @@
   
   //tasks
-  MReadMarsFile read3("Events", pedname);
+  MReadMarsFile read3("Events");
+  read3.DisableAutoScheme();
+  static_cast<MRead&>(read3).AddFiles(pruns);  
+
+  MExtractSignal  sigcalc;
   MCalibrate      photcalc;
+  photcalc.SetCalibrationMode(MCalibrate::kFfactor);
   MPedPhotCalc    pedphotcalc;  
 
-  read3.DisableAutoScheme();
-  
   tlist3.AddToList(&read3);
   tlist3.AddToList(&geomapl);
@@ -233 +440 @@
   tlist3.PrintStatistics();
 
-
   /**********************/
   /* PRODUCE NICE PLOTS */
   /**********************/
 
+  if (usedisplay)
+    {
+      
+      MHCamera disp0(geomcam, "MPedPhotCam;ped", "Pedestals");
+      MHCamera disp1(geomcam, "MPedPhotCam;rms", "Sigma Pedestal");
+      
+      disp0.SetCamContent(pedphotcam, 0);
+      disp0.SetCamError  (pedphotcam, 1);
+      
+      disp1.SetCamContent(pedphotcam, 1);
+      
+      disp0.SetYTitle("Pedestal signal [photons]");
+      disp1.SetYTitle("Pedestal signal RMS [photons]");
+      
+      //
+      // Display data
+      //
+      TCanvas &c1 = display->AddTab("PedPhotCam");
+      c1.Divide(2,3);
+      
+      CamDraw(c1, 1, 2, disp0, 0);
+      CamDraw(c1, 2, 2, disp1, 1);
+      
+    }
+  
+
+#if 0
   const UShort_t npix = 577;
 
@@ -250 +483 @@
   TH1F* calhist    = new TH1F("calhist","Calibration factors",npix,0,npix);
   TH1F* caldist    = new TH1F("caldist","Calibration factors",100,0,1);
+  TH1F* qehist     = new TH1F("qehist", "Quantrum efficiencies",npix,0,npix);
+  TH1F* qedist     = new TH1F("qedist", "Quantrum efficiencies",100,0,1);
 
   // pedestal signals
@@ -260 +495 @@
   for(int i=0;i<npix;i++)
     {
+      MCalibrationChargePix &calpix = (MCalibrationChargePix&)calcam[i];
+      MCalibrationQEPix     &qepix  = (MCalibrationQEPix&)    qecam[i];
+
       const Float_t ped        = pedcam[i].GetPedestal();
       const Float_t pedrms     = pedcam[i].GetPedestalRms();
-      const Float_t cal        = calcam[i].GetMeanConversionBlindPixelMethod();
-      const Float_t calerr     = calcam[i].GetErrorConversionBlindPixelMethod();
+      const Float_t cal        = calpix.GetMeanConvFADC2Phe();
+      const Float_t qe         = qepix .GetQECascadesFFactor();
       const Float_t pedphot    = pedphotcam[i].GetMean();
       const Float_t pedphotrms = pedphotcam[i].GetRms();
@@ -274 +512 @@
       calhist->Fill(i,cal);
       caldist->Fill(cal);
+      qehist->Fill(i,qe);
+      qedist->Fill(qe);
 
       pedphothist->Fill(i,pedphot);
@@ -287 +527 @@
   canvas->SetBorderMode(0);    // Delete the Canvas' border line     
   canvas->cd();
-  canvas->SetBorderMode(0);
-     
-  canvas->Divide(2,4);
+
+  canvas->Divide(2,5);
 
   // draw pedestal histo
@@ -322 +561 @@
   pedrmsdist->Draw("same");
 
-
   TLegend* leg2 = new TLegend(.14,.68,.39,.88);
   leg2->SetHeader("");
@@ -339 +577 @@
   calhist->SetMaximum(1);
   calhist->SetMinimum(0);
-  calhist->GetYaxis()->SetTitle("Calibration factor (photon/ADC count)");
+  calhist->GetYaxis()->SetTitle("Calibration factor (phe/ADC count)");
   calhist->SetStats(kFALSE);
   calhist->Draw();
@@ -347 +585 @@
   gPad->cd();
   gPad->SetBorderMode(0);
-  caldist->GetXaxis()->SetTitle("Calibration factor (photon/ADC count)");
+  caldist->GetXaxis()->SetTitle("Calibration factor (phe/ADC count)");
   caldist->Draw();
 
+  // draw qe histo
+  canvas->cd(5);
+  gPad->cd();
+  gPad->SetBorderMode(0);
+  qehist->GetXaxis()->SetTitle("Pixel SW Id");
+  qehist->SetMaximum(1);
+  qehist->SetMinimum(0);
+  qehist->GetYaxis()->SetTitle("Quantum efficiency for cascades");
+  qehist->SetStats(kFALSE);
+  qehist->Draw();
+
+  // draw qe distribution
+  canvas->cd(6);
+  gPad->cd();
+  gPad->SetBorderMode(0);
+  qedist->GetXaxis()->SetTitle("Quantum efficiency for cascades");
+  qedist->Draw();
+
   // draw pedestal signal histo
-  canvas->cd(5);
+  canvas->cd(7);
   gPad->cd();
   gPad->SetBorderMode(0);
@@ -360 +616 @@
 
   // draw pedestal signal distribution
-  canvas->cd(6);
+  canvas->cd(8);
   gPad->cd();
   gPad->SetBorderMode(0);
@@ -367 +623 @@
 
   // draw pedestal signal rms histo
-  canvas->cd(7);
+  canvas->cd(9);
   gPad->cd();
   gPad->SetBorderMode(0);
@@ -376 +632 @@
 
   // draw pedestal signal rms distribution
-  canvas->cd(8);
+  canvas->cd(10);
   gPad->cd();
   gPad->SetBorderMode(0);
@@ -382 +638 @@
   pedphotrmsdist->Draw();
 
-  return; 
+  canvas->SaveAs("pedphotcalc.root");
+
+#endif
 }
-