Index: trunk/MagicSoft/Mars/mtemp/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/Changelog (revision 3856) +++ trunk/MagicSoft/Mars/mtemp/Changelog (revision 3857) @@ -18,2 +18,15 @@ -*-*- END OF LINE -*-*- + + 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/mifae/IFAEIncl.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/IFAEIncl.h (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/IFAEIncl.h (revision 3857) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mtemp/mifae/IFAELinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/IFAELinkDef.h (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/IFAELinkDef.h (revision 3857) @@ -0,0 +1,9 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MDCA+; + +#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/MDCA.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MDCA.cc (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/MDCA.cc (revision 3857) @@ -0,0 +1,261 @@ + +///////////////////////////////////////////////////////////////////////////// +// // +// MDCA.cc // +// Container to store Hillas parameters and DCA stuff // +// // +// Author(s): S.C. Commichau, L.S. Stark, 7/2003 // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MDCA.h" + +using namespace std; + +ClassImp(MDCA); + +// Default constructor +MDCA::MDCA(const char *name, const char *title) +{ + fName = name ? name : "MDCA"; + fTitle = title ? title : "Storage container for Hillas parameters and DCA of one event"; + + Reset(); + + fEllipse = new TEllipse; + fRefCircle = new TEllipse; + + fLineL = new TLine; + fLineW = new TLine; + fLineX = new TLine; + fLineY = new TLine; + fLineDCA = new TLine; + fLineMean = new TLine; + +} + +// Destructor: Deletes ellipse and lines if they do exist + +MDCA::~MDCA() +{ + Clear(); +} + +// Initialize parameters with default values + +void MDCA::Reset() +{ + + fLength = -1; + fWidth = -1; + fDelta0 = 0; + fMeanX = 0; + fMeanY = 0; + fDelta1 = 0; + fDCA = -1000; + fX1W = 0; + fY1W = 0; + fX2W = 0; + fY2W = 0; + fX1L = 0; + fY1L = 0; + fX2L = 0; + fY2L = 0; + fXDCA = 0; + fYDCA = 0; + +} + + +// Print parameters to *fLog + +void MDCA::Print(Option_t *) const +{ + Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg; + + if (atg<0) + atg += 180; + + *fLog << all; + *fLog << "Basic Image Parameters (" << GetName() << ")" << endl; + *fLog << " - Length [mm] = " << fLength << endl; + *fLog << " - Width [mm] = " << fWidth << endl; + *fLog << " - Delta0 [deg] = " << fDelta0*kRad2Deg << endl; + *fLog << " - Meanx [mm] = " << fMeanX << endl; + *fLog << " - Meany [mm] = " << fMeanY << endl; + *fLog << " - atg(y/x) [deg] = " << atg << endl; + *fLog << " - DCA [mm] = " << fDCA << endl; + *fLog << " - Delta1 [deg] = " << fDelta1*kRad2Deg << endl; + + +} + +void MDCA::Paint(Option_t *opt) +{ + Clear(); + + if (fLength<=0 || fWidth<=0) //fLength<0 || fWidth<0 doesn't look nice... + return; //We get a circle! + + // Length line + fLineL = new TLine(fX1L, fY1L, fX2L, fY2L); + fLineL->SetLineWidth(2); + fLineL->SetLineColor(2); + fLineL->Draw(); + + // Width line + fLineW = new TLine(fX1W, fY1W, fX2W, fY2W); + fLineW->SetLineWidth(2); + fLineW->SetLineColor(2); + fLineW->Draw(); + + // Coordinate system + fLineX = new TLine(-600,fYRef,600,fYRef); + fLineY = new TLine(fXRef,-600,fXRef,600); + fLineX->SetLineWidth(1); + fLineX->SetLineColor(108); + fLineY->SetLineWidth(1); + fLineY->SetLineColor(108); + fLineX->Draw(); + fLineY->Draw(); + + // DCA line + fLineDCA = new TLine(fXRef,fYRef,fXDCA+fXRef,fYDCA+fYRef); + fLineDCA->SetLineWidth(2); + fLineDCA->SetLineColor(2); + fLineDCA->Draw(); + + // COG line + fLineMean = new TLine(fXRef,fYRef,fMeanX,fMeanY); + fLineMean->SetLineWidth(2); + fLineMean->SetLineColor(2); + fLineMean->Draw(); + + // Reference point marker + fRefCircle = new TEllipse(fXRef, fYRef, 10, 10, 0, 360, 0); + fRefCircle->SetLineColor(108); + fRefCircle->SetFillColor(108); + fRefCircle->Draw(); + + // Hillas ellipse + fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta0*kRad2Deg+180); + fEllipse->SetLineWidth(2); + fEllipse->SetLineColor(2); + fEllipse->Draw(); + +} + + +// If an ellipse and lines exist they will be deleted + +void MDCA::Clear(Option_t *) +{ + if (!fEllipse && !fRefCircle && !fLineL && !fLineW && !fLineX && !fLineY && !fLineDCA && !fLineMean) + return; + + delete fEllipse; + delete fRefCircle; + delete fLineL; + delete fLineW; + delete fLineX; + delete fLineY; + delete fLineDCA; + delete fLineMean; + + fLineL = NULL; + fLineX = NULL; + fLineY = NULL; + fLineW = NULL; + fLineDCA = NULL; + fLineMean = NULL; + fEllipse = NULL; + fRefCircle = NULL; +} + + +Int_t MDCA::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil) +{ + // Get basic Hillas parameters from MHillas + fDelta0 = hil.GetDelta(); + fMeanX = hil.GetMeanX(); + fMeanY = hil.GetMeanY(); + fLength = hil.GetLength(); + fWidth = hil.GetWidth(); + + // The Length Line - rotation and shift + fX1L = - (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm] + fY1L = - (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm] + fX2L = (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm] + fY2L = (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm] + + // The Width Line - rotation and shift + fX1W = (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm] + fY1W = - (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm] + fX2W = - (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm] + fY2W = (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm] + + // Vector of orientation of the shower axis + fr1 = fX2L-fX1L; + fr2 = fY2L-fY1L; + + // Determine parameters to calculate coordinates of the DCA vector + flambda = (fr1*(fXRef-fMeanX) + fr2*(fYRef-fMeanY))/(fr1*fr1 + fr2*fr2); + fmu = (fMeanY-fYRef)/fr1 + flambda*fr2/fr1; + + // Components of the DCA vector + fXDCA = -fmu*fr2; + fYDCA = fmu*fr1; + + // Components of vector going from intersection point of the DCA vector + // with the shower axis to the COG + fd1 = fMeanX + fmu*fr2 - fXRef; + fd2 = fMeanY - fmu*fr1 - fYRef; + + // Calculate DCA value + fDCA = sqrt(fXDCA*fXDCA + fYDCA*fYDCA); + + // Calculate angle of the shower axis with respect to the x-axis + fDelta1 = acos(fd1/sqrt(fd1*fd1 + fd2*fd2)); //botch + + // Calculate angle of the shower axis with respect to the y-axis + //fDelta1 = acos(fd2/sqrt(fd1*fd1 + fd2*fd2)); + + // Determine the correct sign of the DCA (cross product of DCA vector and + // vector going from the intersection point of the DCA vector with the shower axis + // to the COG) + if((fmu*(-fr2*(fMeanY-fYRef)-fr1*(fMeanX-fXRef)))<0) + fDCA = -fDCA; + + gRandom->Rannor(gx,gy); + gx = fabs(gx); + + // This is nice but does not remove the systematics in the profile plot... + //if(((1-0.6*gx)*(180-kRad2Deg*fDelta1)>120) || ((1-0.6*gx)*kRad2Deg*fDelta1>120)) + // fDCA = -1; + + // Enlarge the interval of Delta correctly... + if((fMeanY-fYRef-fmu*fr1)<0) + fDelta1 = TwoPi-fDelta1; + + // Enlarge the interval of Delta correctly... (Delta with respect to the y-axis) + // if(-(fMeanX-fXRef+fmu*fr2)<0) + // fDelta1 = TwoPi-fDelta1; + + // This has to be improved... + if(fr1 == 0 || fr2 == 0 || (fr1*fr1+fr2*fr2) == 0 || sqrt(fd1*fd1+fd2*fd2) == 0) + fDCA = -fDCA; + + SetReadyToSave(); + + return 0; +} + +void MDCA::SetRefPoint(const Float_t fXRef0, const Float_t fYRef0) +{ + fXRef = fXRef0; + fYRef = fYRef0; +} + + + + Index: trunk/MagicSoft/Mars/mtemp/mifae/MDCA.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MDCA.h (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/MDCA.h (revision 3857) @@ -0,0 +1,118 @@ +#ifndef MARS_MDCA +#define MARS_MDCA + +///////////////////////////////////////////////////////////////////////////// +// // +// MDCA.cc // +// Container to store the DCA stuff // +// // +// Author(s): S.C. Commichau, L.S. Stark, 7/2003 // +// // +// // +///////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MGeomCam.h" +#include "MGeomPix.h" +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" +#include "MHillas.h" + +// Add offsets to the Width and Length to draw nice lines +#define OffsetW 20.0 +#define OffsetL 440.0 + +#define kmm2Deg 0.00317460317 +#define TwoPi 6.28318530717958623 + +class TArrayF; +class TEllipse; +class MGeoCam; +class MCerPhotEvt; +class MHillas; + +class MDCA : public MParContainer +{ + private: + + Float_t fLength; // [mm] major axis of ellipse + Float_t fWidth; // [mm] minor axis of ellipse + Float_t fDelta0; // [rad] angle of major axis with x-axis (-Pi/2 .. Pi/2) + Float_t fDelta1; // [rad] my angle of major axis with x-axis (0 .. 2*Pi) + Float_t fMeanX; // [mm] x-coordinate of center of ellipse + Float_t fMeanY; // [mm] y-coordinate of center of ellipse + Float_t fX1W; // [mm] x-coordinate of 1st point of Lline + Float_t fY1W; // [mm] y-coordinate of 1st point of Lline + Float_t fX2W; // [mm] x-coordinate of 2nd point of Lline + Float_t fY2W; // [mm] y-coordinate of 2nd point of Lline + Float_t fX1L; // [mm] x-coordinate of 1st point of Wline + Float_t fY1L; // [mm] y-coordinate of 1st point of Wline + Float_t fX2L; // [mm] x-coordinate of 2nd point of Wline + Float_t fY2L; // [mm] y-coordinate of 2nd point of Wline + Double_t fDCA; // [mm] Distance of Cloasest Approach + Float_t fXDCA; // [mm] x-coordinate of 2nd point of the DCA-line + Float_t fYDCA; // [mm] y-coordinate of 2nd point of the DCA-line + Float_t fXRef; // [mm] x-coordinate of reference point + Float_t fYRef; // [mm] y-coordinate of reference point + Float_t fmu; + Float_t flambda; + Float_t fr1, fr2; // [mm] Coordinates of the orientation vector of the shower axis + Float_t fd1, fd2; // [mm] Coordinates of the DCA vector + Float_t gx, gy; + + TEllipse *fEllipse; // Graphical object to draw the ellipse, ROOT! + TEllipse *fRefCircle; // To draw reference point + + TLine *fLineL; // Shower axis + TLine *fLineW; // Line perpendicular to the shower axis + TLine *fLineX; // x-axis of the coordinate system + TLine *fLineY; // y-axis of the coordinate system + TLine *fLineDCA; // DCA line + TLine *fLineMean; // Line to COG of the shower + + public: + MDCA(const char *name=NULL, const char *title=NULL); + + ~MDCA(); + + void Reset(); + + Int_t Calc(const MGeomCam &geom, const MCerPhotEvt &pix, const MHillas &hil); + + void Print(Option_t *opt=NULL) const; + + void Paint(Option_t *opt=NULL); + + void Clear(Option_t *opt=NULL); + + void SetRefPoint(const Float_t fXRef0, const Float_t fYRef0); + + Float_t GetDCA() const { return fDCA; } + + Float_t GetDelta() const { return fDelta1; } + + ClassDef(MDCA, 1) + +}; + +#endif + + + + + Index: trunk/MagicSoft/Mars/mtemp/mifae/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Makefile (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/Makefile (revision 3857) @@ -0,0 +1,72 @@ +################################################################## +# +# makefile +# +# for the MARS IFAE library +# +################################################################## +# @maintitle + +# @code + +# +# please change all system depend values in the +# config.mk.${OSTYPE} file +# +# +include ../../Makefile.conf.$(OSTYPE) +include ../../Makefile.conf.general + +# +CINT = IFAE +# + +PROGRAMS = makeHillas +SOLIB = ../../libmars.so +LIB = mifae.a + +#------------------------------------------------------------------------------ + +INCLUDES = -I. \ + -I../../mbase \ + -I../../mpedestal \ + -I../../mbadpixels \ + -I../../mfileio \ + -I../../mraw \ + -I../../manalysis \ + -I../../mgui \ + -I../../mgeom \ + -I../../msignal \ + -I../../mcalib \ + -I../../mfilter \ + -I../../mhbase \ + -I../../mimage \ + -I../../mpointing \ + -I../../mastro + + +.SUFFIXES: .c .cc .h .o + +SRCFILES = \ + MDCA.cc +# MSrcRotate.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ +all: $(LIB) $(PROGRAMS) + @echo " Done. " + @echo " " + +include ../../Makefile.rules + + +$(PROGRAMS): $(SOLIB) $(LIB) $(PROGRAMS:=.o) + @echo " Linking $@ ..." + $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(SOLIB) $(LIB) $@.o $(MARS_LIB) -o $@ + +mrproper: clean rmbak + +# @endcode Index: trunk/MagicSoft/Mars/mtemp/mifae/makeHillas.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/makeHillas.cc (revision 3857) +++ trunk/MagicSoft/Mars/mtemp/mifae/makeHillas.cc (revision 3857) @@ -0,0 +1,892 @@ +/*********************************/ +/* Compute the hillas parameters */ +/*********************************/ + +#include "TString.h" +#include "TArrayS.h" + +#include "MParList.h" +#include "MTaskList.h" +#include "MPedestalCam.h" +#include "MBadPixelsCam.h" +#include "MReadMarsFile.h" +#include "MGeomApply.h" +#include "MPedCalcPedRun.h" +#include "MEvtLoop.h" +#include "MGeomCamMagic.h" +#include "MExtractedSignalCam.h" +#include "MCalibrationChargeCam.h" +#include "MHCalibrationChargeCam.h" +#include "MHCalibrationRelTimeCam.h" +#include "MExtractSignal.h" +#include "MCalibrationChargeCalc.h" +#include "MFCosmics.h" +#include "MContinue.h" +#include "MFillH.h" +#include "MLog.h" +#include "MCerPhotEvt.h" +#include "MPedPhotCam.h" +#include "MCalibrate.h" +#include "MPedPhotCalc.h" +#include "MHillas.h" +#include "MRawRunHeader.h" +#include "MSrcPosCam.h" +#include "MBlindPixelCalc.h" +#include "MImgCleanStd.h" +#include "MHillasSrcCalc.h" +#include "MHillasCalc.h" +#include "MWriteRootFile.h" +#include "MProgressBar.h" +#include "MArgs.h" + +#include +#include + +using namespace std; + +Bool_t readDatacards(TString& filename); +void makeHillas(); + +// initial and final time slices to be used in signal extraction +const Byte_t hifirst = 1; +const Byte_t hilast = 14; +const Byte_t lofirst = 4; +const Byte_t lolast = 13; + +// declaration of variables read from datacards +TString outname; +TString pedname; +TString calname; +TString* datafile; +ULong_t nmaxevents=999999999; +Short_t calflag=1; +Float_t lcore = 3.0; +Float_t ltail = 1.5; +Int_t nfiles = 0; + +const TString defaultcard="input.datacard"; +static void Usage() +{ + gLog <" << endl << endl; + gLog << " : datacards file name (dafault input.datacards)" << endl; + gLog << " -?/-h: This help" << endl << endl; +} + +int main(int argc, char **argv) +{ + // evaluate arguments + MArgs arg(argc, argv); + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return -1; + } + + TString datacard = arg.GetArgumentStr(0); + if(!datacard.Length()) + datacard = defaultcard; + + if(!readDatacards(datacard)) + { + cout << "Error reading datacards. Stoping" << endl; + return -1; + } + makeHillas(); +} + +/*************************************************************/ +void makeHillas() +{ +#if 0 + MStatusDisplay *display = new MStatusDisplay; + display->SetUpdateTime(3000); + display->Resize(850,700); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); +#endif + + /************************************/ + /* FIRST LOOP: PEDESTAL COMPUTATION */ + /************************************/ + + MParList plist1; + MTaskList tlist1; + plist1.AddToList(&tlist1); + + // containers + MPedestalCam pedcam; + MBadPixelsCam badcam; + // badcam.AsciiRead("badpixels.dat"); + + plist1.AddToList(&pedcam); + plist1.AddToList(&badcam); + + //tasks + MReadMarsFile read("Events", pedname); + read.DisableAutoScheme(); + + MGeomApply geomapl; + MPedCalcPedRun pedcalc; + + tlist1.AddToList(&read); + tlist1.AddToList(&geomapl); + tlist1.AddToList(&pedcalc); + + // Create and setup the eventloop + MEvtLoop pedloop; + pedloop.SetParList(&plist1); + //pedloop.SetDisplay(display); + + cout << "*************************" << endl; + cout << "** COMPUTING PEDESTALS **" << endl; + cout << "*************************" << endl; + + if (!pedloop.Eventloop()) + return; + + tlist1.PrintStatistics(); + + /*****************************/ + /* SECOND LOOP: CALIBRATION */ + /*****************************/ + + MParList plist2; + MTaskList tlist2; + plist2.AddToList(&tlist2); + plist2.AddToList(&pedcam); + plist2.AddToList(&badcam); + + // new containers + MGeomCamMagic geomcam; + MExtractedSignalCam sigcam; + MCalibrationChargeCam calcam; + MHCalibrationChargeCam histcharge; + MHCalibrationRelTimeCam histtime; + + plist2.AddToList(&geomcam); + plist2.AddToList(&sigcam); + plist2.AddToList(&calcam); + plist2.AddToList(&histcharge); + + //tasks + MReadMarsFile read2("Events", calname); + read2.DisableAutoScheme(); + + MExtractSignal sigcalc; + sigcalc.SetRange(hifirst,hilast,lofirst,lolast); + MCalibrationChargeCalc calcalc; + MFCosmics cosmics; + MContinue cont(&cosmics); + + + MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam"); + + tlist2.AddToList(&read2); + tlist2.AddToList(&geomapl); + tlist2.AddToList(&sigcalc); + tlist2.AddToList(&cont); + tlist2.AddToList(&fillcam); + tlist2.AddToList(&calcalc); + + // Create and setup the eventloop + MEvtLoop calloop; + calloop.SetParList(&plist2); + //calloop.SetDisplay(display); + + cout << "***************************" << endl; + cout << "** COMPUTING CALIBRATION **" << endl; + cout << "***************************" << endl; + + if (!calloop.Eventloop()) + return; + + tlist2.PrintStatistics(); + + MLog gauglog; + gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1); + calcam.SetLogStream(&gauglog); + calcam.Print(); + calcam.SetLogStream(&gLog); +#if 0 + // Create histograms to display + MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges"); + MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges"); + MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit"); + MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas"); + MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge"); + MHCamera disp6 (geomcam, "Cal;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)"); + MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)"); + MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)"); + MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)"); + MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)"); + MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)"); + MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)"); + MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)"); + MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)"); + MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded"); + MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted"); + MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp18 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain"); + MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain"); + MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain"); + MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain"); + MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); + MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration"); + MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration"); + MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration"); + + MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals"); + MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS"); + + MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times"); + MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times"); + MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit"); + MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels"); + + MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times"); + MHCamera disp34 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times"); + + // Fitted charge means and sigmas + disp1.SetCamContent(calcam, 0); + disp1.SetCamError( calcam, 1); + disp2.SetCamContent(calcam, 2); + disp2.SetCamError( calcam, 3); + + // Fit probabilities + disp3.SetCamContent(calcam, 4); + + // Reduced Sigmas and reduced sigmas per charge + disp4.SetCamContent(calcam, 5); + disp4.SetCamError( calcam, 6); + disp5.SetCamContent(calcam, 7); + disp5.SetCamError( calcam, 8); + + // F-Factor Method + disp6.SetCamContent(calcam, 9); + disp6.SetCamError( calcam, 10); + disp7.SetCamContent(calcam, 11); + disp7.SetCamError( calcam, 12); + disp8.SetCamContent(calcam, 13); + disp8.SetCamError( calcam, 14); + + // Blind Pixel Method + disp9.SetCamContent(calcam, 15); + disp9.SetCamError( calcam, 16); + disp10.SetCamContent(calcam,17); + disp10.SetCamError( calcam,18); + disp11.SetCamContent(calcam,19); + disp11.SetCamError( calcam,20); + + // PIN Diode Method + disp12.SetCamContent(calcam,21); + disp12.SetCamError( calcam,22); + disp13.SetCamContent(calcam,23); + disp13.SetCamError( calcam,24); + disp14.SetCamContent(calcam,25); + disp14.SetCamError( calcam,26); + + // Pixels with defects + disp15.SetCamContent(calcam,27); + disp16.SetCamContent(calcam,28); + disp17.SetCamContent(badcam,9); + disp18.SetCamContent(badcam,16); + disp19.SetCamContent(badcam,15); + disp20.SetCamContent(calcam,29); + disp21.SetCamContent(calcam,30); + + // Lo Gain calibration + disp22.SetCamContent(calcam,31); + + // Valid flags + disp23.SetCamContent(calcam,32); + disp24.SetCamContent(calcam,33); + disp25.SetCamContent(calcam,34); + + // Pedestals + disp26.SetCamContent(calcam,35); + disp26.SetCamError( calcam,36); + disp27.SetCamContent(calcam,37); + disp27.SetCamError( calcam,38); + + // Relative Times + disp28.SetCamContent(histtime,0); + disp28.SetCamError( histtime,1); + disp29.SetCamContent(histtime,2); + disp29.SetCamError( histtime,3); + disp30.SetCamContent(histtime,4); + disp31.SetCamContent(histtime,5); + disp32.SetCamContent(histtime,6); + + // Absolute Times + disp33.SetCamContent(calcam,39); + disp33.SetCamError( calcam,40); + disp34.SetCamContent(calcam,41); + + disp1.SetYTitle("Charge [FADC units]"); + disp2.SetYTitle("\\sigma_{Charge} [FADC units]"); + disp3.SetYTitle("P_{Charge} [1]"); + + disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]"); + disp5.SetYTitle("Reduced Sigma / Mean Charge [1]"); + + disp6.SetYTitle("Nr. Photo-electrons [1]"); + disp7.SetYTitle("Conversion Factor [Ph/FADC Count]"); + disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] "); + + disp9.SetYTitle("Photon flux [ph/mm^2]"); + disp10.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp12.SetYTitle("Photon flux [ph/mm^2]"); + disp13.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp15.SetYTitle("[1]"); + disp16.SetYTitle("[1]"); + disp17.SetYTitle("[1]"); + disp18.SetYTitle("[1]"); + disp19.SetYTitle("[1]"); + disp20.SetYTitle("[1]"); + disp21.SetYTitle("[1]"); + disp22.SetYTitle("[1]"); + disp23.SetYTitle("[1]"); + disp24.SetYTitle("[1]"); + disp25.SetYTitle("[1]"); + + disp26.SetYTitle("Ped [FADC Counts ]"); + disp27.SetYTitle("RMS_{Ped} [FADC Counts ]"); + + disp28.SetYTitle("Time Offset [ns]"); + disp29.SetYTitle("Timing resolution [ns]"); + disp30.SetYTitle("P_{Time} [1]"); + + disp31.SetYTitle("[1]"); + disp32.SetYTitle("[1]"); + + disp33.SetYTitle("Mean Abs. Time [FADC slice]"); + disp34.SetYTitle("RMS Abs. Time [FADC slices]"); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + // Charges + TCanvas &c1 = display->AddTab("Fit.Charge"); + c1.Divide(2, 3); + + CamDraw(c1, disp1,calcam,1, 2 , 2); + CamDraw(c1, disp2,calcam,2, 2 , 2); + + // Fit Probability + TCanvas &c2 = display->AddTab("Fit.Prob"); + c2.Divide(1,3); + + CamDraw(c2, disp3,calcam,1, 1 , 4); + + // Reduced Sigmas + TCanvas &c3 = display->AddTab("Red.Sigma"); + c3.Divide(2,3); + + CamDraw(c3, disp4,calcam,1, 2 , 2); + CamDraw(c3, disp5,calcam,2, 2 , 2); + + // F-Factor Method + TCanvas &c4 = display->AddTab("F-Factor"); + c4.Divide(3,3); + + CamDraw(c4, disp6,calcam,1, 3 , 2); + CamDraw(c4, disp7,calcam,2, 3 , 2); + CamDraw(c4, disp8,calcam,3, 3 , 2); + + // Blind Pixel Method + TCanvas &c5 = display->AddTab("BlindPix"); + c5.Divide(3, 3); + + CamDraw(c5, disp9,calcam,1, 3 , 9); + CamDraw(c5, disp10,calcam,2, 3 , 2); + CamDraw(c5, disp11,calcam,3, 3 , 2); + + // PIN Diode Method + TCanvas &c6 = display->AddTab("PINDiode"); + c6.Divide(3,3); + + CamDraw(c6, disp12,calcam,1, 3 , 9); + CamDraw(c6, disp13,calcam,2, 3 , 2); + CamDraw(c6, disp14,calcam,3, 3 , 2); + + // Defects + TCanvas &c7 = display->AddTab("Defects"); + c7.Divide(4,2); + + CamDraw(c7, disp15,calcam,1,4, 0); + CamDraw(c7, disp16,calcam,2,4, 0); + CamDraw(c7, disp20,calcam,3,4, 0); + CamDraw(c7, disp21,calcam,4,4, 0); + + // BadCam + TCanvas &c8 = display->AddTab("Defects"); + c8.Divide(3,2); + + CamDraw(c8, disp17,badcam,1,3, 0); + CamDraw(c8, disp18,badcam,2,3, 0); + CamDraw(c8, disp19,badcam,3,3, 0); + + // Valid flags + TCanvas &c9 = display->AddTab("Validity"); + c9.Divide(4,2); + + CamDraw(c9, disp22,calcam,1,4,0); + CamDraw(c9, disp23,calcam,2,4,0); + CamDraw(c9, disp24,calcam,3,4,0); + CamDraw(c9, disp25,calcam,4,4,0); + + // Pedestals + TCanvas &c10 = display->AddTab("Pedestals"); + c10.Divide(2,3); + + CamDraw(c10,disp26,calcam,1,2,1); + CamDraw(c10,disp27,calcam,2,2,2); + + // Rel. Times + TCanvas &c11 = display->AddTab("Fitted Rel. Times"); + c11.Divide(3,3); + + CamDraw(c11,disp28,calcam,1,3,2); + CamDraw(c11,disp29,calcam,2,3,2); + CamDraw(c11,disp30,calcam,3,3,4); + + // Time Defects + TCanvas &c12 = display->AddTab("Time Def."); + c12.Divide(2,2); + + CamDraw(c12, disp31,calcam,1,2, 0); + CamDraw(c12, disp32,calcam,2,2, 0); + + // Abs. Times + TCanvas &c13 = display->AddTab("Abs. Times"); + c13.Divide(2,3); + + CamDraw(c13,disp33,calcam,1,2,2); + CamDraw(c13,disp34,calcam,2,2,2); +#endif + /************************************************************************/ + /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */ + /************************************************************************/ + MParList plist3; + MTaskList tlist3; + plist3.AddToList(&tlist3); + + // containers + MCerPhotEvt nphot; + MPedPhotCam nphotrms; + + plist3.AddToList(&geomcam); + plist3.AddToList(&pedcam); + plist3.AddToList(&calcam); + plist3.AddToList(&badcam); + plist3.AddToList(&sigcam); + plist3.AddToList(&nphot); + plist3.AddToList(&nphotrms); + + + MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault; + if(calflag==0) + calMode=MCalibrate::kNone; + + + //tasks + MReadMarsFile read3("Events", pedname); + read3.DisableAutoScheme(); + + MExtractSignal extsig; + extsig.SetRange(hifirst,hilast,lofirst,lolast); + MCalibrate photcalc(calMode); + MPedPhotCalc photrmscalc; + + tlist3.AddToList(&read3); + tlist3.AddToList(&geomapl); + tlist3.AddToList(&extsig); + tlist3.AddToList(&photcalc); + tlist3.AddToList(&photrmscalc); + + // Create and setup the eventloop + MEvtLoop evtloop3; + evtloop3.SetParList(&plist3); + if (!evtloop3.Eventloop()) + return; + + tlist3.PrintStatistics(); + + /************************************************************************/ + /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */ + /************************************************************************/ + + MParList plist4; + MTaskList tlist4; + plist4.AddToList(&tlist4); + + // containers + MHillas hillas; + MSrcPosCam source; + MRawRunHeader runhead; + + plist4.AddToList(&geomcam); + plist4.AddToList(&pedcam); + plist4.AddToList(&calcam); + plist4.AddToList(&badcam); + plist4.AddToList(&nphot); + plist4.AddToList(&nphotrms); + plist4.AddToList(&source); + plist4.AddToList(&hillas); + plist4.AddToList(&runhead); + + //tasks + MReadMarsFile read4("Events"); + + for(Int_t i=0;iSetBorderMode(0); + MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist"); + // obj1->AddNotify(evt); + + c.cd(i+j); + gPad->SetBorderMode(0); + obj1->Draw(); + ((MHCamera*)obj1)->SetPrettyPalette(); + + if (fit != 0) + { + c.cd(i+2*j); + gPad->SetBorderMode(0); + TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName()); + +// obj2->Sumw2(); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); + const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); + const Double_t integ = obj2->Integral("width")/2.5066283; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + if (rms == 0. || width == 0. ) + return; + + switch (fit) + { + case 1: + TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); + sgaus->SetBit(kCanDelete); + sgaus->SetParNames("Area","#mu","#sigma"); + sgaus->SetParameters(integ/rms,mean,rms); + sgaus->SetParLimits(0,0.,integ); + sgaus->SetParLimits(1,min,max); + sgaus->SetParLimits(2,0,width/1.5); + obj2->Fit("sgaus","QLR"); + obj2->GetFunction("sgaus")->SetLineColor(kYellow); + break; + + case 2: + TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); + dgaus->SetBit(kCanDelete); + dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); + dgaus->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + dgaus->SetParLimits(0,integ-1.5,integ+1.5); + dgaus->SetParLimits(1,min+(width/10.),mean); + dgaus->SetParLimits(2,0,width/2.); + // The right-sided Gauss + dgaus->SetParLimits(3,0,integ); + dgaus->SetParLimits(4,mean,max-(width/10.)); + dgaus->SetParLimits(5,0,width/2.); + obj2->Fit("dgaus","QLRM"); + obj2->GetFunction("dgaus")->SetLineColor(kYellow); + break; + + case 3: + TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); + tgaus->SetBit(kCanDelete); + tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + tgaus->SetParameters(integ,(min+mean)/2,width/4., + integ/width/3.,(max+mean)/2.,width/4., + integ/width/3.,mean,width/2.); + // The left-sided Gauss + tgaus->SetParLimits(0,integ-1.5,integ+1.5); + tgaus->SetParLimits(1,min+(width/10.),mean); + tgaus->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + tgaus->SetParLimits(3,0.,integ); + tgaus->SetParLimits(4,mean,max-(width/10.)); + tgaus->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + tgaus->SetParLimits(6,0.,integ); + tgaus->SetParLimits(7,min,max); + tgaus->SetParLimits(8,width/4.,width/1.5); + obj2->Fit("tgaus","QLRM"); + obj2->GetFunction("tgaus")->SetLineColor(kYellow); + break; + case 4: + obj2->Fit("pol0","Q"); + obj2->GetFunction("pol0")->SetLineColor(kYellow); + break; + case 9: + break; + default: + obj2->Fit("gaus","Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } + + TArrayI s0(3); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + + TArrayI s1(3); + s1[0] = 3; + s1[1] = 4; + s1[2] = 5; + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + // Just to get the right (maximum) binning + TH1D *half[4]; + half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); + half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); + half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); + half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); + + for (int i=0; i<4; i++) + { + half[i]->SetLineColor(kRed+i); + half[i]->SetDirectory(0); + half[i]->SetBit(kCanDelete); + half[i]->Draw("same"); + } + + gPad->Modified(); + gPad->Update(); + + } +} +#endif + +Bool_t readDatacards(TString& filename) +{ + ifstream ifun(filename.Data()); + if(!ifun) + { + cout << "File " << filename << " not found" << endl; + return kFALSE; + } + + TString word; + Int_t ifile = 0; + + while(ifun >> word) + { + // skip comments + if(word[0]=='/' && word[1]=='/') + { + while(ifun.get()!='\n'); // skip line + continue; + } + + // number of events + if(strcmp(word.Data(),"NEVENTS")==0) + ifun >> nmaxevents; + + // pedestal file name + if(strcmp(word.Data(),"PFILE")==0) + { + if(pedname.Length()) + cout << "readDataCards Warning: overriding pedestal file name" << endl; + ifun >> pedname; + } + + // calibration file name + if(strcmp(word.Data(),"CFILE")==0) + { + if(calname.Length()) + cout << "readDataCards Warning: overriding calibration file name" << endl; + ifun >> calname; + } + + // number of data files + if(strcmp(word.Data(),"NFILES")==0) + { + if(nfiles) + { + Int_t dummy; + cout << "readDataCards Warning: trying to set a new value of number of files. Skiping..." << endl; + ifun >> dummy; + } + else + { + ifun >> nfiles; + datafile = new TString[nfiles]; + } + } + + + // data file name + if(strcmp(word.Data(),"DFILE")==0) + { + TString newfile; + ifun >> newfile; + if(ifile> outname; + } + + // calibration flag + if(strcmp(word.Data(),"CALFLAG")==0) + ifun >> calflag; + + // cleaning level + if(strcmp(word.Data(),"CLEANLEVEL")==0) + { + ifun >> lcore; + ifun >> ltail; + } + } + + + // Dump read values + cout << "************************************************" << endl; + cout << "* Datacards read from file " << filename << endl; + cout << "************************************************" << endl; + cout << "Pedestal file: " << pedname << endl; + cout << "Calibration file: " << calname << endl; + cout << "Number of data files: " << nfiles << endl; + cout << "Data files: " << endl; + for(int i=0;i