Changeset 7188 for trunk/MagicSoft/Mars/macros

Timestamp:

07/13/05 19:06:26 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/macros

Files:

• starmc.C (modified) (7 diffs)
• starmcstereo.C (modified) (9 diffs)

trunk/MagicSoft/Mars/macros/starmc.C

@@ -52 +52 @@
   // differences in gain of outer pixels)
   //
-  CalibrationFilename = new TString("/data1/magic/mc_data/root/Period025/gammas_nonoise/Gamma_zbin0_*.root");
+  CalibrationFilename = new TString("nonoise/Gamma_zbin0_0_7_1000to1009_w0.root");
   // File to be used in the calibration (must be a camera file without added noise)
 
-
-  Char_t* AnalysisFilename = "Gamma_zbin1_0_*.root";  // File to be analyzed
-
+  Char_t* AnalysisFilename = "Gamma_zbin0_0_*.root";  // File to be analyzed
 
 
@@ -76 +74 @@
   // USER CHANGE: tail cuts for image analysis
 
-  Float_t CleanLev[2] = {7., 5.}; 
-  MImgCleanStd  clean(CleanLev[0], CleanLev[1]); // Applies tail cuts to image.
-  clean.SetMethod(MImgCleanStd::kAbsolute);      // In absolute units (phot or phe- as chosen below)
-
+  Float_t CleanLev[2] = {10., 5.}; // Tail cuts for the analysis loop
+  MImgCleanStd clean2(CleanLev[0], CleanLev[1]); // Applies tail cuts to image.
+  clean2.SetMethod(MImgCleanStd::kAbsolute); 
 
   //  USER CHANGE: signal extraction
-  //
-  //  MExtractFixedWindowPeakSearch sigextract;
-  //  sigextract.SetWindows(6, 6, 4);
 
   //  MExtractTimeAndChargeDigitalFilter sigextract;
@@ -91 +85 @@
 
   MExtractTimeAndChargeSpline sigextract;
+  sigextract.SetChargeType(MExtractTimeAndChargeSpline::kIntegral);
   sigextract.SetRiseTimeHiGain(0.5);
   sigextract.SetFallTimeHiGain(0.5);
@@ -96 +91 @@
   // USER CHANGE: high to low gain ratio. DEPENDS on EXTRACTOR!!
   // One should calculate somewhere else what this factor is for each extractor!
-  Float_t hi2low = 12.; // tentative value for spline with risetime 0.5, fall time 0.5
+
+  //  Float_t hi2low = 10.83; 
+  // value for digital filter, HG window 4, LG window 6
+
+  Float_t hi2low = 11.28; 
+  // value for spline with risetime 0.5, fall time 0.5, low gain stretch 1.5
 
 
@@ -108 +108 @@
 
 
+  MImgCleanStd  clean(0.,0.); 
+  // All pixels containing any photon will be accepted. This is what we want
+  // for the calibration loop (since no noise is added)
+  clean.SetMethod(MImgCleanStd::kAbsolute); 
+  // In absolute units (phot or phe- as chosen below)
+
   MSrcPosCam src;
   //
-  // ONLY FOR WOBBLE MODE!!
-  //   src.SetX(120.);  // units: mm
+  // ONLY FOR WOBBLE MODE!! : set the rigt source position on camera!
+  //   src.SetX(120.);   // units: mm. Value for MC w+ files
+  //   src.SetX(-120.);  // units: mm. Value for MC w- files
 
   src.SetReadyToSave();
@@ -244 +251 @@
   //
 
+
   MProgressBar bar;
   bar.SetWindowName("Calibrating...");
@@ -273 +281 @@
   tlist.RemoveFromList(&read);
 
+  tlist.AddToListBefore(&clean2, &clean);
+  tlist.RemoveFromList(&clean);
+
   //
   // Analyzed only the desired fraction of events, skip the rest:

trunk/MagicSoft/Mars/macros/starmcstereo.C

@@ -17 +17 @@
 !
 !   Author(s):  Abelardo Moralejo 1/2004 <mailto:moralejo@pd.infn.it>
-!               Thomas Bretz, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -33 +32 @@
 /////////////////////////////////////////////////////////////////////////////
 
-//
-// User change. 
-//
-
-Float_t ctx[7] = {0., 0., 0., 0., 0., 0., 0.}; 
-Float_t cty[7] = {-70., -40., -30., 30., 50., 60., 70.}; // in meters
-//
-// FIXME: unfortunately present version of reflector was not prepared for
-// stereo configurations and keeps no track of CT position. So the positions
-// must be set above by the user, making sure that they correspond to the
-// files one is analysing.
-//
-
-void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 5)
+void starmcstereo(Int_t ct1 = 1, Int_t ct2 = 2)
 {
-  if (ct1 > sizeof(ctx)/sizeof(ctx[0]) || 
-      ct2 > sizeof(ctx)/sizeof(ctx[0]) )
-    {
-      cout << endl << "Wrong CT id number!" << endl;
+  // ------------- user change -----------------
+
+  TString* CalibrationFilename = 0;
+
+  // Calibration file: a file with no added noise. Comment out next line if you
+  // do not want to calibrate the data (means SIZE will be in ADC counts)
+
+  CalibrationFilename = new TString("nonoise/Gamma_20_0_7_200000to200009_XX_w0.root");
+
+  Char_t* AnalysisFilename = "Gamma_20_0_7_*_XX_w0.root";  // File to be analyzed
+  Char_t* OutFileTag      = "gammas";           // Output file tag
+
+  // First open input files to check that the required telescopes 
+  // are in the file, and get telescope coordinates.
+
+  TChain *rh = new TChain("RunHeaders");
+  rh->Add(AnalysisFilename);
+  MMcCorsikaRunHeader *corsrh = new MMcCorsikaRunHeader();
+  rh->SetBranchAddress("MMcCorsikaRunHeader.", &corsrh);
+  rh->GetEvent(0); 
+  // We assume that all the read files will have the same telescopes inside, 
+  // so we look only into the first runheader.
+  Int_t allcts = corsrh->GetNumCT();
+  if (ct1 >  allcts || ct2 >  allcts)
+    {
+      cout << endl << "Wrong CT id number, not contained in input file!" << endl;
       return;
     }
+  // Set telescope coordinates as read from first runheader:
+  Float_t ctx[2];
+  Float_t cty[2];
+  ctx[0] = ((*corsrh)[ct1-1])->GetCTx();
+  cty[0] = ((*corsrh)[ct1-1])->GetCTy();
+  ctx[1] = ((*corsrh)[ct2-1])->GetCTx();
+  cty[1] = ((*corsrh)[ct2-1])->GetCTy();
+
+  // Now find out number of pixels in each camera:
+  MMcConfigRunHeader* confrh1 = new MMcConfigRunHeader();
+  MMcConfigRunHeader* confrh2 = new MMcConfigRunHeader();
+  rh->SetBranchAddress("MMcConfigRunHeader;1.", &confrh1);
+  rh->SetBranchAddress("MMcConfigRunHeader;2.", &confrh2);
+  rh->GetEvent(0);
+  Int_t npix[2];
+  npix[0] = confrh1->GetNumPMTs();
+  npix[1] = confrh2->GetNumPMTs();
+
+  rh->Delete();
+
 
   Int_t CT[2] = {ct1, ct2};  // Only 2-telescope analysis for the moment
-  Int_t NCTs = sizeof(CT)/sizeof(CT[0]);
+  Int_t NCTs = 2;
 
 
   // ------------- user change -----------------
 
-  Char_t* AnalysisFilename = "gam-yXX-00001.root";  // File to be analyzed
-  Char_t* OutFileTag      = "gammas";           // Output file tag
-
-  Float_t CleanLev[2] = {4., 3.}; // Tail cuts for image analysis
-
-  Int_t BinsHigh[2] = {5, 9}; // First and last FADC bin of the range to be integrated,
-  Int_t BinsLow[2]  = {5, 9}; // for high and low gain respectively.
+  Float_t BinsHigh[2] = {0, 79};
+  Float_t BinsLow[2]  = {0, 79};  // FADC slices (2GHz sampling)
+  Float_t CleanLev[2] = {7., 5.}; // Units: phes (absolute cleaning will be used later!) 
+  // Tail cuts for the analysis loop. In the first (calibration) loop they will not 
+  // be used; we run over a noiseless file and we want to accept all pixels with any 
+  // number of phes.
+
+
+  MImgCleanStd**     clean = new MImgCleanStd*[NCTs];  
+
+  MImgCleanStd* clean[0] = new MImgCleanStd(1.,1.);
+  MImgCleanStd* clean[1] = new MImgCleanStd(1.,1.); 
+  // Just dummy levels. Since the calibration file will be a noiseless file, RMS of 
+  // pedestal will be 0, and all events will be accepted, regardless of the cleaning level.
+  // For some reason the above lines do not work if made on a loop! (???)
+  clean[0]->SetSerialNumber(CT[0]);
+  clean[1]->SetSerialNumber(CT[1]);
+
+  MExtractTimeAndChargeSpline* sigextract = new MExtractTimeAndChargeSpline[NCTs];
+  MMcCalibrationUpdate*  mccalibupdate = new MMcCalibrationUpdate[NCTs];
+  MCalibrateData* calib = new MCalibrateData[NCTs];
+  MMcCalibrationCalc* mccalibcalc = new MMcCalibrationCalc[NCTs];
 
   // -------------------------------------------
-
-  //
-  // This is a demonstration program which calculates the image 
-  // parameters from a Magic Monte Carlo root file. Units of Size are
-  // for the moment, FADC counts.
-  //
-  //
   // Create a empty Parameter List and an empty Task List
   // The tasklist is identified in the eventloop by its name
   //
   MParList  plist;
-
   MTaskList tlist;
-
   plist.AddToList(&tlist);
 
@@ -89 +123 @@
   MBadPixelsCam badpix[NCTs];
 
-  for (Int_t ict = 0; ict < NCTs; ict++)
+  Float_t hi2lowratio = 10.0;
+
+  for (Int_t i = 0; i < NCTs; i++)
     {
       TString s = "MSrcPosCam;";
-      s += CT[ict];
-      src[ict].SetName(s);
-      src[ict].SetReadyToSave();
-      plist.AddToList(&(src[ict]));
+      s += CT[i];
+      src[i].SetName(s);
+      src[i].SetReadyToSave();
+      plist.AddToList(&(src[i]));
 
       TString b = "MBadPixelsCam;";
-      b += CT[ict];
-      badpix[ict].SetName(b);
-      badpix[ict].SetReadyToSave();
-      plist.AddToList(&(badpix[ict]));
-    }
+      b += CT[i];
+      badpix[i].SetName(b);
+      badpix[i].InitSize(npix[i]);
+      badpix[i].SetReadyToSave();
+      plist.AddToList(&(badpix[i]));
+
+      sigextract[i].SetSerialNumber(CT[i]);
+      sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]);
+      sigextract[i].SetRiseTimeHiGain(0.5);
+      sigextract[i].SetFallTimeHiGain(0.5);
+      sigextract[i].SetLoGainStretch(1.);
+
+      mccalibupdate[i].SetSerialNumber(CT[i]);
+      mccalibupdate[i].SetUserLow2HiGainFactor(hi2lowratio);
+      mccalibupdate[i].SetSignalType(MCalibrateData::kPhe);
+
+      calib[i].SetSerialNumber(CT[i]);
+      calib[i].SetCalibConvMinLimit(0.);
+      calib[i].SetCalibConvMaxLimit(100.); // Override limits for real data  
+      calib[i].SetCalibrationMode(MCalibrateData::kFfactor);
+      // Do not change CalibrationMode (just indicates where the cal. constants will be stored)
+      calib[i].SetSignalType(mccalibupdate[i].GetSignalType());
+
+      mccalibcalc[i].SetSerialNumber(CT[i]);
+      mccalibcalc[i].SetMinSize(200); 
+      // Minimum SIZE for an event to be used in the calculation of calibration constants.
+      // Units are ADC counts, and value depends on signal extractor!
+    }
+
+
   //
   // Now setup the tasks and tasklist:
@@ -110 +171 @@
   read.DisableAutoScheme();
 
-  read.AddFile(AnalysisFilename);
+  if (CalibrationFilename)
+    read.AddFile(CalibrationFilename->Data());
+
 
   MGeomApply*        apply = new MGeomApply[NCTs];
-
   MMcPedestalCopy*   pcopy = new MMcPedestalCopy[NCTs];
-
-  MExtractSignal*      sigextract = new MExtractSignal[NCTs];
-
-  MMcCalibrationUpdate*  mccalibupdate = new MMcCalibrationUpdate[NCTs];
-  MCalibrate* calib = new MCalibrate[NCTs];
-
-  MImgCleanStd**     clean = new MImgCleanStd*[NCTs];
-
   MHillasCalc*       hcalc = new MHillasCalc[NCTs];
-  MHillasSrcCalc*    scalc = new MHillasSrcCalc[NCTs];
 
   TString outfile = "star_";
   outfile += CT[0];
-  if (NCTs==2)
-    {
-      outfile += "_";
-      outfile += CT[1];
-    }
+  outfile += "_";
+  outfile += CT[1];
 
   //
@@ -144 +194 @@
   outfile = "star_";
   outfile += CT[0];
-  if (NCTs==2)
-    {
-      outfile += "_";
-      outfile += CT[1];
-    }
+  outfile += "_";
+  outfile += CT[1];
 
   outfile += "_";
@@ -159 +206 @@
     {
       apply[i]->SetSerialNumber(CT[i]);
-
       pcopy[i]->SetSerialNumber(CT[i]);
 
-      sigextract[i]->SetSerialNumber(CT[i]);
-      sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]);
-
-      mccalibupdate[i]->SetSerialNumber(CT[i]);
-      calib[i]->SetSerialNumber(CT[i]);
-
-      clean[i] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
-      clean[i]->SetSerialNumber(CT[i]);
-
       hcalc[i]->SetSerialNumber(CT[i]);
-      scalc[i]->SetSerialNumber(CT[i]);
+      hcalc[i].Disable(MHillasCalc::kCalcHillasSrc); 
+      // Source-dependent parameters not needed in the first loop (calibration)
 
       write1.SetSerialNumber(CT[i]);
       write2.SetSerialNumber(CT[i]);
 
-      write1.AddContainer(write1.AddSerialNumber("MMcEvt"),       "Events");
-      write1.AddContainer(write1.AddSerialNumber("MHillas"),      "Events");
-      write1.AddContainer(write1.AddSerialNumber("MHillasExt"),   "Events");
-      write1.AddContainer(write1.AddSerialNumber("MHillasSrc"),   "Events");
-      write1.AddContainer(write1.AddSerialNumber("MNewImagePar"), "Events");
-      write1.AddContainer(write1.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
-      write2.AddContainer(write2.AddSerialNumber("MMcEvt"),       "Events");
-      write2.AddContainer(write2.AddSerialNumber("MHillas"),      "Events");
-      write2.AddContainer(write2.AddSerialNumber("MHillasExt"),   "Events");
-      write2.AddContainer(write2.AddSerialNumber("MHillasSrc"),   "Events");
-      write2.AddContainer(write2.AddSerialNumber("MNewImagePar"), "Events");
-      write2.AddContainer(write2.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
-    }
-
-  if (NCTs==2)
-    {
-      write1.AddContainer("MStereoPar", "Events");
-      write2.AddContainer("MStereoPar", "Events");
-    }
+      write1.AddContainer("MMcEvt",       "Events");
+      write1.AddContainer("MHillas",      "Events");
+      write1.AddContainer("MHillasExt",   "Events");
+      write1.AddContainer("MHillasSrc",   "Events");
+      write1.AddContainer("MNewImagePar", "Events");
+      write1.AddContainer("MSrcPosCam",   "Events");
+      write2.AddContainer("MMcEvt",       "Events");
+      write2.AddContainer("MHillas",      "Events");
+      write2.AddContainer("MHillasExt",   "Events");
+      write2.AddContainer("MHillasSrc",   "Events");
+      write2.AddContainer("MNewImagePar", "Events");
+      write2.AddContainer("MSrcPosCam",   "Events");
+    }
+
+  MStereoPar* mstereo = new MStereoPar;
+  plist.AddToList(mstereo);
+
+  write1.AddContainer(mstereo, "Events");
+  write2.AddContainer(mstereo, "Events");
+  // We use MWriteRootFile::AddContainer(MParContainer* ,...) instead
+  // of using the container name as above, because in the former case the 
+  // serial number tag (indicating the telescope id) is added to the 
+  // container name, which is fine for containers of which there is one
+  // per telescope. However, the container MStereoPar is unique, since it
+  // is filled with information coming from both telescopes.
 
   write1.AddContainer("MRawRunHeader", "RunHeaders");
@@ -204 +248 @@
 
   tlist.AddToList(&read);
+
+  // Skip untriggered events (now camera simulation output contains by default all simulated events)
+  MContinue* trigger = new MContinue("(MMcTrig;1.fNumFirstLevel<1) && (MMcTrig;2.fNumFirstLevel<1)","Events");
+  tlist.AddToList(trigger);
 
   for (i = 0; i < NCTs; i++)
@@ -214 +262 @@
       tlist.AddToList(clean[i]);
       tlist.AddToList(&(hcalc[i]));
-      tlist.AddToList(&(scalc[i]));
-    }
-
+      tlist.AddToList(&(mccalibcalc[i]));
+    }
+
+
+  MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5");
+  MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5");
+  //
+  // ^^^^ Filters to divide output in two: test and train samples.
+  //
+  write1.SetFilter (&filter1);
+  write2.SetFilter (&filter2);
+
+  //
+  // Create and set up the eventloop
+  //
+  MProgressBar bar;
+  bar.SetWindowName("Calibrating");
+
+  MEvtLoop evtloop;
+  evtloop.SetProgressBar(&bar);
+  evtloop.SetParList(&plist);
+
+  //
+  // First loop: calibration loop. Go over MC events simulated with 
+  // no noise, to correlate SIZE with the number of phes and get the
+  // conversion factor (this is done by MMcCalibrationCalc).
+  //
+  if (CalibrationFilename)
+    {
+      if (!evtloop.Eventloop())
+        return;
+    }
+
+  tlist.PrintStatistics();
+
+  ///////////////////////////////////////////////////////////////////////
+
+
+  // Now prepare the second loop: go over the events you want to analyze.
+  // This time the MMcCalibrationUpdate tasks will apply the previously
+  // calculated calibration factors.
+
+  // First substitute the reading task:
+  MReadMarsFile read2("Events");
+  read2.AddFile(AnalysisFilename);  
+  read2.DisableAutoScheme();
+  tlist.AddToListBefore(&read2, &read);
+  tlist.RemoveFromList(&read);
+
+  // Delete cleaning tasks and create new ones with absolute cleaning method:
+    for (Int_t i= 0; i < NCTs; i++ )
+    {
+      tlist.RemoveFromList(clean[i]);
+      delete clean[i];
+    }
+
+  // New cleaning tasks:
+  clean[0] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
+  clean[1] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
+  clean[0]->SetMethod(MImgCleanStd::kAbsolute); 
+  clean[1]->SetMethod(MImgCleanStd::kAbsolute); 
+  clean[0]->SetSerialNumber(CT[0]);
+  clean[1]->SetSerialNumber(CT[1]);
+  tlist.AddToListBefore(clean[0],&(hcalc[0]));
+  tlist.AddToListBefore(clean[1],&(hcalc[1]));
+
+  tlist.RemoveFromList(&(mccalibcalc[0]));
+  tlist.RemoveFromList(&(mccalibcalc[1])); // Remove calibration tasks from list.
+
+  // Now calculate also source-dependent Hillas parameters:
+  for (i = 0; i < NCTs; i++)
+    hcalc[i].Enable(MHillasCalc::kCalcHillasSrc);
+
+  // Add task to calculate stereo parameters:
   MStereoCalc stereocalc;
   stereocalc.SetCTids(CT[0],CT[1]);
-
-  //
-  // FIXME: telescope coordinates must be introduced by the user, since
-  // they are not available nor in the camera file, neither in the reflector 
-  // output.
-  //
-
-  stereocalc.SetCT1coor(ctx[CT[0]-1],cty[CT[0]-1]);
-  stereocalc.SetCT2coor(ctx[CT[1]-1],cty[CT[1]-1]);
-
+  stereocalc.SetCT1coor(ctx[0],cty[0]);
+  stereocalc.SetCT2coor(ctx[1],cty[1]);
   tlist.AddToList(&stereocalc);
 
-
-  MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5");
-  MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5");
-  //
-  // ^^^^ Filters to divide output in two: test and train samples.
-  //
-
-  write1.SetFilter (&filter1);
-  write2.SetFilter (&filter2);
-
+  // Add writing tasks:
   tlist.AddToList(&filter1);
   tlist.AddToList(&write1);
@@ -246 +348 @@
   tlist.AddToList(&write2);
 
-  //
-  // Create and set up the eventloop
-  //
-  MProgressBar bar;
-
-  MEvtLoop evtloop;
-  evtloop.SetProgressBar(&bar);
-  evtloop.SetParList(&plist);
-
-  //
-  // Execute your analysis
-  //
+  bar.SetWindowName("Analyzing");
+
   if (!evtloop.Eventloop())
     return;
 
+  tlist.PrintStatistics();
+
   for (Int_t i= 0; i < NCTs; i++ )
     delete clean[i];
 
-  tlist.PrintStatistics();
+  plist.FindObject("MCalibrationChargeCam;1")->Write();
+  plist.FindObject("MCalibrationChargeCam;2")->Write();
+
+  plist.FindObject("MCalibrationQECam;1")->Write();
+  plist.FindObject("MCalibrationQECam;2")->Write();
+
+  return;
 }