Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 3545)
+++ trunk/MagicSoft/Mars/Changelog	(revision 3546)
@@ -18,4 +18,8 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2004/03/18: Raquel de los Reyes
+
+   * macros/DAQDataCheck.C
+     - macro to check the data from the DAQ system (.raw files).
 
  2004/03/18: Thomas Bretz
Index: trunk/MagicSoft/Mars/macros/DAQDataCheck.C
===================================================================
--- trunk/MagicSoft/Mars/macros/DAQDataCheck.C	(revision 3546)
+++ trunk/MagicSoft/Mars/macros/DAQDataCheck.C	(revision 3546)
@@ -0,0 +1,668 @@
+/* ======================================================================== *\
+!
+! *
+! * 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): Raquel de los Reyes, 03/2004 <mailto:reyes@gae.ucm.es>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
+!
+!
+\* ======================================================================== *//////////////////////////////////////////////////////////////////////////////
+//
+// This macro made the check of the DAQ files (.raw files).
+//
+// WARNING!: This macro only works if the directory is alphabetically ordered
+//
+///////////////////////////////////////////////////////////////////////////
+#include "MAGIC.h"
+
+void DataAnalysis(const MReadMarsFile *pedfiles,const MReadMarsFile *calfiles,const MReadMarsFile *datafiles,TString directory="./")
+{
+
+  TVirtualPad *c;
+
+  Bool_t batchmode = kTRUE;
+
+  if(!batchmode)
+    {
+      MStatusDisplay *d = new MStatusDisplay;
+      d->SetLogStream(&gLog, kTRUE);            // Disables output to stdout  
+    }
+
+  // 
+  // Set up the source run-range
+  //
+  TRegexp run("_[0-9][0-9][0-9][0-9][0-9]_");
+  TRegexp type("_[A-Z]_");
+  TString ped = "", cal = "", data = "", source = "";
+  if(pedfiles->GetEntries()!=0)
+    source = pedfiles->GetFileName();
+  else if(calfiles->GetEntries()!=0)
+    source = calfiles->GetFileName();
+  else
+    source = datafiles->GetFileName();
+
+  ped = pedfiles->GetFileName();
+  if(calfiles->GetEntries()!=0)
+    cal = calfiles->GetFileName();
+  if(datafiles->GetEntries()!=0)
+    data = datafiles->GetFileName();
+  
+  ped = ped(ped.Index(run)+1,5);
+  cal = cal(cal.Index(run)+1,5);
+  data = data(data.Index(run)+1,5);
+
+  source = source(source.Index(type)+3,source.Length());
+  source.Remove(source.Last('_'),source.Length());
+
+//    cout << ped << endl;
+//    cout << cal << endl;
+//    cout << data << endl;
+//    cout << source << endl;
+
+  // **********************************************************************
+  // ***************************** PEDESTALS ******************************
+  // **********************************************************************
+  if(pedfiles->GetEntries()==0)
+    {
+      cout << "Warning, no entries found in pedestal run(s)!!!"<< endl;
+      break;
+    }
+  MParList plist;
+  MTaskList tlist;
+  plist.AddToList(&tlist);
+  
+//    pedfiles->Print("all");
+  //
+  // Now setup the tasks and tasklist for the pedestals:
+  // ---------------------------------------------------
+  //
+  
+  //tasks
+  pedfiles->DisableAutoScheme();
+  tlist.AddToList(pedfiles);
+  
+  // pedestal tasks and containers
+  MGeomCamMagic  geomcam;
+  plist.AddToList(&geomcam);
+  MPedestalCam   pedcam;
+  plist.AddToList(&pedcam);
+  
+  MGeomApply geomapl;
+  tlist.AddToList(&geomapl);
+  MPedCalcPedRun pedcalc;
+  tlist.AddToList(&pedcalc);
+  
+  //
+  // Create and setup the eventloop
+  //
+  MEvtLoop evtloop;
+  evtloop.SetParList(&plist);
+  if(!batchmode)
+    evtloop.SetDisplay(d);
+  
+  //
+  // Execute first analysis
+  //
+  cout << "*************************" << endl;
+  cout << "** COMPUTING PEDESTALS **" << endl;
+  cout << "*************************" << endl;
+  if (!evtloop.Eventloop())
+    return;
+  
+  tlist.PrintStatistics();
+  
+  //
+  // Display the pedestal checking plots
+  //
+  if ((d = evtloop.GetDisplay()))
+    TCanvas &c1 = d.AddTab("PEDESTALS");
+  else
+    TCanvas *c1 = new TCanvas();
+  MHCamera meanped(geomcam,"MPedestalCam;ped","Pedestals");
+  meanped.SetCamContent(pedcam, 0);
+  meanped.SetCamError(pedcam, 1);
+  MHCamera rmsped(geomcam,"MPedestalCam;var","Sigma Pedestal");
+  rmsped.SetCamContent(pedcam, 2);
+  c1->Divide(1,2);
+  gPad->SetBorderMode(0);
+  c1->cd(1);
+  meanped.DrawClone("nonewEPhist");
+  c1->cd(2);
+  gPad->SetBorderMode(0);
+  c = gPad;
+  c->Divide(2,1);
+  c->cd(1);
+  meanped.DrawClone("nonewpixelindex");
+  c->cd(2);
+  rmsped.DrawClone("nonewpixelindex");
+
+  //
+  // Save canvas/display into a postcript file
+  //
+
+  if ((d = evtloop.GetDisplay()))
+    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+  else
+    c1->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps(");
+
+  if(calfiles->GetEntries()==0)
+    c1->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps]");
+
+  // **********************************************************************
+  // ***************************** CALIBRATION ****************************
+  // **********************************************************************
+  //
+  // Create a empty Parameter List and an empty Task List
+  //
+  if(calfiles->GetEntries()==0)
+    {
+      cout << "Warning, no entries found in calibration run(s)!!!"<< endl;
+      break;
+    }
+  MParList  plist2;
+  MTaskList tlist2;
+  plist2.AddToList(&tlist2);
+  
+//    calfiles->Print("all");
+  //
+  // Now setup the new tasks and tasklist for the calibration
+  // ---------------------------------------------------
+  //
+  
+  calfiles->DisableAutoScheme();
+  tlist2.AddToList(calfiles);  
+  
+  MExtractedSignalCam        sigcam;
+  MArrivalTimeCam            timecam;
+  MBadPixelsCam              badcam;
+  MCalibrationChargeCam      calcam;
+  MCalibrationChargePINDiode pindiode;
+  MCalibrationChargeBlindPix blindpix;
+  
+  MHCalibrationRelTimeCam     histtime;
+  MHCalibrationChargeCam      histcharge;
+  MHCalibrationChargePINDiode histpin;
+  MHCalibrationChargeBlindPix histblind;
+  //
+  // As long, as we don't have digital modules,
+  // we have to set the color of the pulser LED by hand
+  //
+  blindpix.SetColor(kCT1);
+  pindiode.SetColor(kCT1);
+  //
+  // Get the previously created MPedestalCam and MGeomCamMagic 
+  // into the new Parameter List
+  //
+  plist2.AddToList(&pedcam);
+  plist2.AddToList(&geomcam);
+
+  plist2.AddToList(&sigcam);
+  plist2.AddToList(&timecam);
+  plist2.AddToList(&badcam);
+  plist2.AddToList(&calcam);
+  plist2.AddToList(&histtime);
+  plist2.AddToList(&histcharge);
+  plist2.AddToList(&histpin);
+  plist2.AddToList(&histblind);
+  
+  //
+  // We saw that the signal jumps between slices,
+  // thus take the sliding window
+  //
+  MExtractSignal2        sigcalc2;
+  MExtractPINDiode       pincalc;
+  MExtractBlindPixel     blindcalc;
+  MArrivalTimeCalc2      timecalc;
+  MCalibrationChargeCalc calcalc;
+  
+  MFillH filltime( "MHCalibrationRelTimeCam"    , "MArrivalTimeCam");
+  MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); 
+  MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
+  MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");
+  
+  // Skip the MFillH draw function
+  filltime.SetBit(MFillH::kDoNotDisplay);
+  fillpin.SetBit(MFillH::kDoNotDisplay);
+  fillblind.SetBit(MFillH::kDoNotDisplay);
+  fillcam.SetBit(MFillH::kDoNotDisplay);
+  
+  //
+  // Skip the HiGain vs. LoGain calibration
+  //
+  calcalc.SkipHiLoGainCalibration();
+  
+  //
+  // Apply a filter against cosmics
+  // (was directly in MCalibrationCalc in earlier versions)
+  //
+  MFCosmics            cosmics;
+  MContinue            cont(&cosmics);
+  //
+  // Get the previously created MGeomApply into the new Task List
+  //  
+  tlist2.AddToList(&geomapl);
+
+  tlist2.AddToList(&sigcalc2);
+  tlist2.AddToList(&pincalc);
+  tlist2.AddToList(&blindcalc);
+  //
+  // In case, you want to skip the cosmics rejection,
+  // uncomment the next line
+  //
+  tlist2.AddToList(&cont);
+  //
+  // In case, you want to skip the somewhat lengthy calculation
+  // of the arrival times using a spline, uncomment the next two lines
+  //
+  tlist2.AddToList(&timecalc);
+  tlist2.AddToList(&filltime);
+
+  tlist2.AddToList(&fillpin);
+  tlist2.AddToList(&fillblind);
+  tlist2.AddToList(&fillcam);
+  tlist2.AddToList(&calcalc);
+  
+  //
+  // Create and setup the eventloop
+  //
+  MEvtLoop evtloop2;
+  evtloop2.SetParList(&plist2);
+  if(!batchmode)
+    evtloop2.SetDisplay(d);
+  
+  //
+  // Execute second analysis
+  //
+  cout << "***********************************" << endl;
+  cout << "** COMPUTING CALIBRATION FACTORS **" << endl;
+  cout << "***********************************" << endl;
+  if (!evtloop2.Eventloop())
+    return;
+  
+  tlist2.PrintStatistics();
+  
+  //
+  // print the most important results of all pixels to a file
+  //
+  //    MLog gauglog;
+  //    gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
+  //    calcam.SetLogStream(&gauglog);
+  //    calcam.Print();
+  //    calcam.SetLogStream(&gLog);
+
+  TH1 *hist;
+
+  //
+  // Display the calibration checking plots
+  //      
+  // ---------------------- Conversion factor ---------------------------
+  if ((d = evtloop2.GetDisplay()))
+    TCanvas &c2 = d.AddTab("CONV FACTOR");
+  else
+    TCanvas *c2 = new TCanvas();
+  c2->Divide(1,2);
+  c2->cd(1);
+  gPad->SetBorderMode(0);
+  c = gPad;
+  c->Divide(2,1);
+  gPad->SetBorderMode(0);
+  c->cd(1);
+  gPad->SetBorderMode(0);
+  MHCamera ConvFfactor(geomcam,"Cal;FFactorConv","Conversion Factor to photons (F-Factor method)");
+  ConvFfactor.SetCamContent(calcam, 11);
+  ConvFfactor.DrawClone("nonewpixelindex");
+  c->cd(2);
+  MHCamera ConvBlindPix(geomcam,"Cal;BlindPixConv","Conversion Factor to photons (Blind Pixel method)");
+  ConvBlindPix.SetCamContent(calcam, 17);
+  ConvBlindPix.DrawClone("nonewpixelindex");
+  c2->cd(2);
+  gPad->SetBorderMode(0);
+  MHCamera ConvPINDiode(geomcam,"Cal;PINDiodeConv","Conversion Factor to photons (PINDiode method)");
+  ConvPINDiode.SetCamContent(calcam, 23);
+  ConvPINDiode.DrawClone("nonewpixelindex");
+  // ----------------- Pixels with defects ------------------------------
+  if ((d = evtloop2.GetDisplay()))
+    TCanvas &c3 = d.AddTab("DEFECT PIXELS");
+  else
+    TCanvas *c3 = new TCanvas();
+  c3->Divide(3,2);
+  c3->cd(1);
+  MHCamera PixExc(geomcam,"Cal;Excluded","Pixels previously excluded");
+  PixExc.SetCamContent(calcam, 27);
+  PixExc.DrawClone("nonewpixelindex");
+  c3->cd(2);
+  MHCamera PixNoFit(geomcam,"Cal;NotFitted","Pixels that could not be fitted");
+  PixNoFit.SetCamContent(calcam, 28);
+  PixNoFit.DrawClone("nonewpixelindex");
+  c3->cd(3);
+  MHCamera PixNoValid(geomcam,"Cal;NotValidFit","Pixels with not valid fit results");
+  PixNoValid.SetCamContent(calcam, 29);
+  PixNoValid.DrawClone("nonewpixelindex");
+  c3->cd(4);
+  MHCamera PixFFactor(geomcam,"Cal;FFactorValid","Pixels with valid F-Factor calibration");
+  PixFFactor.SetCamContent(calcam, 35);
+  PixFFactor.DrawClone("nonewpixelindex");
+  c3->cd(5);
+  MHCamera PixBlind(geomcam,"Cal;BlindPixelValid","Pixels with valid BlindPixel calibration");
+  PixBlind.SetCamContent(calcam, 36);
+  PixBlind.DrawClone("nonewpixelindex");
+  c3->cd(6);
+  MHCamera PixPINDiode(geomcam,"Cal;PINdiodeFFactorValid","Pixels with valid PINDiode calibration");
+  PixPINDiode.SetCamContent(calcam, 37);
+  PixPINDiode.DrawClone("nonewpixelindex");
+  // ------------------------- Fitted mean charges ----------------------------
+  if ((d = evtloop2.GetDisplay()))
+    TCanvas &c4 = d.AddTab("FIT CHARGES");
+  else
+    TCanvas *c4 = new TCanvas();
+  c4->Divide(1,2);
+  MHCamera fitmeancharge(geomcam,"Cal;Charge","Fitted Mean Charges");
+  fitmeancharge.SetCamContent(calcam, 0);
+  fitmeancharge.SetCamError(calcam, 1);
+  c4->cd(1);
+  gPad->SetBorderMode(0);
+  fitmeancharge.DrawClone("nonewEPhist");
+  c4->cd(2);
+  c = gPad;
+  gPad->SetBorderMode(0);
+  c->Divide(2,1);
+  c->cd(1);
+  gPad->SetBorderMode(0);
+  fitmeancharge.DrawClone("nonewpixelindex");
+  c->cd(2);
+  gPad->SetBorderMode(0);
+  hist = (TH1*)fitmeancharge->Projection("Charge");
+  hist->SetXTitle("Charge");
+  hist->DrawCopy("hist");
+  // ----------------------- Reduced Sigma per Charge -------------------------
+  if ((d = evtloop2.GetDisplay()))
+    TCanvas &c5 = d.AddTab("RED SIGMA");
+  else
+    TCanvas *c5 = new TCanvas();
+  c5->Divide(1,2);
+  MHCamera redsigmacharge(geomcam,"Cal;RSigmaCharge","Reduced Sigma per Charge");
+  redsigmacharge.SetCamContent(calcam, 7);
+  redsigmacharge.SetCamError(calcam, 8);
+  c5->cd(1);
+  gPad->SetBorderMode(0);
+  redsigmacharge.DrawClone("nonewEPhist");
+  c5->cd(2);
+  c = gPad;
+  gPad->SetBorderMode(0);
+  c->Divide(2,1);
+  c->cd(1);
+  gPad->SetBorderMode(0);
+  redsigmacharge.DrawClone("nonewpixelindex");
+  c->cd(2);
+  gPad->SetBorderMode(0);
+  hist = redsigmacharge->Projection("RSigmaCharge");
+  hist->SetXTitle("Reduced Sigma per Charge");
+  hist->Draw("hist");
+  // ----------------------- Absolute arrival times -------------------------
+  if ((d = evtloop2.GetDisplay()))
+    TCanvas &c6 = d.AddTab("ABS TIMES");
+  else
+    TCanvas *c6 = new TCanvas();
+  c6->Divide(1,2);
+  MHCamera abstime(geomcam,"Cal;AbsTimeMean","Absolute Arrival Times");
+  abstime.SetCamContent(calcam, 42);
+  abstime.SetCamError(calcam, 43);
+  c6->cd(1);
+  gPad->SetBorderMode(0);
+  abstime.DrawClone("nonewEPhist");
+  c6->cd(2);
+  c = gPad;
+  gPad->SetBorderMode(0);
+  c->Divide(2,1);
+  c->cd(1);
+  gPad->SetBorderMode(0);
+  abstime.DrawClone("nonewpixelindex");
+  c->cd(2);
+  gPad->SetBorderMode(0);
+  hist = (TH1*)abstime->Projection("AbsTimeMean");
+  hist->SetXTitle("Absolute arrival time (time slice)");
+  hist->DrawCopy("hist");
+
+  //
+  // Save canvas/display into a postcript file
+  //
+
+  if ((d = evtloop2.GetDisplay()))
+    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+  else
+    {
+      c2->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+      c3->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+      c4->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+      c5->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+      c6->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+    }
+
+  if(calfiles->GetEntries()!=0)
+    c6->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps]");
+
+  // ********************************************************************
+  // ***************************** DATA *********************************
+  // ********************************************************************
+
+  if(datafiles->GetEntries()==0)
+    {
+      cout << "Warning, no entries found in data run(s)!!!"<< endl;
+      break;
+    }
+  //
+  // Create an empty Parameter List and an empty Task List
+  //
+  MParList  plist3;
+  MTaskList tlist3;
+  plist3.AddToList(&tlist3);
+  
+//    datafiles->Print("all");
+  //
+  // Now setup the new tasks and tasklist for the data
+  // ---------------------------------------------------
+  //
+  
+  datafiles->DisableAutoScheme();
+  tlist3.AddToList(datafiles);  
+  
+  // Task containers
+  MExtractedSignalCam sigcam;
+  //
+  // Get the previously created MPedestalCam and MGeomCamMagic
+  // into the new Parameter List
+  //
+  plist3.AddToList(&geomcam);
+  plist3.AddToList(&pedcam);
+
+  plist3.AddToList(&sigcam);
+
+  // Display containers
+  MHTimeDiffTime difftime("DiffTime","Differential time between events");
+  plist3.AddToList(&difftime);
+  MBinning binsdifftime("BinningTimeDiff");
+  binsdifftime.SetEdges(200, 0., 0.2);
+  plist3.AddToList(&binsdifftime);
+  MBinning binstime("BinningTime");
+  binstime.SetEdges(10000, 0., 1e10);
+  plist3.AddToList(&binstime);
+
+  // Data tasks
+  MExtractSignal2 signal2;
+
+  MFillH fillDiffTime(&difftime,"MRawEvtData");
+  fillDiffTime.SetBit(MFillH::kDoNotDisplay);
+  tlist3.AddToList(&fillDiffTime);
+
+  //
+  // Get the previously created MGeomApply into the new Task List
+  //
+  tlist3.AddToList(&geomapl);
+  tlist3.AddToList(&signal2);
+
+  //
+  // Create and setup the eventloop
+  //
+  MEvtLoop evtloop3;
+  evtloop3.SetParList(&plist3);
+  if(!batchmode)
+    evtloop3.SetDisplay(d);
+  
+  //
+  // Execute third analysis
+  //
+  cout << "********************" << endl;
+  cout << "** COMPUTING DATA **" << endl;
+  cout << "********************" << endl;
+  if (!evtloop3.Eventloop())
+    return;
+  
+  tlist3.PrintStatistics();
+
+  TH1 *hist1;
+  TH2 *hist2;
+
+  //
+  // Display the data checking canvas
+  //
+  // -------------------- Diffential time of events ---------------------
+  if ((d = evtloop3.GetDisplay()))
+    TCanvas &c9 = d.AddTab("TIME");
+  else
+    TCanvas *c9 = new TCanvas();
+  gStyle->SetPadGridX(kTRUE);
+  gStyle->SetPadGridY(kTRUE);
+  gStyle->SetOptFit();
+  c9->Divide(2,2);
+  hist2 = difftime.GetHist();
+  c9->cd(1);
+  gPad->SetLogy();
+  hist1 = hist2->ProjectionX("ProjX_sumtime", -1, 9999, "E");  
+  hist1->SetTitle("Distribution of \\Delta t [s]");
+  hist1->SetXTitle("\\Delta t [s]");
+  hist1->SetYTitle("Counts");
+  hist1->GetXaxis()->SetRange(0,hist1->GetXaxis()->GetNbins());
+  TF1 *f1 = new TF1("f1","[0]*exp(x*[1])",hist1->GetMaximumBin(),hist1->GetXaxis()->GetXmax());
+  hist1->Fit(f1,"");
+  hist1->DrawCopy();
+  c9->cd(2);
+  gPad->SetLogy();
+  hist1 = hist2->ProjectionX("ProjX_sumtime", -1, 9999, "E");  
+  hist1->SetTitle("Distribution of \\Delta t [s] (deadtime)");
+  hist1->SetXTitle("\\Delta t [s]");
+  hist1->SetYTitle("Counts");
+//    cout << hist1->GetXaxis()->GetBinCenter(hist1->GetMaximumBin())<< endl;
+  hist1->GetXaxis()->SetRange(0,hist1->GetMaximumBin());
+  hist1->DrawCopy();
+  c9->cd(3);
+  hist1 = difftime.GetHist()->ProjectionY("ProjY_sumtime", -1, 9999, "E");  
+  hist1->SetTitle("Distribution of time [s]");
+  hist1->SetXTitle("Time t [s]");
+  hist1->SetYTitle("Counts");
+  hist1->DrawCopy();
+  c9->cd(4);
+
+  //
+  // Save canvas/display into a postcript file
+  //
+
+  if ((d = evtloop3.GetDisplay()))
+    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
+  else
+    {
+      c9->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps)");
+    }
+  
+}
+
+void DAQDataCheck(const TString directory="/remote/data2/data/Magic-DATA/Period014/rootdata/2004_02_10/")
+{
+
+  MDirIter Next;
+  Next.AddDirectory(directory,"*_E.root",-1);
+
+//    Next.Print("all");
+
+  TString fName="file.root";
+
+  MReadMarsFile *pedfiles;
+  MReadMarsFile *calfiles;
+  MReadMarsFile *datafiles;
+
+  TString source="";
+
+  fName = Next();  // absolut path
+
+  while(!fName.IsNull())
+    {
+
+      TString file = fName(fName.Last('/')+1,fName.Length()); // root file name
+      file = file(0,file.Last('_'));
+
+//        TString run = file(file.First('_')+1,5);
+//        TString type = file(file.Index('_',file.First('_')+1)+1,1);
+      source = file(file.Last('_')+1,file.Length());
+
+      // Pedestal runs
+      if(fName.Contains("_P_"))
+	pedfiles = new MReadMarsFile("Events");
+      else
+	pedfiles->Rewind();
+
+      while(fName.Contains("_P_")&&fName.Contains(source))
+	{
+	  pedfiles->AddFile(fName);
+	  fName = Next();  
+	}
+
+      // Calibration runs
+      if(fName.Contains("_C_")||(!fName.Contains(source)))
+	calfiles = new MReadMarsFile("Events");
+      else
+	calfiles->Rewind();
+
+      while(fName.Contains("_C_")&&fName.Contains(source))
+	{
+	  calfiles->AddFile(fName);
+	  fName = Next();  
+	}
+
+      // Data runs
+      if(fName.Contains("_D_")||(!fName.Contains(source)))
+	datafiles = new MReadMarsFile("Events");
+//        else
+//  	break;
+      while(fName.Contains("_D_")&&fName.Contains(source))
+	{
+	  datafiles->AddFile(fName);
+	  fName = Next();  
+	}
+
+//        pedfiles->Print();
+//        calfiles->Print("all");
+//        datafiles->Print("all");
+
+      DataAnalysis(pedfiles,calfiles,datafiles,directory);
+
+      cout << "----------------------------------------------" << endl;
+
+    }
+
+}
+
+
+