Changeset 6196 for trunk/MagicSoft

Timestamp:

02/02/05 12:00:43 (20 years ago)

Author:

blanch

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C (modified) (14 diffs)

trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillasMC.C

@@ -50 +50 @@
   TString* OutFilename2;
 
+  TString  rffilename;
+  TString  dispfilename;
+
   // ------------- USER CHANGE -----------------
   // Comment line starting "CalibrationFileName" to disable calibration. 
@@ -57 +60 @@
 
   // File to be used in the calibration (must be a camera file without added noise)
-  CalibrationFilename = new TString("/discpepe/root_0.73mirror/wuerzburg/gammas_nonoise/Gamma_zbin0_0_7_1000to1009_w0.root");
+  CalibrationFilename = new TString("/mnt/local_wdjrico/jrico/mc/Gammas/Gamma_zbin3_90_7_1290to1299_w0_nonoise.root");
   // File to be analyzed
-  AnalysisFilename = new TString("/discpepe/root_0.73mirror/wuerzburg/gammas/Gamma_zbin*.root");
+  AnalysisFilename = new TString("/mnt/users/blanch/magic/TestSample/file*.root");
 
   // Change output file names as desired. 
   // If you want only one output (not dividing the events in Train and Test samples),
   // comment the initialization of OutFilename2.
-  OutFilename1 = new TString("/mnt/users/domingo/MAGIC/Disp05/All_gammas_0.73mirror_zbin0to12_cleaned_4035_Spline.root");
-  //  OutFilename1 = new TString("star_train.root");   // Output file name 1 (test)
+  OutFilename1 = new TString("/mnt/users/blanch/Temp/Prova.root");
   //  OutFilename2 = new TString("star_test.root");    // Output file name 2 (train)
+
+  // File to read RandomForest
+  rffilename="/mnt/users/blanch/magic/Mars-All/Mars_Standard06/mtemp/mifae/programs/RFstd.root";
+  // File to read Disp
+  dispfilename="/mnt/users/blanch/magic/Mars-All/Mars_Standard06/mtemp/mifae/programs/DISPstd.root";
 
   // Fraction of events (taken at random) which one wants to process from the 
@@ -115 +122 @@
   plist.AddToList(&badpix);
 
+  MCerPhotEvt     nphot;
+  // Stores number of phe
+  plist.AddToList(&nphot);
   
   // Now setup the tasks and tasklist:
@@ -150 +160 @@
   MMcCalibrationCalc mccalibcalc; 
   // Calculates calibration constants to convert from ADC counts to photons.
+
+  MAddGainFluctuation gainfluc;
+  //gainfluc.FillHist(1,0.5);
+  gainfluc.FillHist(0,0.1); // defaul line to not add fluctuations
+  // Adds Gain fluctuations
 
   MImgCleanStd  clean(CleanLev[0], CleanLev[1]); 
@@ -189 +204 @@
   write1.AddContainer("MMcFadcHeader",       "RunHeaders");
   write1.AddContainer("MMcTrigHeader",       "RunHeaders");
+  write1.AddContainer("MGainFluctuationCam", "RunHeaders");
 
   write1.AddContainer("MRawEvtHeader",   "Events");
@@ -201 +217 @@
   write1.AddContainer("MTopology",       "Events");
   write1.AddContainer("MPointingPos",    "Events");
+  write1.AddContainer("MHadronness",     "Events");
+  write1.AddContainer("MImageParDisp"  , "Events");
   //write1.AddContainer("MArrivalTimeCam", "Events");
   //write1.AddContainer("MCerPhotEvt",     "Events");
@@ -215 +233 @@
       write2.AddContainer("MMcFadcHeader",       "RunHeaders");
       write2.AddContainer("MMcTrigHeader",       "RunHeaders");
+      write2.AddContainer("MGainFluctuationCam", "RunHeaders");
 
       write2.AddContainer("MRawEvtHeader",   "Events");
@@ -227 +246 @@
       write2.AddContainer("MTopology",       "Events");
       write2.AddContainer("MPointingPos",    "Events");
+      write2.AddContainer("MHadronness",     "Events");
+      write2.AddContainer("MImageParDisp"  , "Events");
       //write2.AddContainer("MArrivalTimeCam", "Events");
       //write2.AddContainer("MCerPhotEvt",     "Events");
@@ -260 +281 @@
     }
 
-
-  //
-  // SECOND LOOP: Analysis loop (process file with noise)
+  //
+  // SECOND LOOP: TO READ THE RANDOM FOREST FILE(S) 
+  //
+
+  MParList   plistrf;
+  MTaskList  tlistrf;
+  MRanForest ranforest;
+  MRanTree   rantree;
+
+  if(rffilename.Length())
+    {
+      plistrf.AddToList(&tlistrf);
+      plistrf.AddToList(&ranforest);
+      plistrf.AddToList(&rantree);
+
+      MReadTree readrf("Tree",rffilename);
+      readrf.DisableAutoScheme();
+      
+      MRanForestFill rffill;
+      rffill.SetNumTrees(100);
+      
+      tlistrf.AddToList(&readrf);
+      tlistrf.AddToList(&rffill);
+      
+      MEvtLoop evtlooprf;
+      evtlooprf.SetParList(&plistrf);
+      if (!evtlooprf.Eventloop())
+        return;
+      
+      tlistrf.PrintStatistics();
+    }
+  
+
+  //
+  // THIRD LOOP: Analysis loop (process file with noise)
   //
   bar.SetWindowName("Analyzing...");
@@ -272 +325 @@
   plist.AddToList(&timecam);
   plist.AddToList(&topology);
-
+  plist.AddToList(&rantree);
+  plist.AddToList(&ranforest);
+ 
   MArrivalTimeCalc2  timecalc;
   // Calculates the arrival time as the mean time of the fWindowSize 
@@ -290 +345 @@
 
   // Change the read task by another one which reads the file we want to analyze:
-  MReadMarsFile read2("Events");
-  read2.AddFile(AnalysisFilename->Data());
+  MReadMarsFile read2("Events","/mnt/users/blanch/magic/TestSample/file*.root");
+  //  MReadMarsFile read2("Events","/mnt/users/blanch/magic/TestSample/file53.root");
+  //  read2.AddFile(AnalysisFilename->Data());
   read2.DisableAutoScheme();
   tlist.AddToListBefore(&read2, &read);
   tlist.RemoveFromList(&read);
+
 
   // Add new tasks (Islands and Topology calculations)
@@ -301 +358 @@
   //  tlist.AddToListBefore(&islclean,&hcalc);
   tlist.AddToListBefore(&topcalc,&hcalc);
+
+  // Add Task to Add Gain Fluctuations
+  tlist.AddToListBefore(&gainfluc, &hcalc);
+
+  MGeomCamMagic geomcam;// = new MGeomCam();
+  plist.AddToList(&geomcam);
+
+
+  MHillasDisplay*  disphillas=NULL;
+  disphillas = new MHillasDisplay(&nphot,&geomcam);
 
   // Analyze only the desired fraction of events, skip the rest:
@@ -312 +379 @@
   tlist.RemoveFromList(&mccalibcalc); 
 
+  // disp
+  // (read in optimum Disp parameter values)  
+  MDispParameters dispin;
+  TArrayD dispPar;
+  if(dispfilename.Length())
+    {
+      TFile inparam(dispfilename);
+      dispin.Read("MDispParameters");
+      cout << "Optimum parameter values taken for calculating Disp : " << endl;
+
+      dispPar =  dispin.GetParameters();
+      for (Int_t i=0; i<dispPar.GetSize(); i++)
+        cout << dispPar[i] << ",  ";
+      cout << endl;
+      
+      inparam.Close();
+    }
+  // Disp results container
+  MImageParDisp imagepardisp;
+  // Disp histograms
+  MHDisp hdisp;
+  hdisp.SetName("MHDispAsym");
+  hdisp.SetSelectedPos(4);
+  MFillH filldisp("MHDispAsym[MHDisp]", "");
+  
+  // Add RandomForest and Disp task computation
+  MRanForestCalc hadrcalc;
+  MDispCalc         dispcalc;  
+
+  tlist.AddToList(&dispcalc);
+  tlist.AddToList(&filldisp);
+  if(rffilename.Length())
+    tlist.AddToList(&hadrcalc);
+ 
   // Add tasks to write output:
   if (OutFilename2)
@@ -321 +422 @@
   tlist.AddToList(&write1); 
 
+  // Add display image with hillas
+  disphillas->SetPSFile();
+  disphillas->SetPSFileName("provaGF.ps");
+  disphillas->SetPause(kFALSE); 
+  //tlist.AddToList(disphillas);
+  
 
   if (!evtloop.Eventloop())