Index: /trunk/MagicSoft/Mars/mtemp/mifae/Changelog
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mifae/Changelog	(revision 4821)
+++ /trunk/MagicSoft/Mars/mtemp/mifae/Changelog	(revision 4822)
@@ -18,4 +18,8 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+  2004/08/31 Javier Rico
+    * macros/mergeClean.C, mmcCleaning, cleanComp.C
+     - Added macros for the cleaning comparison analysis
 
   2004/08/27 Javier Rico
Index: /trunk/MagicSoft/Mars/mtemp/mifae/macros/cleanComp.C
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mifae/macros/cleanComp.C	(revision 4822)
+++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/cleanComp.C	(revision 4822)
@@ -0,0 +1,445 @@
+/* ======================================================================== *\!
+!   Author(s): Javier Rico
+\* ======================================================================== */
+
+void cleanComp()
+{    
+  // general settings
+  gROOT->Reset();
+  gStyle->SetCanvasColor(0);
+  gStyle->SetCanvasBorderMode(0);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetFrameBorderMode(0);
+  gStyle->SetOptTitle(1);
+  //  gStyle->SetTitleOffset(1.7,"y");
+  //  gStyle->SetPadLeftMargin(0.15);
+  gStyle->SetOptStat(111110);
+  gStyle->SetStatColor(0);
+  gStyle->SetStatBorderSize(1);
+  gStyle->SetStatW(0.2);
+  gStyle->SetStatH(0.1);
+  gStyle->SetStatX(0.9);
+  gStyle->SetStatY(0.9);
+  
+  // define files, chain, trees, canvas...
+  TChain* ton;
+  ton= new TChain("Events");
+  
+  // data file(s)
+  ton->Add("/local_disk/jrico/mc/cleanComp0-3.root");
+  
+  // define aliases
+  ton->SetAlias("run","MRawRunHeader.fRunNumber");
+  ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
+  ton->SetAlias("energy","MMcEvt.fEnergy");
+
+  // alpha0 is for the file without noise
+  // alphaNOCL is for the file with noise and no cleaning
+  ton->SetAlias("alpha0","MHillasSrc0.fAlpha");
+  ton->SetAlias("alphaNOCL","MHillasSrcnoclean.fAlpha");
+  ton->SetAlias("alpha3015","MHillasSrc3015.fAlpha");
+  ton->SetAlias("alpha4025","MHillasSrc4025.fAlpha");
+  ton->SetAlias("alpha5035","MHillasSrc5035.fAlpha");
+  ton->SetAlias("alpha6045","MHillasSrc6045.fAlpha");
+  
+  //  ton->Scan("run:event","alpha0<-30 && energy>2500 && energy<5000");
+  
+  //  plotAlphaWidth(ton);
+   plotAlphaWidthVsEnergy(ton);
+  // plotEfficiencies(ton);
+}
+
+void plotAlphaWidth(TChain* ton)
+{
+  // define and fill histos
+  const Int_t nbins=90;
+  TH1F* halpha3015 = new TH1F("Alpha3015" ,"Alpha3015",nbins,-180,180);
+  TH1F* halpha4025 = new TH1F("Alpha4025" ,"Alpha4025",nbins,-180,180);
+  TH1F* halpha5035 = new TH1F("Alpha5035" ,"Alpha5035",nbins,-180,180);
+  TH1F* halpha6045 = new TH1F("Alpha6045" ,"Alpha6045",nbins,-180,180);
+  
+  halpha3015->SetXTitle("alpha-alpha0 (deg)");
+  halpha3015->SetYTitle("Entries");
+  
+  halpha3015->SetLineColor(1);
+  halpha4025->SetLineColor(2);
+  halpha5035->SetLineColor(3);
+  halpha6045->SetLineColor(4);
+  
+  ton->Draw("alpha3015-alpha0>>Alpha3015","alpha3015!=0");
+  ton->Draw("alpha4025-alpha0>>Alpha4025","alpha4025!=0");
+  ton->Draw("alpha5035-alpha0>>Alpha5035","alpha5035!=0");
+  ton->Draw("alpha6045-alpha0>>Alpha6045","alpha6045!=0");
+  
+  // plot
+  TCanvas *c1 = new TCanvas("c1","c1",800,500);
+  c1->cd(1);
+  
+  halpha3015->Draw("HIST");
+  halpha4025->Draw("HISTSAMES");
+  halpha5035->Draw("HISTSAMES");
+  halpha6045->Draw("HISTSAMES");
+  
+  gPad->Update();
+  TPaveStats* pvstat4025  = (TPaveStats*) halpha4025->GetListOfFunctions()->FindObject("stats");
+  TPaveStats* pvstat5035  = (TPaveStats*) halpha5035->GetListOfFunctions()->FindObject("stats");
+  TPaveStats* pvstat6045  = (TPaveStats*) halpha6045->GetListOfFunctions()->FindObject("stats");
+  
+  Float_t shifty = pvstat4025->GetY2NDC()-pvstat4025->GetY1NDC();
+  
+  pvstat4025->SetY1NDC(pvstat4025->GetY1NDC()-shifty);
+  pvstat4025->SetY2NDC(pvstat4025->GetY2NDC()-shifty);
+  pvstat5035->SetY1NDC(pvstat5035->GetY1NDC()-2*shifty);
+  pvstat5035->SetY2NDC(pvstat5035->GetY2NDC()-2*shifty);
+  pvstat6045->SetY1NDC(pvstat6045->GetY1NDC()-3*shifty);
+  pvstat6045->SetY2NDC(pvstat6045->GetY2NDC()-3*shifty);
+  
+  gPad->Modified();
+  gPad->Update();
+  
+  c1->Print("cleanComp.ps");
+}
+
+/*******************************************************************************/
+void plotAlphaWidthVsEnergy(TChain* ton)
+{
+  const Int_t nbins = 8;
+  Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
+  
+  // define the arrays that will feed the graphs
+  // value of rms of the alpha-alpha0 distribution, and errors
+  Float_t alpha0[nbins];
+  Float_t alpha3015[nbins];
+  Float_t alpha4025[nbins];
+  Float_t alpha5035[nbins];
+  Float_t alpha6045[nbins];
+
+  Float_t ealpha0[nbins];
+  Float_t ealpha3015[nbins];
+  Float_t ealpha4025[nbins];
+  Float_t ealpha5035[nbins];
+  Float_t ealpha6045[nbins];
+
+  // x axis and error
+  Float_t x[nbins];
+  Float_t ex[nbins];
+
+  TH1F* hp[nbins];
+
+  // define the function that will be used to fit the histograms
+  TF1 f2("f2","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
+  f2.SetParameters(0,100,0,02);
+  f2.FixParameter(0,0);
+  f2.SetLineWidth(1);
+
+  // fill the arrays
+  for(Int_t i=0;i<nbins;i++)
+    {
+      //      cout << "Bin " << i << endl;
+      TH1F dummy0("dummy0" ,"dummy0",100,-90,90);
+      TH1F dummy3015("dummy3015" ,"dummy3015",100,-180,180);
+      TH1F dummy4025("dummy4025" ,"dummy4025",100,-180,180);
+      TH1F dummy5035("dummy5035" ,"dummy5035",100,-180,180);
+      TH1F dummy6045("dummy6045" ,"dummy6045",100,-180,180);
+
+      TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);
+
+      // define the cuts for each of the cleaning levels
+      char cut[100];
+      sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
+      char cut3015[200];
+      sprintf(cut3015,"%s && alpha3015!=0",cut);
+      char cut4025[200];
+      sprintf(cut4025,"%s && alpha4025!=0",cut);
+      char cut5035[200];
+      sprintf(cut5035,"%s && alpha5035!=0",cut);
+      char cut6045[200];
+      sprintf(cut6045,"%s && alpha6045!=0",cut);
+
+      TF1 f1("f1","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
+      f1.SetParameters(0,100,0,02);
+      f1.FixParameter(0,0);
+
+
+      // fill the histograms to extract later the information from rms, entries, means, etc
+      ton->Draw("alpha0>>dummy0",cut);
+      ton->Draw("alpha3015-alpha0>>dummy3015",cut3015);
+      ton->Draw("alpha4025-alpha0>>dummy4025",cut4025);
+      ton->Draw("alpha5035-alpha0>>dummy5035",cut5035);
+      ton->Draw("alpha6045-alpha0>>dummy6045",cut6045);
+
+      ton->Draw("energy>>xdummy",cut);
+
+      // fit to a gaussian and get the sigma to be plotted later
+      dummy0.Fit("f1");
+      alpha0[i]   =f1.GetParameter(3);
+      dummy3015.Fit("f1");
+      alpha3015[i]=f1.GetParameter(3);
+      dummy4025.Fit("f1");
+      alpha4025[i]=f1.GetParameter(3);
+      dummy5035.Fit("f1");
+      alpha5035[i]=f1.GetParameter(3);
+      dummy6045.Fit("f1");
+      alpha6045[i]=f1.GetParameter(3);
+
+      Float_t nevt0    = dummy0.GetEntries();
+      Float_t nevt3015 = dummy3015.GetEntries();
+      Float_t nevt4025 = dummy4025.GetEntries();
+      Float_t nevt5035 = dummy5035.GetEntries();
+      Float_t nevt6045 = dummy6045.GetEntries();
+
+      // since gaussian is not so good aproximation, compute error like follows
+      ealpha0[i]   = alpha0[i]/TMath::Sqrt(2*nevt0);
+      ealpha3015[i]= alpha3015[i]/TMath::Sqrt(2*nevt3015);
+      ealpha4025[i]= alpha4025[i]/TMath::Sqrt(2*nevt4025);
+      ealpha5035[i]= alpha5035[i]/TMath::Sqrt(2*nevt5035);
+      ealpha6045[i]= alpha6045[i]/TMath::Sqrt(2*nevt6045);
+
+      // copy some histos to be plotted later
+      Char_t nom[10];
+      sprintf(nom,"hp[%02d]",i);
+      hp[i] = new TH1F(nom,nom,100,-180,180);
+      hp[i]->GetXaxis()->SetTitle("alpha (deg)");
+      hp[i]->Add(&dummy6045);
+      hp[i]->Fit("f2");
+
+      //      cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
+      x[i]  = xdummy.GetMean();            
+      //      Float_t rms = xdummy.GetRMS();
+      ex[i] = 0;//rms*rms/xdummy.GetEntries();
+    }
+
+  //#if 0
+  TCanvas* cdum = new TCanvas("cdum","cdum",600,900);
+  cdum->Divide(2,nbins/2+1);
+  for(Int_t i=0;i<nbins;i++)
+    {
+      cdum->cd(i+1);
+      hp[i]->Draw();
+      gPad->Modified();
+      gPad->Update();
+    }
+  //#endif
+  
+  // create the graphs
+  TGraphErrors* gr0    = new TGraphErrors(nbins,x,alpha0,ex,ealpha0);
+  TGraphErrors* gr3015 = new TGraphErrors(nbins,x,alpha3015,ex,ealpha3015);
+  TGraphErrors* gr4025 = new TGraphErrors(nbins,x,alpha4025,ex,ealpha4025);
+  TGraphErrors* gr5035 = new TGraphErrors(nbins,x,alpha5035,ex,ealpha5035);
+  TGraphErrors* gr6045 = new TGraphErrors(nbins,x,alpha6045,ex,ealpha6045);  
+
+  // marker/line colors
+  gr0->SetMarkerStyle(20);
+  gr3015->SetMarkerStyle(24);
+  gr4025->SetMarkerStyle(21);
+  gr5035->SetMarkerStyle(25);
+  gr6045->SetMarkerStyle(22);
+
+  gr0->SetMarkerColor(1);
+  gr3015->SetMarkerColor(2);
+  gr4025->SetMarkerColor(3);
+  gr5035->SetMarkerColor(4);
+  gr6045->SetMarkerColor(6);
+
+  gr0->SetLineColor(1);
+  gr3015->SetLineColor(2);
+  gr4025->SetLineColor(3);
+  gr5035->SetLineColor(4);
+  gr6045->SetLineColor(6);
+
+  //////////////////////
+  // plot alpha width //
+  //////////////////////
+  TCanvas *c1 = new TCanvas("c1","c1",800,500);
+  c1->cd(1);
+  gPad->SetGridx();  
+  gPad->SetGridy();  
+  gPad->SetLogx();
+
+  TGraphErrors* grfirst = gr4025;
+  grfirst->Draw("APL");
+  gr0->Draw("PL");
+  gr3015->Draw("PL");
+  gr4025->Draw("PL");
+  gr5035->Draw("PL");
+  gr6045->Draw("PL");
+
+  grfirst->GetHistogram()->SetTitle("Cleaning comparison (alpha)");
+  grfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
+  grfirst->GetHistogram()->GetYaxis()->SetTitle("#sigma (alpha-alpha_{0}) [deg]");
+  grfirst->GetHistogram()->SetMaximum(80);
+
+  TPaveLabel* label = new TPaveLabel(0.25,0.75,0.6,0.85,"MC Gammas - zbins 0-3","NDC");
+  label->SetBorderSize(0);
+  label->SetFillColor(0);
+  label->Draw();
+
+  // legend
+  TLegend* leg = new TLegend(.6,.5,.8,.8);
+  leg->SetHeader("");
+  leg->SetFillColor(0);
+  leg->SetLineColor(0);
+  leg->SetBorderSize(0);
+  leg->AddEntry(gr0,"No noise","p");
+  leg->AddEntry(gr3015,"(3.0,1.5)","p");
+  leg->AddEntry(gr4025,"(4.0,2.5)","p");
+  leg->AddEntry(gr5035,"(5.0,3.5)","p");
+  leg->AddEntry(gr6045,"(6.0,4.5)","p");
+  leg->Draw();
+
+  gPad->Modified();
+  gPad->Update();
+  // print to file
+  c1->Print("cleanAlphaWidthVSEnergy.ps");
+}
+
+/*******************************************************************************/
+void plotEfficiencies(TChain* ton)
+{
+  const Int_t nbins = 8;
+  Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
+  
+  // define the arrays that will feed the graphs
+  // gamma detection efficiency, and errors  
+  Float_t effNOCL[nbins];
+  Float_t eff3015[nbins];
+  Float_t eff4025[nbins];
+  Float_t eff5035[nbins];
+  Float_t eff6045[nbins];
+
+  Float_t eeffNOCL[nbins];
+  Float_t eeff3015[nbins];
+  Float_t eeff4025[nbins];
+  Float_t eeff5035[nbins];
+  Float_t eeff6045[nbins];
+
+  // x axis and error
+  Float_t x[nbins];
+  Float_t ex[nbins];
+
+  // fill the arrays
+  for(Int_t i=0;i<nbins;i++)
+    {
+      //      cout << "Bin " << i << endl;
+      TH1F dummyNOCL("dummyNOCL" ,"dummyNOCL",100,-90,90);
+      TH1F dummy3015("dummy3015" ,"dummy3015",100,-90,90);
+      TH1F dummy4025("dummy4025" ,"dummy4025",100,-90,90);
+      TH1F dummy5035("dummy5035" ,"dummy5035",100,-90,90);
+      TH1F dummy6045("dummy6045" ,"dummy6045",100,-90,90);
+
+      TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);
+
+      // define the cuts for each of the cleaning levels
+      char cut[100];
+      sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
+      char cutNOCL[200];
+      sprintf(cutNOCL,"%s && alphaNOCL!=0",cut);
+      char cut3015[200];
+      sprintf(cut3015,"%s && alpha3015!=0",cut);
+      char cut4025[200];
+      sprintf(cut4025,"%s && alpha4025!=0",cut);
+      char cut5035[200];
+      sprintf(cut5035,"%s && alpha5035!=0",cut);
+      char cut6045[200];
+      sprintf(cut6045,"%s && alpha6045!=0",cut);
+
+
+      // fill the histograms to extract later the information from rms, entries, means, etc
+      ton->Draw("alpha0>>dummyNOCL",cutNOCL);
+      ton->Draw("alpha0>>dummy3015",cut3015);
+      ton->Draw("alpha0>>dummy4025",cut4025);
+      ton->Draw("alpha0>>dummy5035",cut5035);
+      ton->Draw("alpha0>>dummy6045",cut6045);
+
+      ton->Draw("energy>>xdummy",cut);
+
+      Float_t nevt0    = dummyNOCL.GetEntries();
+      Float_t nevt3015 = dummy3015.GetEntries();
+      Float_t nevt4025 = dummy4025.GetEntries();
+      Float_t nevt5035 = dummy5035.GetEntries();
+      Float_t nevt6045 = dummy6045.GetEntries();
+
+      // fill the array of efficiencies
+      eff3015[i] = nevt3015/nevt0*100;
+      eff4025[i] = nevt4025/nevt0*100;
+      eff5035[i] = nevt5035/nevt0*100;
+      eff6045[i] = nevt6045/nevt0*100;
+
+      eeff3015[i] = TMath::Sqrt(nevt3015)/nevt0*100;
+      eeff4025[i] = TMath::Sqrt(nevt4025)/nevt0*100;
+      eeff5035[i] = TMath::Sqrt(nevt5035)/nevt0*100;
+      eeff6045[i] = TMath::Sqrt(nevt6045)/nevt0*100;
+
+      //      cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
+      x[i]  = xdummy.GetMean();            
+      //      Float_t rms = xdummy.GetRMS();
+      ex[i] = 0;//rms*rms/xdummy.GetEntries();
+    }
+
+  // create the graphs
+  TGraphErrors* efgr3015 = new TGraphErrors(nbins,x,eff3015,ex,eeff3015);
+  TGraphErrors* efgr4025 = new TGraphErrors(nbins,x,eff4025,ex,eeff4025);
+  TGraphErrors* efgr5035 = new TGraphErrors(nbins,x,eff5035,ex,eeff5035);
+  TGraphErrors* efgr6045 = new TGraphErrors(nbins,x,eff6045,ex,eeff6045);  
+
+  // marker/line colors
+  efgr3015->SetMarkerStyle(24);
+  efgr4025->SetMarkerStyle(21);
+  efgr5035->SetMarkerStyle(25);
+  efgr6045->SetMarkerStyle(22);
+
+  efgr3015->SetMarkerColor(2);
+  efgr4025->SetMarkerColor(3);
+  efgr5035->SetMarkerColor(4);
+  efgr6045->SetMarkerColor(6);
+
+  efgr3015->SetLineColor(2);
+  efgr4025->SetLineColor(3);
+  efgr5035->SetLineColor(4);
+  efgr6045->SetLineColor(6);
+
+  ///////////////////////
+  // plot efficiencies //
+  ///////////////////////
+  TCanvas *c1 = new TCanvas("c1","c1",800,500);
+  c1->cd(1);
+  gPad->SetGridx();  
+  gPad->SetGridy();  
+  gPad->SetLogx();
+
+  TGraphErrors* efgrfirst = efgr6045;
+  efgrfirst->Draw("APL");
+  efgr3015->Draw("PL");
+  efgr4025->Draw("PL");
+  efgr5035->Draw("PL");
+  efgr6045->Draw("PL");
+
+  efgrfirst->GetHistogram()->SetTitle("Cleaning comparison (efficiency)");
+  efgrfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
+  efgrfirst->GetHistogram()->GetYaxis()->SetTitle("Detection efficiency (%)");
+
+  // label
+  TPaveLabel* eflabel = new TPaveLabel(0.5,0.49,0.85,0.59,"MC Gammas - zbins 0-3","NDC");
+  eflabel->SetBorderSize(0);
+  eflabel->SetFillColor(0);
+  eflabel->Draw();
+ 
+  // legend
+  TLegend* efleg = new TLegend(.6,.2,.8,.5);
+  efleg->SetHeader("");
+  efleg->SetFillColor(0);
+  efleg->SetLineColor(0);
+  efleg->SetBorderSize(0);
+  efleg->AddEntry(efgr3015,"(3.0,1.5)","p");
+  efleg->AddEntry(efgr4025,"(4.0,2.5)","p");
+  efleg->AddEntry(efgr5035,"(5.0,3.5)","p");
+  efleg->AddEntry(efgr6045,"(6.0,4.5)","p");
+  efleg->Draw();
+
+  gPad->Modified();
+  gPad->Update();
+
+  // print to file
+  c1->Print("cleanEfficiency.ps");
+}
+
Index: /trunk/MagicSoft/Mars/mtemp/mifae/macros/mergeClean.C
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mifae/macros/mergeClean.C	(revision 4822)
+++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/mergeClean.C	(revision 4822)
@@ -0,0 +1,411 @@
+/* ======================================================================== *\
+!   Take different MC input files and combine them into single one
+!   Author(s): Oscar Blanch & Javier Rico
+\* ======================================================================== */
+
+void mergeClean()
+{
+  const TString fileout = "/local_disk/jrico/mc/prueba.root";
+
+  const Int_t nfiles = 11;
+  const Int_t nlevels= 5;
+
+  const Char_t* levelname[nlevels]={"noclean","3015","4025","5035","6045"};
+
+  // input files -- no noise
+  const Char_t* NonoiseOutFilename[nfiles] = 
+    {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_nonoise.root",
+     "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_nonoise.root"
+    };
+
+  // input files -- different cleaning levels
+  const Char_t* NoiseOutFilename[nlevels][nfiles]   = 
+    {{"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_noclean.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_noclean.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_noclean.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_noclean.root"},
+     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_3015.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_3015.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_3015.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_3015.root"},
+     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_4025.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_4025.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_4025.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_4025.root"},   
+     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_5035.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_5035.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_5035.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_5035.root"},
+     {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_6045.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_6045.root",   
+      "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_6045.root",
+      "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_6045.root"}
+    };   
+
+
+  // Getting input braches 
+  TChain* tno  = new TChain("Events");
+  TChain* tyes[nlevels];
+  for(Int_t i=0;i<nlevels;i++)
+    tyes[i] = new TChain("Events");
+
+  for(Int_t i=0;i<nfiles;i++)
+    {
+      tno->Add(NonoiseOutFilename[i]);
+      for(Int_t j=0;j<nlevels;j++)
+	tyes[j]->Add(NoiseOutFilename[j][i]);
+    }
+
+  // input braches for no noise file
+  MCerPhotEvt*   nphotNo;
+  MHillas*       hillasNo;
+  MHillasSrc*    hillassrcNo;
+  MMcEvt*        mcevtNo;
+  MRawEvtHeader* evtheaderNo;
+  MRawRunHeader* runheaderNo;
+
+  tno->SetBranchAddress("MCerPhotEvt.",&nphotNo);
+  tno->SetBranchAddress("MHillas.",&hillasNo);
+  tno->SetBranchAddress("MHillasSrc.",&hillassrcNo);
+  tno->SetBranchAddress("MMcEvt.",&mcevtNo);
+  tno->SetBranchAddress("MRawEvtHeader.",&evtheaderNo);
+  tno->SetBranchAddress("MRawRunHeader.",&runheaderNo);
+
+  
+  // input branches for noise files
+  MCerPhotEvt*   nphotYes[nlevels];
+  MHillas*       hillasYes[nlevels];
+  MHillasSrc*    hillassrcYes[nlevels];
+  MMcEvt*        mcevtYes[nlevels];
+  MRawEvtHeader* evtheaderYes[nlevels];
+  MRawRunHeader* runheaderYes[nlevels];
+
+  for(Int_t i=0;i<nlevels;i++)
+    {      
+      tyes[i]->SetBranchAddress("MCerPhotEvt.",&nphotYes[i]);
+      tyes[i]->SetBranchAddress("MHillas.",&hillasYes[i]);
+      tyes[i]->SetBranchAddress("MHillasSrc.",&hillassrcYes[i]);
+      tyes[i]->SetBranchAddress("MMcEvt.",&mcevtYes[i]);
+      tyes[i]->SetBranchAddress("MRawEvtHeader.",&evtheaderYes[i]);
+      tyes[i]->SetBranchAddress("MRawRunHeader.",&runheaderYes[i]);
+    }
+
+  // Setting ouput branches
+  TFile* fout = new TFile(fileout,"RECREATE");
+  TTree* to   = new TTree("Events","Events");   
+
+  // common containers
+  TBranch* bmcevt        = to->Branch("MMcEvt.",       "MMcEvt",       &mcevtNo);
+  TBranch* brawevtheader = to->Branch("MRawEvtHeader.","MRawEvtHeader",&evtheaderNo);
+  TBranch* brawrunheader = to->Branch("MRawRunHeader.","MRawRunHeader",&runheaderNo);
+
+  // noise/cleaning dependent containers
+  to->Branch("MCerPhotEvt0.",  "MCerPhotEvt",  &nphotNo);
+  to->Branch("MHillas0.",      "MHillas",      &hillasNo);
+  to->Branch("MHillasSrc0.",   "MHillasSrc",   &hillassrcNo);
+
+  //  cout << "llego1" << endl;
+  for(Int_t i=0;i<nlevels;i++)
+    {
+      //      cout << "i: " << i << endl;
+      
+      char mcerphot[100];
+      char mhillas[100];
+      char mhillassrc[100];
+
+
+      sprintf(mcerphot,  "MCerPhotEvt%s.",levelname[i]);
+      sprintf(mhillas,   "MHillas%s."    ,levelname[i]);
+      sprintf(mhillassrc,"MHillasSrc%s." ,levelname[i]);
+
+      to->Branch(mcerphot,  "MCerPhotEvt",  &nphotYes[i]);
+      to->Branch(mhillas,   "MHillas",      &hillasYes[i]);
+      to->Branch(mhillassrc,"MHillasSrc",   &hillassrcYes[i]);      
+    }
+  //  cout << "llego2" << endl;
+  
+      
+  Int_t nentries1 = (Int_t)tno->GetEntries();
+  cout << "Entries: No noise file: " << nentries1 << endl;
+  Int_t nentries2[nlevels];
+  Int_t latestk[nlevels];
+  Int_t latestevt[nlevels];
+  Int_t nfound[nlevels];
+  for(Int_t i=0;i<nlevels;i++)
+    {
+      nfound[i]=0;
+      latestk[i]=0;
+      latestevt[i]=99999999;
+      nentries2[i] = (Int_t)tyes[i]->GetEntries();
+      cout << "Entries: cleaning level [" << levelname[i] << "] :" << nentries2[i] << endl;
+    }
+
+  // loop over all events in no-noise file and look for them in the other files
+  for(Int_t i=0;i<nentries1;i++)
+    //  for(Int_t i=0;i<100;i++)
+    {      
+      tno->GetEntry(i);
+      Int_t runno = runheaderNo->GetRunNumber();
+      Int_t evtno = evtheaderNo->GetDAQEvtNumber();
+
+      cout << "Checking run " << runno<< ", event " << evtno << endl;
+
+      //check if the no-noise event is present in all of the other files
+      for(Int_t j=0;j<nlevels;j++)
+	{
+	  Int_t runyes;
+	  Int_t evtyes;
+	  cout << "Level " << levelname[j] << endl;
+	  Int_t k=latestk[j];	  
+	  while(k<nentries2[j])
+	    {
+	      //	      cout << k << endl;
+	      tyes[j]->GetEntry(k++);
+	      runyes = runheaderYes[j]->GetRunNumber();
+	      evtyes = evtheaderYes[j]->GetDAQEvtNumber();
+	      if(runyes==runno && evtyes==evtno)
+		break;
+	      // next condition speeds up selection when events are grouped by runs and
+	      // in event increasing order
+	      if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
+		break;
+	      // next condition speeds up selection when events are grouped by runs
+	      // and in run increasing order
+	      //	      if(runyes>runno)
+	      //		break;
+	    }
+
+	  if(k>=nentries2[j])
+	    {
+	      k=0;
+	      while(k<latestk[j])
+		{
+		  //		  cout << k << endl;
+		  tyes[j]->GetEntry(k++);
+		  runyes = runheaderYes[j]->GetRunNumber();
+		  evtyes = evtheaderYes[j]->GetDAQEvtNumber();
+		  if(runyes==runno && evtyes==evtno)
+		    break;
+		  // next condition speeds up selection when events are grouped by runs and
+		  // in event increasing order
+		  if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
+		    break;
+		  // next condition speeds up selection when events are grouped by runs
+		  // and in run increasing order
+		  //		  if(runyes>runno)
+		  //		    break;
+		}
+	    }
+	  
+	  // the event has not been found, assign dummy containers
+	  if(runyes=!runno || evtyes!=evtno)
+	    {
+	      cout << "Not Found!" << endl;	      
+	      for(Int_t l=j;l<nlevels;l++)
+		{
+		  nphotYes[l]->Reset();
+		  hillasYes[l]->Reset();
+		  hillassrcYes[l]->Reset();		  
+		}
+	      break;
+	    }
+	  else
+	    {
+	      cout << "Found!" << endl;
+	      latestevt[j] = evtyes;
+	      latestk[j]=k;	  
+	      nfound[j]++;
+	    }
+	}
+      to->Fill();
+    }
+
+  // loop over all events in first noise file and look for them in the no-noise file in order to save also purely noisy events  
+  cout << "***************************************" << endl;
+  cout << "*** Looking for purely noisy events ***" << endl;
+  cout << "***************************************" << endl;
+  Int_t nfound2[nlevels];
+  for(Int_t i=0;i<nlevels;i++)
+    {
+      nfound2[i]=0;
+      latestk[i]=0;
+      latestevt[i]=99999999;
+    }
+  
+  // readdress the common container addresses
+  bmcevt->SetAddress(&mcevtYes[0]);
+  brawevtheader->SetAddress(&evtheaderYes[0]);
+  brawrunheader->SetAddress(&runheaderYes[0]);
+
+  for(Int_t i=0;i<nentries2[0];i++)
+  //  for(Int_t i=0;i<100;i++)
+    {      
+      tyes[0]->GetEntry(i);
+      Int_t run = runheaderYes[0]->GetRunNumber();
+      Int_t evt = evtheaderYes[0]->GetDAQEvtNumber();
+      Int_t runno;
+      Int_t evtno;
+
+      cout << "Checking run " << run << ", event " << evt << endl;
+
+      //check if the event is present in the no-noise file
+      Int_t kk=latestk[0];	  
+      while(kk<nentries1)
+	{
+	  //	  cout << kk << endl;
+	  tno->GetEntry(kk++);
+	  runno = runheaderNo->GetRunNumber();
+	  evtno = evtheaderNo->GetDAQEvtNumber();
+	  if(run==runno && evt==evtno)
+	    break;
+	  if(run==runno && latestevt[0]<evt && evt<evtno)
+	    break;
+	  //	  if(runno>run)
+	  //	    break;
+	}
+
+      if(kk>=nentries1)
+	{
+	  kk=0;
+	  while(kk<latestk[0])
+	    {
+	      //		  cout << kk << endl;
+	      tno->GetEntry(kk++);
+	      runno = runheaderNo->GetRunNumber();
+	      evtno = evtheaderNo->GetDAQEvtNumber();
+	      if(run==runno && evt==evtno)
+		break;
+	      if(run==runno && latestevt[0]<evt && evt<evtno)
+		break;
+	      //	      if(runno>run)
+	      //		break;
+	    }
+	}
+
+      // the event is already included because it is also a no-noise event, dismiss it
+      if(run==runno && evt==evtno)
+	{
+	  cout << "Found already in the no-noise sample, continuing... " << endl;
+	  latestevt[0] = evtno;
+	  latestk[0]=kk;	  
+	  continue;
+	}
+
+      nfound2[0]++;
+      nphotNo->Reset();
+      hillasNo->Reset();
+      hillassrcNo->Reset();
+
+      // look for the event in the rest of the noisy samples
+      for(Int_t j=1;j<nlevels;j++)
+	{
+	  Int_t runyes;
+	  Int_t evtyes;
+	  cout << "Level " << levelname[j] << endl;
+	  Int_t k=latestk[j];	  
+	  while(k<nentries2[j])
+	    {
+	      //	      cout << k << endl;
+	      tyes[j]->GetEntry(k++);
+	      runyes = runheaderYes[j]->GetRunNumber();
+	      evtyes = evtheaderYes[j]->GetDAQEvtNumber();
+	      if(runyes==run && evtyes==evt)
+		break;
+	      if(runyes==run && latestevt[j]<evt && evt<evtyes)
+		break;
+	      //	      if(runyes>run)
+	      //		break;
+	    }
+
+	  if(k>=nentries2[j])
+	    {
+	      k=0;
+	      while(k<latestk[j])
+		{
+		  //		  cout << k << endl;
+		  tyes[j]->GetEntry(k++);
+		  runyes = runheaderYes[j]->GetRunNumber();
+		  evtyes = evtheaderYes[j]->GetDAQEvtNumber();
+		  if(runyes==run && evtyes==evt)
+		    break;
+		  if(runyes==run && latestevt[j]<evt && evt<evtyes)
+		    break;
+		  //		  if(runyes>run)
+		  //		    break;
+		}
+	    }
+	  
+	  // the event has not been found, assign dummy containers
+	  if(runyes=!run || evtyes!=evt)
+	    {
+	      cout << "Not Found!" << endl;	      
+	      for(Int_t l=j;l<nlevels;l++)
+		{
+		  nphotYes[l]->Reset();
+		  hillasYes[l]->Reset();
+		  hillassrcYes[l]->Reset();		  
+		}
+	      break;
+	    }
+	  else
+	    {
+	      cout << "Found!" << endl;
+	      latestevt[j] = evtyes;
+	      latestk[j]=k;	  
+	      nfound2[j]++;
+	    }
+	}
+      to->Fill();    
+    }
+
+  fout->Write();
+  fout->Close();
+
+  for(Int_t i=0;i<nlevels;i++)
+    cout << "Found " << nfound[i] << "+" << nfound2[i] << " for cleaning level " << levelname[i] << endl;
+
+}
Index: /trunk/MagicSoft/Mars/mtemp/mifae/macros/mmcCleaning.C
===================================================================
--- /trunk/MagicSoft/Mars/mtemp/mifae/macros/mmcCleaning.C	(revision 4822)
+++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/mmcCleaning.C	(revision 4822)
@@ -0,0 +1,198 @@
+/* ======================================================================== *\
+   !
+   ! *
+   ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+   ! * Software. It is distributed to you in the hope that it can be a useful
+   ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+   ! * It is distributed WITHOUT ANY WARRANTY.
+   ! *
+   ! * Permission to use, copy, modify and distribute this software and its
+   ! * documentation for any purpose is hereby granted without fee,
+   ! * provided that the above copyright notice appear in all copies and
+   ! * that both that copyright notice and this permission notice appear
+   ! * in supporting documentation. It is provided "as is" without express
+   ! * or implied warranty.
+   ! *
+   !
+   !
+   !   Author(s):
+   !
+   !   Copyright: MAGIC Software Development, 2000-2004
+   !
+   !
+   \* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//
+//  mmcCleaning - study different cleaning levels with mc
+//
+//
+/////////////////////////////////////////////////////////////////////////////
+
+
+void mmcCleaning()
+{
+  //
+  // This is a demonstration program which reads in MC camera files
+  // and produces and output with calibrated events (signal in photons
+  // for all pixels, in MCerPhotEvt containers).
+  //
+
+  // ------------- user change -----------------
+  Char_t*  NonoiseFilename = "/local_disk/jrico/mc/Gamma_zbin12_0_7_1250to1259_w0_nonoise.root";  // File to be used for the calibration (must be a camera file without added noise)
+  Char_t*  NoiseFilename   = "/local_disk/jrico/mc/Gamma_zbin12_0_7_1250to1259_w0.root";  // File to be analyzed
+
+  Char_t* NonoiseOutFilename = "/local_disk/jrico/mc/Gamma_zbin12_0_7_1250to1259_w0_cleaned_nonoise.root";  // Output file name
+  Char_t* NoiseOutFilename   = "/local_disk/jrico/mc/Gamma_zbin12_0_7_1250to1259_w0_cleaned_6045.root";     // Output file name
+
+
+  MExtractSignal    sigextract;  
+  // (other extraction methods can be used)
+
+  sigextract.SetSaturationLimit(240);
+  // Defines when to switch to low gain
+
+  // Define FADC slices to be integrated in high and low gain:
+  sigextract.SetRange(1, 14, 3, 14);
+
+  // ---------------------------------------------------------------------
+  //
+  // 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);
+
+  MPedestalCam  pedcam;
+  MGeomCamMagic geomcam;
+  MCerPhotEvt   nphot;
+
+  plist.AddToList(&pedcam);
+  plist.AddToList(&geomcam);
+  plist.AddToList(&nphot);
+
+  //
+  // Now setup the tasks and tasklist:
+  // ---------------------------------
+  //
+  MReadMarsFile read("Events");
+  read.AddFile(NonoiseFilename);
+  read.DisableAutoScheme();
+
+  MGeomApply geom; 
+  // Reads in geometry from MC file and sets the right sizes for
+  // several parameter containers.
+
+  MMcPedestalCopy   pcopy; 
+  // Copies pedestal data from the MC file run fadc header to the 
+  // MPedestalCam container.
+
+  MDisplay       display(&nphot,&geomcam);
+  MHillasDisplay display2(&nphot,&geomcam);
+
+  MPointingPosCalc pointcalc;
+  // Creates MPointingPos object and fill it with the telescope orientation
+  // information taken from MMcEvt.
+
+  MMcCalibrationUpdate  mccalibupdate;
+
+  MCalibrate calib; 
+  // MCalibrate transforms signals from ADC counts into photons. In the first
+  // loop it applies a "dummy" calibration supplied by MMcCalibrationUpdate, just 
+  // to equalize inner and outer pixels. At the end of the first loop, in the
+  // PostProcess of MMcCalibrationCalc (see below) the true calibration constants
+  // are calculated.
+
+  calib.SetCalibrationMode(MCalibrate::kFfactor);
+
+  MImgCleanStd clean(6.0,4.5);
+  clean.SetRemoveSingle(kFALSE);
+  //
+  // Applies tail cuts to image. Since the calibration is performed on 
+  // noiseless camera files, the precise values of the cleaning levels 
+  // are unimportant (in any case, only pixels without any C-photon will
+  // be rejected).
+  //
+
+  MHillasCalc hcalc; // Calculates Hillas parameters not dependent on source position.
+  MHillasSrcCalc hsrccalc;
+
+  MMcCalibrationCalc mccalibcalc; 
+  // Calculates calibration constants to convert from ADC counts to photons.
+  
+  MWriteRootFile write(NonoiseOutFilename); // Writes output
+  write.AddContainer("MRawRunHeader", "Events");
+  write.AddContainer("MMcEvt",        "Events");
+  write.AddContainer("MRawEvtHeader", "Events");
+  write.AddContainer("MCerPhotEvt",   "Events");
+  write.AddContainer("MHillas",       "Events");
+  write.AddContainer("MHillasSrc",    "Events");
+
+  tlist.AddToList(&read);
+  tlist.AddToList(&geom);
+  tlist.AddToList(&pcopy);
+  tlist.AddToList(&pointcalc);
+  tlist.AddToList(&sigextract);
+  tlist.AddToList(&mccalibupdate);
+  tlist.AddToList(&calib);
+  tlist.AddToList(&clean);
+  tlist.AddToList(&hcalc);
+  tlist.AddToList(&hsrccalc);
+  //  tlist.AddToList(&display2);
+  tlist.AddToList(&mccalibcalc);
+  tlist.AddToList(&write);
+
+  //
+  // First loop: No noise
+  //
+
+  MProgressBar bar;
+  bar.SetWindowName("No noise...");
+
+  MEvtLoop evtloop;
+  evtloop.SetProgressBar(&bar);
+  evtloop.SetParList(&plist);
+
+  if (!evtloop.Eventloop())
+    return;
+  mccalibcalc.GetHistADC2PhotEl()->Write();
+  mccalibcalc.GetHistPhot2PhotEl()->Write();
+  // Writes out the histograms used for calibration.
+
+
+  //
+  // Second loop: process file with noise
+  //
+  MReadMarsFile read2("Events");
+  read2.AddFile(NoiseFilename);
+  read2.DisableAutoScheme();
+  tlist.AddToListBefore(&read2, &read, "All");
+  tlist.RemoveFromList(&read);
+
+  MWriteRootFile write2(NoiseOutFilename); // Writes output
+  write2.AddContainer("MRawRunHeader", "Events");
+  write2.AddContainer("MMcEvt",        "Events");
+  write2.AddContainer("MRawEvtHeader", "Events");
+  write2.AddContainer("MCerPhotEvt",   "Events");
+  write2.AddContainer("MHillas",       "Events");
+  write2.AddContainer("MHillasSrc",    "Events");
+
+  //  tlist.AddToListBefore(&display2,&clean, "All");
+  tlist.AddToListBefore(&write2, &write, "All");
+  tlist.RemoveFromList(&write);
+
+
+
+  bar.SetWindowName("Cleaning 3015...");
+  clean.SetRemoveSingle();
+
+  //  tlist.RemoveFromList(&clean);
+  //  tlist.RemoveFromList(&hcalc);
+  tlist.RemoveFromList(&mccalibcalc);
+
+  if (!evtloop.Eventloop())
+    return;
+
+  tlist.PrintStatistics();
+}