Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 4353) +++ 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(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(×); + + // Tasks +/* MReadMarsFile read4("Events"); + static_cast(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, 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 <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)); + +} + + + + + + + + +