Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 4415) +++ (revision ) @@ -1,183 +1,0 @@ -Please do not write behind the end of line tag so that it is possible -to send a Changelog file to a line printer directly. Thanks. - -For your entries please follow this style: - -_yyyy/mm/dd:_[Your Name] -_ -___*_filename: -_____-_description -_ -___*_filename2, filename3: -_____-_description -_ -_ -_ - -While an underscore is a placeholder for a white-space or an empty line. - - -*-*- END OF LINE -*-*- - - 2004/07/22: Sabrina Stark - - * /mtemp/meth/MFindStars.[h,cc] - - new fit function, incl. constant background and rotation of - 2d gauss - - initialization of MAstroCamera - - new functions: - - CorrSourcePos -> calculates a corrected - source position using the images of two stars - - DistBetweenStars -> calculates the distances between - two stars and the source - - * /mtemp/meth/MAstroCamera.[h,cc] - - TList *fCatList: list with star positions from catalog - - new functions: - - StarPosInCamera (similar to FillStarList, but - incl. aberration) - - GetCatList -> Access function to fCatList - - * /mtemp/meth/MStarLocalPos.[h,cc] - - in function SetFitValue -> background included - - - * /mtemp/meth/CorrectSrcPos.C - - macro to calculate a corrected source position using - the images of two stars, these positions are written - to an output file - - - 2004/07/20: David Paneque - - * /mtemp/mmpi/SupercutsONOFFClasses directory containing the - classes needed to optimize supercuts (and extract gamma signal) - using the ON and OFF data. These classes are listed below: - - MFRandomSplit.cc - MSupercutsCalcONOFF.cc - MTSupercutsApplied.cc - MHFindSignificanceONOFF.cc - MFindSupercutsONOFFThetaLoop.cc - MFindSupercutsONOFF.cc - - In order to add (to the compiled mars version) these classes you must add the - following line to the SUBDIRS field of the general Makefile - - mtemp/mmpi/SupercutsONOFFClasses - - - And that's it. It should compile without problems... there is no interference - with the other classes. - - I also add few macros (in the directory /mtemp/mmpi/macros/) which should help - the user to use (correctly) the mentioned classes. - - SuperCutsONOFFMacro.C - SCDynamicalSupercutsApplied.C - SCNexSignificanceVSAlphaCut.C - - - The whole analysis is set and executed by running the macro SuperCutsONOFFMacro.C. - I wrote quite some comments on it; it should be self-explanatory. - The macro SCDynamicalSupercutsApplied.C can be used to make a quick - inspection to the cuts applied to the events that survived the cuts, and the - macro SCNexSignificanceVSAlphaCut.C plots the significance and the excess events - vs the cut in the ALPHA parameter. - - I also included the directory mtemp/mmpi/asciifiles containing ascii files with - some initial supercuts values and the optimized supercuts for Mkn421 Feb 15, - that are used as input for the supercuts ON-OFF analysis. - - - 2004/06/29: Javier Rico - - * MIslands.[h,cc], MIslandCalc.[h,cc] - - Remove obsolete versions of these classes whose up-to-date ver- - sions are available at mifae/library - - - 2004/06/26: Eva Domingo - Pepe Flix - - * mifae/library/MDisp.[h,cc] - * mifae/library/MDispCalc.[h,cc] - - Added classes for evaluation of DISP method. - - 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 - - * mpisa/macros/AlphaPlot.C - -Added a new macro to plot eccess in alpha - - 2004/05/26: Nicola Galante - - * mpisa/macros/timedist3.C - - Added a new macro to plot signal time distributions - - 2004/05/26: Nicola Galante - - * mpisa/macros/production3.C - - Added a new production macro. It is based on - bootcampstandardanalisys.C - - 2004/05/23: Thomas Hengstebeck - * mberlin/macros/ScaleHillas.C , SrcCorrect.C , Convert2Matrix.C , - CutOptim.C , AlphaPlot.C , IOMkn421.h - - Added analysis macros - - - 2004/05/18: Daniel Mazin - * mmpi/macros/hillasread.C - - Removed a small bug - - * mmpi/macros/callalphaplot.C - - Added. macro to produce an alpha plot for a chosen position in the sky map - - Derotation is possible if theta and phi for each event are available - - Using OFF data to estimate significance is optional - - 2004/05/18: Sebastian Commichau - * meth/MHillasSrc.cc/h - - Added. This is a modified version, containing now - the DCA parameter. - - * meth/MDCACalc.cc/h, MDCA.cc/h - - Added. Based on MHillasCalc/MHillas, it calculates - the DCA parameter and provides a nice draw function - - 2004/05/17: Daniel Mazin - * mmpi/macros - - Added to put analysis macros there - - * mmpi/macros/callcleaning.C , cleaninghillas.C - - Added. macros to make cleaning and calculate Hillas parameters - - * mmpi/macros/mtools.C - -Added. some simple mathematical functions - - * mmpi/macros/falsesourcemethod.C, skyplot.h, hillasread.C, gridloop.C - -Added. macros to produce 2D plots of False Source Method - - 2004/04/28: Javier Rico - * mifae/Changelog - - Added. Next changes in directory mifae will be reported - exclusively in there - - - 2004/04/27: Javier Rico - * mifae - - Add new directory for IFAE temporal stuff - - * mifae/Makefile,IFAEIncl.h,IFAELinkDef.h - - Add basic stuff to run the directory - - * mifae/MDCA.[h,cc] - - Add Commichau & Starks's MDCA (provisional) - - * mifae/makeHillas.cc,makehillas.datacard - - Add program to generate hillas parameters' file Index: trunk/MagicSoft/Mars/mtemp/meth/CorrectSrcPos.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/CorrectSrcPos.C (revision 4415) +++ (revision ) @@ -1,359 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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): Luisa Sabrina Stark Schneebeli, 07/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! This macro is partially derived from findstars.C -\* ======================================================================== */ - - -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; -} - - -void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0) -{ - - // - // 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); - - - MGeomCamMagic geomcam; - MCameraDC dccam; - MStarLocalCam starcam; - - plist.AddToList(&geomcam); - plist.AddToList(&dccam); - plist.AddToList(&starcam); - - // - // Now setup the tasks and tasklist: - // --------------------------------- - // - - // Reads the trees of the root file and the analysed branches - MReadReports read; - read.AddTree("Currents"); - read.AddFile(directory+filename); // after the reading of the trees!!! - read.AddToBranchList("MReportCurrents.*"); - - MGeomApply geomapl; - TString continuoslightfile = - "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; - MCalibrateDC dccal; - dccal.SetFileName(continuoslightfile); - - - const Int_t numblind = 12; - const Short_t x[numblind] = { 8, 27, - 507, 508, 509, 510, 511, - 543, 559, 560, 561, 567}; - const TArrayS blindpixels(numblind,(Short_t*)x); - Float_t ringinterest = 100; //[mm] - Float_t tailcut = 3.5; - UInt_t integratedevents = 5; - - MFindStars findstars; - // findstars.SetBlindPixels(blindpixels); - findstars.SetRingInterest(ringinterest); - findstars.SetDCTailCut(tailcut); - findstars.SetNumIntegratedEvents(integratedevents); - findstars.SetMinuitPrintOutLevel(0); - - findstars.SetSourceRaH(11); // Mrk421 - findstars.SetSourceRaM(4); - findstars.SetSourceRaS(26.0); - findstars.SetSourceDecD(38); - findstars.SetSourceDecM(12); - findstars.SetSourceDecS(36.0); - - - // prints - MPrint pdc("MCameraDC"); - MPrint pstar("MStarLocalCam"); - - tlist.AddToList(&geomapl); - tlist.AddToList(&read); - tlist.AddToList(&dccal); - // tlist.AddToList(&pdc, "Currents"); - tlist.AddToList(&findstars, "Currents"); - // tlist.AddToList(&pstar, "Currents"); - - // - // Create and setup the eventloop - // - MEvtLoop evtloop; - evtloop.SetParList(&plist); - - const Int_t numpar = 3; - Float_t starpos[10*numpar]; - Float_t starposcat[10*numpar]; - - if (numEvents > 0) - { - if (!evtloop.Eventloop(numEvents)) - return; - } - else - { - if (!evtloop.PreProcess()) - return; - - MHCamera display(geomcam); - display.SetPrettyPalette(); - - TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600); - - TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2); - h0->GetXaxis()->SetTitle("(deg)"); - h0->GetYaxis()->SetTitle("(deg)"); - h0->Draw(); - - display.Draw(); - gPad->cd(1); - //display.Draw("same"); - starcam.Draw(); - - UShort_t year; - Byte_t month, day; - - MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents"); - - TList *starcatalog = new TList; - TList *stardc = new TList; - TVector3 *res = new TVector3; - - Int_t maxdist = 100; - UInt_t numevents=0; - Int_t numEvt=0; - const Int_t maxnumEvt = 10000; - Int_t col = 2; - const Double_t kDeg2mm = 189/0.6; - - Double_t resx[maxnumEvt]; - Double_t resy[maxnumEvt]; - Double_t time[maxnumEvt]; - - -// -// Find in the file name the date, run and project name - - Size_t pos = filename.Last('/'); - TString iRun = TString(filename(pos+24,5)); - TString iYear = TString(filename(pos+4,4)); - TString iMonth = TString(filename(pos+9,2)); - TString iDay = TString(filename(pos+12,2)); - - TString iHour = TString(filename(pos+15,2)); - TString iMin = TString(filename(pos+18,2)); - TString iSec = TString(filename(pos+21,2)); - - Size_t poslast = filename.Last('.'); - Size_t posfirst = poslast-1; - while (filename[posfirst] != '_') - posfirst--; - - TString iSource = TString(filename(posfirst+1,poslast-posfirst-1)); - char str[100]; - char str2[100]; - sprintf(str,"Date %s %s %s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); - sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); - - UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec; - UShort_t evtmsec; - - ofstream fposout1; - // fposout1.open (str2, ofstream::app); - fposout1.open (str2, ofstream::trunc); - - while (tlist.Process()) - { - if (numevents >= integratedevents) - { - evttime->GetTime(evthour, evtmin, evtsec,evtmsec); - evttime->Print(); - evthour = evttime->Hour(); - evtmin = evttime->Min(); - evtsec = evttime->Sec(); - - // display.SetCamContent(dccam); - display.SetCamContent(findstars.GetDisplay()); - gPad->Modified(); - gPad->Update(); - Int_t numStars = starcam.GetNumStars(); - starcam.Print("bkmaxposrotsizechi"); - - res = findstars.GetCorrSrcPos(); - cout<<"Corrected source position = "<X()<<" , "<Y()<X(); - resy[numEvt] = res->Y(); - time[numEvt] = numEvt; - numEvt ++; - - // write the corrected source position and the star positions to file - if (!fposout1) - cout << "Could not open output file to write the center of the fit." << endl; - else{ - for ( Int_t i = 0; i < numStars; i++) - { - if (i == 0){ - fposout1 << str<<" "<X()<< " " << res->Y(); - } - fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY(); - } - } - fposout1 << endl; - numevents = 0; - - // Remove the comments if you want to go through the file - // event-by-event: - //if (!HandleInput()) - // break; - } - numevents++; - } - fposout1.close(); - - // Print source position - gROOT->Reset(); - gStyle->SetCanvasColor(0); - gStyle->SetCanvasBorderMode(0); - gStyle->SetPadBorderMode(0); - gStyle->SetFrameBorderMode(0); - gStyle->SetOptStat(00000000); - - TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200); - c1->Divide(1,3); - - TMath math; - - Double_t minresx, maxresx; - minresx = resx[math.LocMin(numEvt,resx)]; - maxresx = resx[math.LocMax(numEvt,resx)]; - - Double_t minresy, maxresy; - minresy = resy[math.LocMin(numEvt,resy)]; - maxresy = resy[math.LocMax(numEvt,resy)]; - - Double_t minres, maxres; - minres = math.Min(minresx,minresy); - maxres = math.Max(maxresx,maxresy); - - Double_t diff; - diff = maxres - minres; - diff = 0.1*diff; - minres = minres - diff; - maxres = maxres + diff; - - Double_t mintime, maxtime; - mintime = time[math.LocMin(numEvt,time)]; - maxtime = time[math.LocMax(numEvt,time)]; - - c1->cd(1); - TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres); - h1->GetXaxis()->SetTitle("(evt.nr)"); - h1->GetYaxis()->SetTitle("(mm)"); - h1->Draw(); - - TGraph *grtimeevolx = new TGraph(numEvt,time,resx); - grtimeevolx->SetMarkerColor(215); - grtimeevolx->SetMarkerStyle(20); - grtimeevolx->SetMarkerSize (0.4); - grtimeevolx->Draw("P"); - - TGraph *grtimeevoly = new TGraph(numEvt,time,resy); - grtimeevoly->SetMarkerColor(6); - grtimeevoly->SetMarkerStyle(24); - grtimeevoly->SetMarkerSize (0.4); - grtimeevoly->Draw("P"); - - legxy = new TLegend(0.8,0.85,0.95,0.95); - legxy.SetTextSize(0.03); - legxy.AddEntry(grtimeevolx,"src pos x","P"); - legxy.AddEntry(grtimeevoly,"src pos y","P"); - legxy.Draw(); - - c1->cd(2); - TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx); - h2->GetXaxis()->SetTitle("(evt.nr)"); - h2->GetYaxis()->SetTitle("(mm)"); - h2->Draw(); - - TGraph *grtimeevolx = new TGraph(numEvt,time,resx); - grtimeevolx->SetMarkerColor(215); - grtimeevolx->SetMarkerStyle(20); - grtimeevolx->SetMarkerSize (0.4); - grtimeevolx->Draw("P"); - - legxy = new TLegend(0.8,0.85,0.95,0.95); - legxy.SetTextSize(0.03); - legxy.AddEntry(grtimeevolx,"src pos x","P"); - legxy.Draw(); - - c1->cd(3); - TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy); - h3->GetXaxis()->SetTitle("(evt.nr)"); - h3->GetYaxis()->SetTitle("(mm)"); - h3->Draw(); - - TGraph *grtimeevoly = new TGraph(numEvt,time,resy); - grtimeevoly->SetMarkerColor(6); - grtimeevoly->SetMarkerStyle(24); - grtimeevoly->SetMarkerSize (0.4); - grtimeevoly->Draw("P"); - - legxy = new TLegend(0.8,0.85,0.95,0.95); - legxy.SetTextSize(0.03); - legxy.AddEntry(grtimeevoly,"src pos y","P"); - legxy.Draw(); - - evtloop.PostProcess(); - } - - tlist.PrintStatistics(); - -} - - Index: trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C (revision 4416) +++ trunk/MagicSoft/Mars/mtemp/meth/macros/CorrectSrcPos.C (revision 4416) @@ -0,0 +1,359 @@ +/* ======================================================================== *\ +! +! * +! * 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): Luisa Sabrina Stark Schneebeli, 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! This macro is partially derived from findstars.C +\* ======================================================================== */ + + +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; +} + + +void CorrectSrcPos(const TString filename="dc_2004_04_21_22_*_Mrk421.root", const TString directory="/scratch/Period016/cacodata/2004_04_22/", const UInt_t numEvents = 0) +{ + + // + // 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); + + + MGeomCamMagic geomcam; + MCameraDC dccam; + MStarLocalCam starcam; + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&starcam); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddFile(directory+filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + + MGeomApply geomapl; + TString continuoslightfile = + "/scratch/Period016/cacodata/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; + MCalibrateDC dccal; + dccal.SetFileName(continuoslightfile); + + + const Int_t numblind = 12; + const Short_t x[numblind] = { 8, 27, + 507, 508, 509, 510, 511, + 543, 559, 560, 561, 567}; + const TArrayS blindpixels(numblind,(Short_t*)x); + Float_t ringinterest = 100; //[mm] + Float_t tailcut = 3.5; + UInt_t integratedevents = 5; + + MFindStars findstars; + // findstars.SetBlindPixels(blindpixels); + findstars.SetRingInterest(ringinterest); + findstars.SetDCTailCut(tailcut); + findstars.SetNumIntegratedEvents(integratedevents); + findstars.SetMinuitPrintOutLevel(0); + + findstars.SetSourceRaH(11); // Mrk421 + findstars.SetSourceRaM(4); + findstars.SetSourceRaS(26.0); + findstars.SetSourceDecD(38); + findstars.SetSourceDecM(12); + findstars.SetSourceDecS(36.0); + + + // prints + MPrint pdc("MCameraDC"); + MPrint pstar("MStarLocalCam"); + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&dccal); + // tlist.AddToList(&pdc, "Currents"); + tlist.AddToList(&findstars, "Currents"); + // tlist.AddToList(&pstar, "Currents"); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + const Int_t numpar = 3; + Float_t starpos[10*numpar]; + Float_t starposcat[10*numpar]; + + if (numEvents > 0) + { + if (!evtloop.Eventloop(numEvents)) + return; + } + else + { + if (!evtloop.PreProcess()) + return; + + MHCamera display(geomcam); + display.SetPrettyPalette(); + + TCanvas *c = new TCanvas("c", "FindStars", 100, 100,600,600); + + TH2D *h0 = new TH2D("h0","",1,-2,2,1,-2,2); + h0->GetXaxis()->SetTitle("(deg)"); + h0->GetYaxis()->SetTitle("(deg)"); + h0->Draw(); + + display.Draw(); + gPad->cd(1); + //display.Draw("same"); + starcam.Draw(); + + UShort_t year; + Byte_t month, day; + + MTime *evttime = (MTime*)plist->FindObject("MTimeCurrents"); + + TList *starcatalog = new TList; + TList *stardc = new TList; + TVector3 *res = new TVector3; + + Int_t maxdist = 100; + UInt_t numevents=0; + Int_t numEvt=0; + const Int_t maxnumEvt = 10000; + Int_t col = 2; + const Double_t kDeg2mm = 189/0.6; + + Double_t resx[maxnumEvt]; + Double_t resy[maxnumEvt]; + Double_t time[maxnumEvt]; + + +// +// Find in the file name the date, run and project name + + Size_t pos = filename.Last('/'); + TString iRun = TString(filename(pos+24,5)); + TString iYear = TString(filename(pos+4,4)); + TString iMonth = TString(filename(pos+9,2)); + TString iDay = TString(filename(pos+12,2)); + + TString iHour = TString(filename(pos+15,2)); + TString iMin = TString(filename(pos+18,2)); + TString iSec = TString(filename(pos+21,2)); + + Size_t poslast = filename.Last('.'); + Size_t posfirst = poslast-1; + while (filename[posfirst] != '_') + posfirst--; + + TString iSource = TString(filename(posfirst+1,poslast-posfirst-1)); + char str[100]; + char str2[100]; + sprintf(str,"Date %s %s %s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + sprintf(str2,"CorrectSrcPos_%s%s%s_%s_%s.txt",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data()); + + UInt_t evtyear, evtmonth, evtday, evthour, evtmin, evtsec; + UShort_t evtmsec; + + ofstream fposout1; + // fposout1.open (str2, ofstream::app); + fposout1.open (str2, ofstream::trunc); + + while (tlist.Process()) + { + if (numevents >= integratedevents) + { + evttime->GetTime(evthour, evtmin, evtsec,evtmsec); + evttime->Print(); + evthour = evttime->Hour(); + evtmin = evttime->Min(); + evtsec = evttime->Sec(); + + // display.SetCamContent(dccam); + display.SetCamContent(findstars.GetDisplay()); + gPad->Modified(); + gPad->Update(); + Int_t numStars = starcam.GetNumStars(); + starcam.Print("bkmaxposrotsizechi"); + + res = findstars.GetCorrSrcPos(); + cout<<"Corrected source position = "<X()<<" , "<Y()<X(); + resy[numEvt] = res->Y(); + time[numEvt] = numEvt; + numEvt ++; + + // write the corrected source position and the star positions to file + if (!fposout1) + cout << "Could not open output file to write the center of the fit." << endl; + else{ + for ( Int_t i = 0; i < numStars; i++) + { + if (i == 0){ + fposout1 << str<<" "<X()<< " " << res->Y(); + } + fposout1 <<" "<< starcam[i].GetMeanX()<< " " << starcam[i].GetMeanY(); + } + } + fposout1 << endl; + numevents = 0; + + // Remove the comments if you want to go through the file + // event-by-event: + //if (!HandleInput()) + // break; + } + numevents++; + } + fposout1.close(); + + // Print source position + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptStat(00000000); + + TCanvas *c1 = new TCanvas("c1", "SourcePosition", 200, 200,600,1200); + c1->Divide(1,3); + + TMath math; + + Double_t minresx, maxresx; + minresx = resx[math.LocMin(numEvt,resx)]; + maxresx = resx[math.LocMax(numEvt,resx)]; + + Double_t minresy, maxresy; + minresy = resy[math.LocMin(numEvt,resy)]; + maxresy = resy[math.LocMax(numEvt,resy)]; + + Double_t minres, maxres; + minres = math.Min(minresx,minresy); + maxres = math.Max(maxresx,maxresy); + + Double_t diff; + diff = maxres - minres; + diff = 0.1*diff; + minres = minres - diff; + maxres = maxres + diff; + + Double_t mintime, maxtime; + mintime = time[math.LocMin(numEvt,time)]; + maxtime = time[math.LocMax(numEvt,time)]; + + c1->cd(1); + TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minres,maxres); + h1->GetXaxis()->SetTitle("(evt.nr)"); + h1->GetYaxis()->SetTitle("(mm)"); + h1->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + c1->cd(2); + TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minresx,maxresx); + h2->GetXaxis()->SetTitle("(evt.nr)"); + h2->GetYaxis()->SetTitle("(mm)"); + h2->Draw(); + + TGraph *grtimeevolx = new TGraph(numEvt,time,resx); + grtimeevolx->SetMarkerColor(215); + grtimeevolx->SetMarkerStyle(20); + grtimeevolx->SetMarkerSize (0.4); + grtimeevolx->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevolx,"src pos x","P"); + legxy.Draw(); + + c1->cd(3); + TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minresy,maxresy); + h3->GetXaxis()->SetTitle("(evt.nr)"); + h3->GetYaxis()->SetTitle("(mm)"); + h3->Draw(); + + TGraph *grtimeevoly = new TGraph(numEvt,time,resy); + grtimeevoly->SetMarkerColor(6); + grtimeevoly->SetMarkerStyle(24); + grtimeevoly->SetMarkerSize (0.4); + grtimeevoly->Draw("P"); + + legxy = new TLegend(0.8,0.85,0.95,0.95); + legxy.SetTextSize(0.03); + legxy.AddEntry(grtimeevoly,"src pos y","P"); + legxy.Draw(); + + evtloop.PostProcess(); + } + + tlist.PrintStatistics(); + +} + +