Index: trunk/MagicSoft/Mars/mtemp/Changelog
===================================================================
--- trunk/MagicSoft/Mars/mtemp/Changelog	(revision 4350)
+++ trunk/MagicSoft/Mars/mtemp/Changelog	(revision 4354)
@@ -18,4 +18,12 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2004/06/25: Sebastian Commichau
+
+   * meth/macros/analysis.C 
+      -Added macro to analyze data
+
+   * meth/macros/analysis_read.C 
+      -Added macro that reads the output of analysis.C	
 
  2004/06/24: Nicola Galante
Index: trunk/MagicSoft/Mars/mtemp/meth/macros/analysis.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/meth/macros/analysis.C	(revision 4354)
+++ trunk/MagicSoft/Mars/mtemp/meth/macros/analysis.C	(revision 4354)
@@ -0,0 +1,529 @@
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// analysis.C                                                              //
+//                                                                         //
+// Author(s): S.C. Commichau, 6/2004                                       //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+
+// This macro calculates source independent and dependent image parameter 
+// distributions from data and writes them into a file called hillas.root. 
+// To look at the distributions use the macro analysis_read.C
+// Up to now no derotation and mispointing correction is taken into account!
+
+const TString inpathcal = "/ihp/MAGIC/Period016/rootdata2/2004_04_22/";
+//const TString inpath = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/";
+const TString inpath = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/";
+
+
+const TString dname = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/20040423_2375*_D_OffMrk421-6_E.root";
+
+//const TString dname = "/ihp/MAGIC/Period016/rootdata2/2004_04_23/20040423_2343*_D_Mrk421_E.root";
+
+
+const Int_t  pedrun  = 23745;//23427;  // Pedestal file
+const Int_t  calrun1 = 23205;  // Calibration files
+const Int_t  calrun2 = 0;
+const Int_t  datrun1 = 23429;  // Data files
+
+const Bool_t usedisplay = kFALSE;
+
+// Range of slices that are taken into account
+const Byte_t hifirst = 1;
+const Byte_t hilast  = 14;
+const Byte_t lofirst = 3;
+const Byte_t lolast  = 14;
+
+
+void analysis()                             
+{
+
+    /*******************************************
+     * Define some parameters for the analysis *
+     *******************************************/
+        
+    // Change the following line toggle pedestal display 
+    Bool_t DrawPedestals = kFALSE;//kTRUE; 
+    
+    // Set calibration mode
+    // kDefault, kNone, kDummy, kFfactor, kPinDiode,...
+    MCalibrate::CalibrationMode_t CalMode=MCalibrate::kDefault;
+
+    // Cleaning: Core and Tailcut
+    TArrayF clvl(2);
+    clvl[0] = 3.5;
+    clvl[1] = 3.0;	
+
+    // Set bad pixels
+    const Short_t x[16] = {330,395,329,396,389,
+			   323,388,322,384,385,
+			   386,387,321,320,319,
+			   394};
+
+    // Set output filename
+    TString oname = "hillas.root";
+   
+    
+    /***************************************
+     * Define general tasks and containers *
+     ***************************************/
+
+    MCalibrationQECam         qecam;
+    MBadPixelsCam             badcam;
+    MGeomCamMagic             geomcam;
+    MGeomApply                geomapl;
+    MExtractFixedWindow       extractor; 
+    extractor.SetRange(hifirst, hilast, lofirst, lolast);
+    //MExtractFixedWindowPeakSearch extractor;
+    //extractor.SetWindows(8,8);
+    //MExtractSlidingWindow         extractor;
+
+
+    // Iterate over runs
+    MRunIter pruns;
+    MRunIter cruns;
+    MRunIter druns;
+  
+    MReadMarsFile   dread("Events", dname);
+    dread.DisableAutoScheme();
+
+    // Pedestal runs
+    pruns.AddRun(pedrun,inpath);
+    
+    // Calibration runs
+    if (calrun2==0)
+        cruns.AddRun(calrun1,inpathcal);
+    else
+        cruns.AddRuns(calrun1,calrun2,inpath);
+
+    // Data 
+/*    if (1)//datrun2==0)
+	druns.AddRun(datrun1,inpath);
+    else
+	druns.AddRuns(datrun1,datrun2,inpath);//,datrun3,datrun4,datrun5,datrun6,inpath);
+*/
+
+    // Look for the pulser colors 
+    MCalibrationCam::PulserColor_t color;
+ 
+    if (calrun1 < 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("Couldn't find the pulser colour, type 'q' to exit, "
+				    "'green', 'blue', 'uv' or 'ct1' to go on: ");
+	    timer.TurnOff();
+	    
+	    if (input=="q\n")
+		return;
+	    
+	    if (input=="green\n")
+	    {
+		color = MCalibrationCam::kGREEN;
+		break;
+	    }
+	    if (input=="blue\n")
+	    {
+		color = MCalibrationCam::kBLUE;
+		break;
+	    }
+	    if (input=="uv\n" || input=="UV\n")
+	    {
+		color = MCalibrationCam::kUV;
+		break;
+	    }
+	    if (input=="ct1\n" || input=="CT1\n")
+	    {
+		color = MCalibrationCam::kCT1;
+		break;
+	    }
+	}
+    }
+    
+ 
+    /**********************************
+     * 1st Loop: Pedestal Computation *
+     **********************************/
+    
+    // If pixels have to be excluded from the beginning, read 
+    // an ASCII-file with the corresponding pixel numbers (MBadPixelsCam)
+    badcam.AsciiRead("badpixels.txt");
+    
+    MJPedestal                pedloop;
+    pedloop.SetInput(&pruns);   
+    pedloop.SetExtractor(&extractor);
+    pedloop.SetBadPixels(badcam);
+    
+    if (!pedloop.Process())
+	return;
+        
+    // Print the content of MPedestalCam
+    MPedestalCam *pedcam = (MPedestalCam*)pedloop.GetPedestalCam();
+    pedcam->Print(); 
+    
+    // Draw Pedestals and RMS into a camera - ADC counts
+    if (DrawPedestals)
+    {
+	MHCamera *disp0 = new MHCamera(geomcam, "Pedestals", "");
+	MHCamera *disp1 = new MHCamera(geomcam, "PedestalRMS", "");
+	
+	disp0->SetMinimum(0);
+	disp1->SetMinimum(0);
+	
+	TCanvas *c1 = new TCanvas("Pedestals", "Pedestals", 50, 50, 400, 400);
+	TCanvas *c2 = new TCanvas("RMS", "RMS", 460, 50, 400, 400);
+	
+	c1->SetBorderMode(0);
+	c1->SetFillColor(0);
+	
+	c2->SetBorderMode(0);
+	c2->SetFillColor(0);
+	
+	gStyle->SetOptStat(1101);
+	gStyle->SetStatColor(0);
+	gStyle->SetStatBorderSize(1);
+	
+	gPad->SetBorderMode(0);
+	
+	disp0.SetPrettyPalette();
+	disp1.SetPrettyPalette();
+	
+	c1->cd();
+	disp0.Draw();
+	
+	gPad->cd(1);
+	
+	disp0.SetCamContent(*pedcam,0);
+	disp0.SetCamError(*pedcam,1);
+	
+	c2->cd();
+	disp1.Draw();
+	
+	disp1.SetCamContent(*pedcam,2);
+	disp1.SetCamError(*pedcam,3);
+    }
+    
+   
+    /*************************
+     * 2nd Loop: Calibration *
+     *************************/
+
+    MJCalibration     calloop;
+    //calloop.SetColor(color);
+    calloop.SetInput(&cruns);
+    calloop.SetExtractor(&extractor);
+    //calloop.SetQECam(qecam);
+    calloop.SetBadPixels(pedloop.GetBadPixels());
+    
+    // Choose the arrival time Extractor:
+    //MExtractTimeHighestIntegral timeext;
+    //MExtractTimeSpline timeext;
+    // Set Ranges or Windows
+    //timeext.SetRange(0,14,6,14);
+    //calloop.SetTimeExtractor(&timeext);
+    
+    // To calibrate times enable the following line
+    //Bool_t UseTimes = kFALSE;
+    // Bool_t calibrateTimes = kFALSE;
+    //calloop.SetRelTimeCalibration(UseTimes);
+    
+    if (!calloop.Process(pedloop.GetPedestalCam()))
+    return;
+    
+    /**********************************
+     * 3rd Loop: Pedestal Calibration *
+     **********************************/
+
+    // Create a Parameter and a Task List
+    MParList                  plist3;
+    MTaskList                 tlist3;
+    plist3.AddToList(&tlist3);
+
+    // Parameter Containers
+    MCerPhotEvt               evtphot;
+    MPedPhotCam               pedphot;
+    MExtractedSignalCam       sigcam;
+    
+    // Add also some containers from previous loops
+    plist3.AddToList(&geomcam);
+    plist3.AddToList(&pedloop.GetPedestalCam());
+
+    plist3.AddToList(&calloop.GetCalibrationCam());
+    plist3.AddToList(&calloop.GetQECam());          // I have to check this! 
+    //plist3.AddToList(&calloop.GetRelTimeCam());
+    //plist3.AddToList(&calloop.GetBadPixels());
+    plist3.AddToList(&badcam);
+    plist3.AddToList(&sigcam);
+    plist3.AddToList(&evtphot);
+    plist3.AddToList(&pedphot);
+
+    
+    // Tasks
+    MReadMarsFile read3("Events");
+    static_cast<MRead&>(read3).AddFiles(pruns);  
+    read3.DisableAutoScheme();
+    
+    //MExtractSignal            extsig;
+    //extsig.SetRange(hifirst, hilast, lofirst, lolast);
+    MCalibrate                calphot;
+    calphot.SetCalibrationMode(CalMode);
+    MPedPhotCalc              pedphotcalc;  
+    
+    tlist3.AddToList(&read3);
+    tlist3.AddToList(&geomapl);
+    tlist3.AddToList(&extractor);
+    //tlist3.AddToList(&extsig);
+    tlist3.AddToList(&calphot);
+    tlist3.AddToList(&pedphotcalc);
+  
+    
+    // Create and setup the eventloop
+    MProgressBar              bar3;
+    bar3.SetWindowName("Calibrating Pedestals...");
+    
+    MEvtLoop                  evtloop3;
+    evtloop3.SetProgressBar(&bar3);
+    evtloop3.SetParList(&plist3);
+    
+    if (!evtloop3.Eventloop())
+	return;
+    
+    tlist3.PrintStatistics();
+
+
+    /****************************
+     * 4th Loop: Photonize Data *
+     ****************************/
+
+    // Create a Parameter and a Task List
+    MParList                  plist4;
+    MTaskList                 tlist4;
+    plist4.AddToList(&tlist4);
+
+    // Parameter Containers
+    MHillas                   hillas;
+    //MDCA                      dca;
+    MHillasSrc                hillassrc;
+    MSrcPosCam                source;
+    source.SetX(0.0*189/0.6);
+    source.SetY(0.0*189/0.6);
+    MRawRunHeader             runhead;
+    MArrivalTimeCam           timecam;
+    
+    // False Source stuff
+    MTime                     time;
+    MObservatory              obsv;
+    MPointingPos              poin;
+    
+    hillassrc.SetSrcPos(&source);
+    
+
+    // Add also some containers from previous loops
+    //plist4.AddToList(&timecam);
+    plist4.AddToList(&geomcam);
+    plist4.AddToList(&pedloop.GetPedestalCam());
+    plist4.AddToList(&calloop.GetCalibrationCam());
+    //plist4.AddToList(&badcam);
+    plist4.AddToList(&calloop.GetQECam());
+    //plist4.AddToList(&calloop.GetRelTimeCam());
+    //plist4.AddToList(&calloop.GetBadPixels());
+    plist4.AddToList(&source);
+    plist4.AddToList(&hillas);
+    plist4.AddToList(&hillassrc);
+    //plist4.AddToList(&dca);    
+    plist4.AddToList(&evtphot);
+    plist4.AddToList(&pedphot);
+
+    plist4.AddToList(&time);
+    plist4.AddToList(&obsv);
+    plist4.AddToList(&poin);
+    
+    
+    //MArrivalTime times;
+    //plist4.AddToList(&times);
+    
+    // Tasks
+/*  MReadMarsFile read4("Events");
+    static_cast<MRead&>(read4).AddFiles(druns);  
+    read4.DisableAutoScheme();
+*/
+    // Set bad pixels
+    MBlindPixelCalc           blind;
+
+    const TArrayS bp(16,(Short_t*)x);
+    blind.SetPixelIndices(bp);
+
+    MImgCleanStd              clean(clvl[0],clvl[1]);
+    MHillasCalc               hilcalc;
+    MHillasSrcCalc            hilsrc;
+    //MDCACalc                  dcacalc;
+    
+    
+    
+
+    MWriteRootFile write(oname,"RECREATE"); //Alternative options: NEW, UPDATE
+
+    //write.AddContainer("MDCA"         , "Parameters");
+    write.AddContainer("MHillas"        , "Parameters");
+    write.AddContainer("MHillasSrc"     , "Parameters");
+    write.AddContainer("MHillasExt"     , "Parameters");
+    write.AddContainer("MNewImagePar"   , "Parameters");
+    write.AddContainer("MRawEvtHeader"  , "Parameters");
+    write.AddContainer("MRawRunHeader"  , "Parameters");
+    write.AddContainer("MConcentration" , "Parameters");
+    write.AddContainer("MSrcPosCam"     , "Parameters");
+
+    
+    // False Source stuff
+    MHFalseSource             fshist;
+    fshist.SetAlphaCut(12.5);
+    //fshist.SetBgMean(35);
+    MFillH                   hfill(&fshist, "MHillas");
+    //MFillH                   hfill("MHFalseSource", "MHillas");
+
+    //MFalseSrcCalc     falsecalc;
+    //falsecalc.SetAlphaCut(12.5);
+    //falsecalc.SetDistCuts(.25,1.0);
+    //falsecalc.SetSizeCut(1000);
+    //falsecalc.SetWidthCuts(.04,.14);
+    //falsecalc.SetLengthCuts(.14,.26);
+    //falsecalc.SetOutput("ON.root");//outname);
+    //falsecalc.SetRefPoint(0,0);
+
+    //tlist4.AddToList(&read4);
+    tlist4.AddToList(&dread);
+    tlist4.AddToList(&geomapl);
+    tlist4.AddToList(&extractor);
+    //tlist4.AddToList(&blind);
+    //tlist4.AddToList(&timecalc);
+    tlist4.AddToList(&calphot);
+    //tlist4.AddToList(&timecal);
+    tlist4.AddToList(&clean);
+    tlist4.AddToList(&hilcalc);
+    //tlist4.AddToList(&dcacalc);
+    tlist4.AddToList(&hilsrc);
+//    tlist4.AddToList(&falsecalc);
+    //tlist4.AddToList(&hfill);
+    tlist4.AddToList(&write);  
+
+    // Create and setup the eventloop
+    MProgressBar              bar4;
+    bar4.SetWindowName("Analyzing Data...");
+
+    MEvtLoop                  evtloop4;
+    evtloop4.SetProgressBar(&bar4);
+    evtloop4.SetParList(&plist4);
+
+    if (!evtloop4.Eventloop())
+	return;
+
+    tlist4.PrintStatistics();
+
+    //fshist.Draw();
+
+
+}
+
+
+Bool_t HandleInput()
+{
+    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
+
+    while (1)
+    {
+        // While reading the input process GUI events asynchronously
+        timer.TurnOn();
+        TString input = Getline("Type 'q' to exit, <return> to go on: ");
+        timer.TurnOff();
+
+        if (input=="q\n")
+            return kFALSE;
+
+        if (input=="\n")
+            return kTRUE;
+    };
+
+    return kFALSE;
+}
+
+
+// Find the correct pulser color, done by M. Gaug
+MCalibrationCam::PulserColor_t FindColor(MDirIter* run) 
+{
+  
+  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;      
+}
Index: trunk/MagicSoft/Mars/mtemp/meth/macros/analysis_read.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/meth/macros/analysis_read.C	(revision 4354)
+++ trunk/MagicSoft/Mars/mtemp/meth/macros/analysis_read.C	(revision 4354)
@@ -0,0 +1,568 @@
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// analysis_read.C                                                         //
+//                                                                         //
+// Author(s): S.C. Commichau, 3/2004                                       //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+
+// This macro reads the output of analysis.C and plots some of the
+// image parameters. It also produces an ALPHA plot and calculates
+// the significance
+
+void analysis_read(){
+
+    gROOT->Reset();
+    
+    // Some parameters for the ALPHA-plot
+    const Int_t inibin = 9;
+    const Int_t nbins  = 18;
+
+    // ON data chain
+    TChain* CON  = new TChain("Parameters");
+
+    // OFF data chain
+    TChain* COFF = new TChain("Parameters");
+
+    // Add the root file(s) to the chain and treat it like a tree
+    CON->Add("~/Mars/on1.root");
+    COFF->Add("~/Mars/off1.root");
+
+    //CON->Add("~/no_back/hillas_on_all23.root");
+    //COFF->Add("~/no_back/hillas_off_all23.root");
+
+    TString title = "Mrk421 - 04/23/2004";
+
+    // Set some global style options
+    //gStyle->SetOptDate(10);
+    //gStyle->SetDateX(.1);
+    //gStyle->SetDateY(.905);
+    gStyle->SetOptFit(0);
+    gStyle->SetOptStat(110011);
+    gStyle->SetFrameBorderMode(0);
+    //gStyle->SetPalette(1);
+    gStyle->SetPalette(1,0);
+    gStyle->SetFrameBorderSize(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetFrameFillColor(0);
+    gStyle->SetTitleFillColor(0);
+    gStyle->SetTitleBorderSize(0);
+    gStyle->SetStatColor(0);
+    gStyle->SetStatBorderSize(1);
+    // Put the statistics box into the right corner!
+    gStyle->SetStatX(.9);
+    gStyle->SetStatY(.9);
+    //gStyle->SetStatTextColor(2);
+
+    // Set alias for each quantity to simplify the handling
+    CON->SetAlias("length","MHillas.fLength*0.6/189.0");
+    CON->SetAlias("width","MHillas.fWidth*0.6/189.0");
+    CON->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
+    CON->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
+    CON->SetAlias("size","MHillas.fSize");
+    CON->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
+    CON->SetAlias("alpha","MHillasSrc.fAlpha");
+    CON->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");
+    //CON->SetAlias("delta","MDCA.fDelta1*180/3.1415");
+    //CON->SetAlias("dca","MDCA.fDCA*0.6/189.0");
+
+    COFF->SetAlias("length","MHillas.fLength*0.6/189.0");
+    COFF->SetAlias("width","MHillas.fWidth*0.6/189.0");
+    COFF->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
+    COFF->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
+    COFF->SetAlias("size","MHillas.fSize");
+    COFF->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
+    COFF->SetAlias("alpha","MHillasSrc.fAlpha");
+    COFF->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");
+
+
+    // Cuts...
+    // Cut for ALPHA
+    TString cut1 = "(dist>0.2) && (dist<.7) && (width>0.04) && (width<0.14) && (length>0.14) && (length<0.26) && (size>500)";
+        
+    TString cut2 = "";//size>1000 && size<2000";
+
+    TString cut3 = "";
+
+    TString cut4 = "";
+
+    TString cut5 = "";
+    
+
+    // And now the draw area
+    
+    //************ Width vs Length ************
+
+    TCanvas* CWL = new TCanvas("CWL","Canvas",0,0,300,300);
+    CWL->cd();
+    CWL->SetBorderMode(0); 
+    gPad->SetBorderMode(0);
+    gStyle->SetOptStat(0);
+            
+    TH2F* hWL = new TH2F("hWL","",200,0,0.7,200,0,1);
+    
+    // Leaf, Cuts, Draw option, use >>+ to append data to an ex. histogram
+    CON->Draw("length:width>>hWL",cut2);
+
+    hWL->SetTitle(title); 
+    hWL->GetXaxis()->SetTitle("WIDTH [#circ]");
+    hWL->GetYaxis()->SetTitle("LENGTH [#circ]");
+    hWL->GetYaxis()->SetTitleOffset(1.2);
+    hWL->SetMarkerColor(2);
+     
+    hWL->SetMaximum(150);
+    hWL->Draw("COLZ");//CONT
+
+
+    //************ Width and Length ************
+
+    // Create Canvas and set Canvas some options
+    TCanvas* Cwl = new TCanvas("Cwl","Canvas",325,0,300,300);
+    Cwl->cd();
+    Cwl->SetBorderMode(0);   
+    gPad->SetBorderMode(0);
+    gStyle->SetOptStat(110011);
+
+    TH1F* hLength = new TH1F("hLength","",100,0,1.1);
+        
+    CON->Draw("length>>hLength",cut5,"");
+
+    hLength->SetTitle(title); 
+    hLength->GetXaxis()->SetTitle("LENGTH and WIDTH [#circ]");
+    hLength->GetYaxis()->SetTitle("Entries");
+    hLength->GetYaxis()->SetTitleOffset(1.6);
+    hLength->SetFillColor(50);
+    hLength->SetLineColor(50);
+    hLength->SetLineWidth(2);
+    hLength->SetFillStyle(3005);
+    //hLength->SetMaximum(55000);
+    hLength->SetName("LENGTH");
+
+    TH1F* hWidth = new TH1F("hWidth","",100,0,1.1);
+
+    CON->Draw("width>>hWidth",cut5,"");
+ 
+    gStyle->SetStatTextColor(4);
+    hWidth->SetFillColor(4);
+    hWidth->SetLineColor(4);
+    hWidth->SetFillStyle(3004);
+    hWidth->SetLineWidth(2);
+    hWidth->SetName("WIDTH");
+
+    TLegend* leg = new TLegend(.59,.75,.8,.84);
+    leg->AddEntry(hLength,"LENGTH","F");
+    leg->AddEntry(hWidth,"WIDTH","F");
+    leg->SetFillColor(0);
+    leg->SetLineColor(1);
+    leg->SetBorderSize(0);
+
+    hLength->Draw();//DrawNormalized or DrawCopy are possible options
+    hWidth->Draw("sames");
+    //leg->Draw("same");
+
+    // Beautify the statistic boxes!
+    gPad->Update();
+    TPaveStats* pavstat = (TPaveStats*) hLength->GetListOfFunctions()->FindObject("stats");
+    if(pavstat)
+    {
+	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
+	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
+	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
+	pavstat->SetTextColor(50);
+    }
+
+    gPad->Modified();
+    gPad->Update();
+
+
+    //************ Dist ************
+    
+    TCanvas* CDist = new TCanvas("CDist","Canvas",650,0,300,300);
+    CDist->cd();
+    CDist->SetBorderMode(0); 
+    gStyle->SetOptStat(110011);
+    //CDist->SetLogy();
+
+    TH1F* hDist = new TH1F("hDist","",90,0,1.8);
+
+    CON->Draw("dist>>hDist");
+
+    hDist->GetXaxis()->SetTitle("DIST [#circ]");
+    hDist->GetYaxis()->SetTitle("Entries");
+    hDist->GetYaxis()->SetTitleOffset(1.6);
+    //hDist->GetXaxis()->SetTitleOffset(1.3);
+    hDist->SetLineColor(4);
+    hDist->SetMarkerColor(4);
+    hDist->SetFillColor(4);
+    hDist->SetFillStyle(3004);
+    hDist->SetLineWidth(2);
+    hDist->SetTitle(title);
+    hDist->SetName("DIST");
+    
+    hDist->Draw();
+
+
+    //************ MeanY vs MeanX ************
+    
+    TCanvas* CXY = new TCanvas("CXY","Canvas",975,0,300,300);
+    CXY->cd();
+    CXY->SetBorderMode(0);  
+    gStyle->SetOptStat(0);
+    
+    TH2F* hXY = new TH2F("hXY","",500,-1.5,1.5,500,-1.5,1.5);
+
+    CON->Draw("meany:meanx>>hXY",cut3,"");
+
+    TString titlexy = TString(title);
+
+    hXY->SetTitle(titlexy); 
+    hXY->GetXaxis()->SetTitle("MeanX [#circ]");
+    hXY->GetYaxis()->SetTitle("MeanY [#circ]");
+    hXY->GetYaxis()->SetTitleOffset(1.2);
+    hXY->SetMaximum(15);
+     
+    hXY->Draw("COLZ");
+
+
+    //************ DCA ************
+    /*     
+    TCanvas* CDCA = new TCanvas("CDCA","Canvas",10,400,300,300);
+    CDCA->cd();
+    CDCA->SetBorderMode(0); 
+    gStyle->SetOptStat(110011);
+    
+    TH1F* hDCA = new TH1F("hDCA","",100,-1.5,1.5);
+
+    CON->Draw("dca>>hDCA");
+
+    hDCA->GetXaxis()->SetTitle("DCA [#circ]");
+    hDCA->GetYaxis()->SetTitle("Entries");
+    hDCA->GetYaxis()->SetTitleOffset(2);
+    hDCA->SetLineColor(4);
+    hDCA->SetFillColor(4);
+    hDCA->SetFillStyle(3004);
+    hDCA->SetTitle(title);
+    hDCA->SetName("DCA");
+
+    hDCA->Draw();
+    */    
+
+
+    //************ Delta ************
+     
+    TCanvas* CDelta = new TCanvas("CDelta","Canvas",0,335,300,300);
+    CDelta->cd();
+    CDelta->SetBorderMode(0); 
+    gStyle->SetOptStat(110011);
+
+    TH1F* hDelta = new TH1F("hDelta","",100,0,90);
+
+    CON->Draw("delta>>hDelta",cut3);
+    
+    hDelta->GetXaxis()->SetTitle("#delta [#circ]");
+    hDelta->GetYaxis()->SetTitle("Entries");
+    hDelta->GetYaxis()->SetTitleOffset(1.4);
+    hDelta->SetLineColor(4);
+    hDelta->SetFillColor(4);
+    hDelta->SetFillStyle(3004);
+    hDelta->SetLineWidth(2);
+    hDelta->SetTitle(title);
+    hDelta->SetName("Delta");
+
+    hDelta->Draw();
+
+
+    //************ Size ************
+     
+    TCanvas* CSize = new TCanvas("CSize","Canvas",325,335,300,300);
+    CSize->cd();
+    CSize->SetBorderMode(0); 
+    gStyle->SetOptStat(110011);
+
+    TH1F* hSize = new TH1F("hSize","",100,0.1,5000);
+
+    CON->Draw("size>>hSize");
+    
+    hSize->GetXaxis()->SetTitle("Size [Photons]");
+    hSize->GetYaxis()->SetTitle("Entries");
+    hSize->GetYaxis()->SetTitleOffset(1.6);
+    hSize->GetXaxis()->SetTitleOffset(1.2);
+    hSize->SetLineColor(4);
+    hSize->SetFillColor(4);
+    hSize->SetLineWidth(2);
+    hSize->SetFillStyle(3004);
+    hSize->SetTitle(title);
+    hSize->SetName("SIZE");
+    
+    hSize->Draw();
+
+
+    //************ Alpha ************
+       
+    TCanvas* CAlpha = new TCanvas("CAlpha","Canvas",650,335,300,300);
+    CAlpha->cd();
+    CAlpha->SetBorderMode(0); 
+    CAlpha->SetGridx();
+    CAlpha->SetGridy();
+    gStyle->SetOptStat(11);
+    
+    TH1F *hAlpha = (TH1F*)new TH1F("hAlpha","Alpha ON",nbins,0,90);
+    TH1F *hAlphaoff = (TH1F*)new TH1F("hAlphaoff","Alpha OFF",nbins,0,90);
+
+    CON->Draw("abs(alpha)>>hAlpha",cut1);
+
+    COFF->Draw("abs(alpha)>>hAlphaoff",cut1);
+
+
+    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+    // Normalize ON and OFF data
+    cout << "******************************************************" << endl;
+    cout << " Bin error before scaling: " << hAlphaoff->GetBinError(2) << endl;
+
+    // If TH1::Sumw2 has been called before filling, the sum of squares of
+    // weights is also stored - I do not understand this!!
+    hAlphaoff->Sumw2();
+    
+    Float_t level = 0, leveloff = 0;
+    
+    // Calculate the mean of entries beside the signal region
+    for(Int_t ibin = inibin; ibin<=nbins; ibin++)
+	level += hAlpha->GetBinContent(ibin);
+    level /= (nbins-inibin+1);
+    
+    cout << " Mean beside signal: " << level << endl;
+
+    for(Int_t ibin = inibin; ibin<=nbins; ibin++)
+	leveloff += hAlphaoff->GetBinContent(ibin);
+    leveloff /= (nbins-inibin+1);
+
+    // Compare the mean values of both histograms
+    if(leveloff>0)
+	const Float_t norm = level/leveloff;//*nbins/hAlphaoff->GetEntries();
+    cout << " Normalizing by factor " << norm <<endl;
+    hAlphaoff->Scale(norm);
+    
+    cout << " Bin error after scaling: " << hAlphaoff->GetBinError(2) << endl;
+
+    cout << " Cut on (ALPHA): " << cut1 << endl;
+    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+
+   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+    // Significance calculation
+    Float_t significance = 0;
+    const Int_t signbins = 3;
+    
+    Int_t Non  = 0;
+    Int_t Noff = 0;
+
+    cout << " Significance Bins: " << signbins << endl;
+
+    // Determine number of excess events and off events in the same region 
+    for(Int_t ibin = 1; ibin<=signbins;ibin++)
+    {
+	Non +=hAlpha->GetBinContent(ibin);
+	Noff+=hAlphaoff->GetBinContent(ibin);
+    }
+
+    cout << " Non = " << Non << " Noff = " << Noff << endl;
+
+    significance = Significance(Non, Noff, 1.0, 0); //norm
+
+    cout << " Significance: " << significance << endl;
+    cout << "******************************************************" << endl;
+    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
+
+    // Set some draw options and draw the histograms
+    hAlphaoff->GetXaxis()->SetTitle("ALPHA [#circ]");
+    hAlphaoff->GetYaxis()->SetTitle("Entries");
+    hAlphaoff->GetYaxis()->SetTitleOffset(1.4);
+    hAlphaoff->GetXaxis()->SetTitleOffset(1.1);
+
+    char legentry[128];
+    sprintf(legentry,"  (Sign.: %2.3f #sigma,  Norm.: %2.3f,  |ALPHA| #leq %d#circ)", significance, norm, 90/nbins*signbins);
+
+    hAlphaoff->SetTitle(title+legentry);
+    hAlphaoff->SetName("OFF source");
+    hAlphaoff->SetLineColor(4);
+    hAlphaoff->SetLineWidth(2);
+    hAlphaoff->SetMarkerColor(8);
+    hAlphaoff->SetMarkerStyle(21);
+    hAlphaoff->SetMarkerSize(0);
+    hAlphaoff->SetFillColor(18);
+    hAlphaoff->SetMaximum(hAlpha->GetMaximum()+200);
+
+    hAlpha->SetLineColor(1);
+    hAlpha->SetLineWidth(2);
+    hAlpha->SetMarkerColor(50);
+    hAlpha->SetMarkerStyle(20);
+    hAlpha->SetMarkerSize(0.9);
+    hAlpha->SetName("ON source");
+    //hAlpha->SetMaximum(hAlpha->GetMaximum());
+
+    hAlphaoff->Draw("eH");
+    hAlpha->Draw("e1sames");
+
+    // Beautify the statistic boxes!
+    gPad->Update();
+    TPaveStats* pavstat = (TPaveStats*) hAlphaoff->GetListOfFunctions()->FindObject("stats");
+    if(pavstat)
+    {
+	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
+	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
+	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
+	pavstat->SetTextColor(4);
+    }
+
+    TPaveStats* pavstaton = (TPaveStats*) hAlpha->GetListOfFunctions()->FindObject("stats");
+    if(pavstaton)
+    {
+	pavstaton->SetTextColor(50);
+    }
+
+    gPad->Modified();
+    gPad->Update();
+
+
+    //************ Length/Size ************
+    
+    TCanvas* CLSiz = new TCanvas("CLSiz","Canvas",975,335,300,300);
+    CLSiz->cd();
+    CLSiz->SetBorderMode(0); 
+    gStyle->SetOptStat(110011);
+
+    TH1F* hLenSiz = new TH1F("hLenSiz","",100,0,0.005);
+
+    CON->Draw("length/size>>hLenSiz",cut2,"");
+
+    hLenSiz->GetXaxis()->SetTitle("LENGTH/SIZE [#circ/Photons]");
+    hLenSiz->GetYaxis()->SetTitle("Entries");
+    hLenSiz->GetYaxis()->SetTitleOffset(1.5);
+    hLenSiz->GetXaxis()->SetTitleOffset(1.1);
+    hLenSiz->SetLineColor(4);
+    hLenSiz->SetFillColor(4);
+    hLenSiz->SetLineWidth(2);
+    hLenSiz->SetFillStyle(3004);
+    hLenSiz->SetNdivisions(505);
+
+    hLenSiz->SetTitle(title);
+    hLenSiz->SetName("LENGTH/SIZE");
+
+    hLenSiz->Draw();
+
+
+    // Uncomment the following lines for MC only!
+
+    //************ Energy ************
+    /*
+   
+    TH1F *hEnergy = (TH1F*)gDirectory->Get("hEnergy");
+    
+    TCanvas* CEnergy = new TCanvas("CEnergy","Canvas",0,670,300,300);
+    CEnergy->cd();
+    CEnergy->SetBorderMode(0); 
+    //CEnergy->SetLogy();
+    gStyle->SetOptStat(110011);
+
+    TH1F* hEnergy = new TH1F("hEnergy","",100,0,1000);
+
+    CON->Draw("energy>>hEnergy","","");
+     
+    hEnergy->GetXaxis()->SetTitle("Energy [GeV]");
+    hEnergy->GetYaxis()->SetTitle("Entries");
+    hEnergy->GetYaxis()->SetTitleOffset(1.2);
+    hEnergy->GetXaxis()->SetTitleOffset(1.1);
+    hEnergy->SetLineColor(4);
+    hEnergy->SetFillColor(4);
+    hEnergy->SetLineWidth(2);
+    hEnergy->SetFillStyle(3004);
+    hEnergy->SetNdivisions(505);
+
+    hEnergy->SetTitle(title);
+    hEnergy->SetName("Energy");
+
+    hEnergy->Draw();
+    */
+
+
+    //************ MeanX and MeanY ************
+
+    // Create Canvas and set Canvas some options
+    TCanvas* Cxy = new TCanvas("Cxy","Canvas",325,670,300,300);
+    Cxy->cd();    
+    Cxy->SetBorderMode(0);   
+    gPad->SetBorderMode(0);
+    gStyle->SetOptStat(111);
+
+    TH1F* hMeanX = new TH1F("hMeanX","",100,-1.5,1.5);
+    
+    CON->Draw("meanx>>hMeanX",cut3,"");
+
+    TH1F* hMeanY = new TH1F("hMeanY","",100,-1.5,1.5);
+
+    CON->Draw("meany>>hMeanY",cut3,"");
+    
+    hMeanX->SetTitle(title); 
+    hMeanX->GetXaxis()->SetTitle("MeanX and MeanY [#circ]");
+    hMeanX->GetYaxis()->SetTitle("Entries");
+    hMeanX->GetYaxis()->SetTitleOffset(1.6);
+    hMeanX->SetFillColor(50);
+    hMeanX->SetLineColor(50);
+    hMeanX->SetLineWidth(2);
+    hMeanX->SetFillStyle(3005);
+    //hMeanX->SetMaximum(16000);
+    //hMeanX->SetMaximum(6500);
+    hMeanX->SetName("MeanX");
+ 
+    gStyle->SetStatTextColor(4);
+    hMeanY->SetFillColor(4);
+    hMeanY->SetLineColor(4);
+    hMeanY->SetFillStyle(3004);
+    hMeanY->SetLineWidth(2);
+    hMeanY->SetName("MeanY");
+
+    hMeanX->Draw();//DrawNormalized or DrawCopy are possible options
+    hMeanY->Draw("sames");//esames
+
+    // Beautify the statistic boxes!
+    gPad->Update();
+    TPaveStats* pavstat = (TPaveStats*) hMeanX->GetListOfFunctions()->FindObject("stats");
+    if(pavstat)
+    {
+	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
+	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
+	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
+	pavstat->SetTextColor(50);
+    }
+
+    gPad->Modified();
+    gPad->Update();
+
+    //delete hWL;
+    //delete hXY;
+}
+
+
+Double_t Significance(Float_t Non, Float_t Noff, Float_t alpha, Int_t type = 0)
+{
+    // Calculate significance after Li and Ma
+    
+    // Formula Nr. 17
+    if(type == 0)
+	return TMath::Sqrt(2*(Non*TMath::Log(((1+alpha)/alpha)*(Non/(Non+Noff))) + Noff*TMath::Log((1+alpha)*Noff/(Non+Noff))));
+
+    // Formula Nr. 9
+    if(type == 1)
+	return (Non-alpha*Noff)/TMath::Sqrt(alpha*(Non+Noff));
+
+}
+
+
+
+
+
+
+
+
+
