Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2745) +++ trunk/MagicSoft/Mars/Changelog (revision 2746) @@ -4,4 +4,31 @@ -*-*- END OF LINE -*-*- + + 2003/12/23: Wolfgang Wittek + + * macros/ONAnalysis.C + - new macro for the MAGIC analysis, corresponding to CT1Analysis.C + + * macros/AnalyseCT1.C + - deleted because outdated + + * manalysis/MPadOn.[h,cc] + - new class for the MAGIC analysis, corresponding to MCT1PadSchweizer + + * manalysis/MPadOnOFF.[h,cc] + - minor changes (printouts) + + * manalysis/MPedestalWorkaround.[h,cc] + - put zenith angle into MMcEvt + + * manalysis/Makefile + AnalysisLinkDef.h + - added MPadON + + * mmc/MMcEvt.hxx + - add member function SetTelescopeTheta(), SetTelescopePhi() + + + 2003/12/22: Thomas Bretz @@ -451,4 +478,6 @@ 0.7 in which the array containing pixel ratios was not initialized. + + Index: trunk/MagicSoft/Mars/macros/ONAnalysis.C =================================================================== --- trunk/MagicSoft/Mars/macros/ONAnalysis.C (revision 2746) +++ trunk/MagicSoft/Mars/macros/ONAnalysis.C (revision 2746) @@ -0,0 +1,3452 @@ + +//#include "MagicEgyEst.C" + +void InitBinnings(MParList *plist) +{ + gLog << "InitBinnings" << endl; + + //-------------------------------------------- + MBinning *binse = new MBinning("BinningE"); + //binse->SetEdgesLog(30, 1.0e2, 1.0e5); + + //This is Daniel's binning in energy: + binse->SetEdgesLog(14, 296.296, 86497.6); + plist->AddToList(binse); + + //-------------------------------------------- + + MBinning *binssize = new MBinning("BinningSize"); + binssize->SetEdgesLog(50, 10, 1.0e5); + plist->AddToList(binssize); + + MBinning *binsdistc = new MBinning("BinningDist"); + binsdistc->SetEdges(50, 0, 1.4); + plist->AddToList(binsdistc); + + MBinning *binswidth = new MBinning("BinningWidth"); + binswidth->SetEdges(50, 0, 1.0); + plist->AddToList(binswidth); + + MBinning *binslength = new MBinning("BinningLength"); + binslength->SetEdges(50, 0, 1.0); + plist->AddToList(binslength); + + MBinning *binsalpha = new MBinning("BinningAlpha"); + binsalpha->SetEdges(100, -100, 100); + plist->AddToList(binsalpha); + + MBinning *binsasym = new MBinning("BinningAsym"); + binsasym->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsasym); + + MBinning *binsm3l = new MBinning("BinningM3Long"); + binsm3l->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsm3l); + + MBinning *binsm3t = new MBinning("BinningM3Trans"); + binsm3t->SetEdges(50, -1.5, 1.5); + plist->AddToList(binsm3t); + + + //..... + MBinning *binsb = new MBinning("BinningSigmabar"); + binsb->SetEdges( 100, 0.0, 5.0); + plist->AddToList(binsb); + + MBinning *binth = new MBinning("BinningTheta"); + // this is Daniel's binning in theta + //Double_t yedge[8] = + // {9.41, 16.22, 22.68, 28.64, 34.03, 38.84, 43.08, 44.99}; + // this is our binning + Double_t yedge[9] = + {0.0, 17.5, 23.5, 29.5, 35.5, 42., 50., 60., 70.}; + TArrayD yed; + yed.Set(9,yedge); + binth->SetEdges(yed); + plist->AddToList(binth); + + MBinning *bincosth = new MBinning("BinningCosTheta"); + Double_t zedge[9]; + for (Int_t i=0; i<9; i++) + { + zedge[8-i] = cos(yedge[i]/kRad2Deg); + } + TArrayD zed; + zed.Set(9,zedge); + bincosth->SetEdges(zed); + plist->AddToList(bincosth); + + MBinning *binsdiff = new MBinning("BinningDiffsigma2"); + binsdiff->SetEdges(100, -5.0, 20.0); + plist->AddToList(binsdiff); + + // robert ---------------------------------------------- + MBinning *binsalphaf = new MBinning("BinningAlphaFlux"); + binsalphaf->SetEdges(100, -100, 100); + plist->AddToList(binsalphaf); + + MBinning *binsdifftime = new MBinning("BinningTimeDiff"); + binsdifftime->SetEdges(50, 0, 10); + plist->AddToList(binsdifftime); + + MBinning *binstime = new MBinning("BinningTime"); + binstime->SetEdges(50, 44500, 61000); + plist->AddToList(binstime); +} + + +void DeleteBinnings(MParList *plist) +{ + gLog << "DeleteBinnings" << endl; + + TObject *bin; + + //-------------------------------------------- + bin = plist->FindObject("BinningE"); + if (bin) delete bin; + + //-------------------------------------------- + + bin = plist->FindObject("BinningSize"); + if (bin) delete bin; + + bin = plist->FindObject("BinningDist"); + if (bin) delete bin; + + bin = plist->FindObject("BinningWidth"); + if (bin) delete bin; + + bin = plist->FindObject("BinningLength"); + if (bin) delete bin; + + bin = plist->FindObject("BinningAlpha"); + if (bin) delete bin; + + bin = plist->FindObject("BinningAsym"); + if (bin) delete bin; + + bin = plist->FindObject("BinningM3Long"); + if (bin) delete bin; + + bin = plist->FindObject("BinningM3Trans"); + if (bin) delete bin; + + //..... + bin = plist->FindObject("BinningSigmabar"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTheta"); + if (bin) delete bin; + + bin = plist->FindObject("BinningCosTheta"); + if (bin) delete bin; + + bin = plist->FindObject("BinningDiffsigma2"); + if (bin) delete bin; + + + // robert ---------------------------------------------- + bin = plist->FindObject("BinningAlphaFlux"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTimeDiff"); + if (bin) delete bin; + + bin = plist->FindObject("BinningTime"); + if (bin) delete bin; +} + + +//************************************************************************ +void ONAnalysis() +{ + + gLog << "Entry ONAnalysis()" << endl; + + gLog.SetNoColors(); + + if (gRandom) + delete gRandom; + gRandom = new TRandom3(0); + + //----------------------------------------------- + const char *offfile = + "/.magic/data14a/crab_2003Nov/rootdata/2003_11_29/*OffCrab1*.root"; + + const char *onfile = + "/.magic/data14a/crab_2003Nov/rootdata/2003_11_29/*Crab-Nebula*.root"; + const char *onfile1 = + "/.magic/data14a/crab_2003Nov/rootdata/2003_11_30/*CrabNebula*.root"; + + const char *mcfile = + "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root"; + //----------------------------------------------- + + // path for input for Mars + TString inPath = "/.magic/data14a/wittek/ONAnalysis/"; + + // path for output from Mars + TString outPath = "/.magic/data14a/wittek/ONAnalysis/"; + + //----------------------------------------------- + + //TEnv env("macros/CT1env.rc"); + //Bool_t printEnv = kFALSE; + + //************************************************************************ + + // Job A_ON : read ON data + // - generate sigmabar vs. Theta plot; + // - write root file for ON data (ON1.root); + + Bool_t JobA_ON = kTRUE; + Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? + Bool_t WON1 = kFALSE; // write out root file ON1.root ? + + + // Job A_MC : read MC gamma data, + // - read sigmabar vs. Theta plot from ON data + // - do padding; + // - write root file for MC gammas (MC1.root); + + Bool_t JobA_MC = kFALSE; + Bool_t WMC1 = kFALSE; // write out root file MC1.root ? + + + // Job B_RF_UP : read ON1.root (or MC1.root) file + // - if RTrainRF = TRUE : read in training matrices for hadrons and gammas + // - if CTrainRF = TRUE : create training matrices for hadrons and gammas + // - if RTree = TRUE : read in trees, otherwise create trees + // - calculate hadroness for method of RANDOM FOREST + // - update the input files with the hadroness (==>ON2.root or MC2.root) + + Bool_t JobB_RF_UP = kFALSE; + Bool_t CTrainRF = kFALSE; // create matrices of training events + Bool_t RTrainRF = kFALSE; // read in matrices of training events + Bool_t RTree = kFALSE; // read in trees (otherwise grow them) + Bool_t WRF = kFALSE; // update input root file ? + + + + + // Job B_SC_UP : read ON2.root (or MC2.root) file + // - depending on WParSC : create (or read in) supercuts parameter values + // - calculate hadroness for the SUPERCUTS + // - update the input files with the hadroness (==>ON3.root or MC3.root) + + Bool_t JobB_SC_UP = kFALSE; + Bool_t CMatrix = kFALSE; // create training and test matrices + Bool_t RMatrix = kFALSE; // read training and test matrices from file + Bool_t WOptimize = kFALSE; // do optimization using the training sample + // and write supercuts parameter values + // onto the file parSCfile + Bool_t RTest = kFALSE; // test the supercuts using the test matrix + Bool_t WSC = kFALSE; // update input root file ? + + + + // Job C: + // - read ON3.root and MC3.root files + // which should have been updated to contain the hadronnesses + // for the method of + // RF + // SUPERCUTS and + // - produce Neyman-Pearson plots + + Bool_t JobC = kFALSE; + + + // Job D : + // - select g/h separation method XX + // - read ON3 (or MC3) root file + // - apply cuts in hadronness + // - make plots + + Bool_t JobD = kFALSE; + + + + + // Job E_XX : extended version of E_XX (including flux plots) + // - select g/h separation method XX + // - read MC root file +   // - calculate eff. collection area + // - optimize energy estimator + // - read ON root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts + + + Bool_t JobE_XX = kFALSE; + Bool_t CCollArea= kFALSE; // calculate eff. collection areas + Bool_t OEEst = kFALSE; // optimize energy estimator + Bool_t WEX = kFALSE; // update root file ? + Bool_t WRobert = kFALSE; // write out Robert's file ? + + + + //************************************************************************ + + + //--------------------------------------------------------------------- + // Job A_ON + //========= + // read ON data file + + // - produce the 2D-histogram "sigmabar versus Theta" + // (SigmaTheta_ON.root) for ON data + // (to be used for the padding of the MC gamma data) + + // - write a file of ON events (ON1.root) + // (after the standard cuts, before the g/h separation) + // (to be used together with the corresponding MC gamma file (MC1.root) + // for the optimization of the g/h separation) + + + if (JobA_ON) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job A_ON" << endl; + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobA_ON, WHistON, WON1 = " + << (JobA_ON ? "kTRUE" : "kFALSE") << ", " + << (WHistON ? "kTRUE" : "kFALSE") << ", " + << (WON1 ? "kTRUE" : "kFALSE") << endl; + + + // name of input root file + TString filenamein(onfile); + TString filenamein1(onfile1); + + // name of output root file + TString outNameImage = outPath; + outNameImage += "ON"; + outNameImage += "1.root"; + + //-------------------------------------------------- + // use for padding sigmabar vs. Theta from ON data + TString typeHist = "ON"; + gLog << "typeHist = " << typeHist << endl; + + // name of file to conatin the histograms for the padding + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta_"; + outNameSigTh += typeHist; + outNameSigTh += ".root"; + + + //----------------------------------------------------------- + MTaskList tliston; + MParList pliston; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + TFile test("/.magic/data14a/crab_2003Nov/rootdata/2003_11_29/20031128_03127_P_Crab-Nebula_E.root"); + test.ls(); + + //MCT1ReadPreProc read(filenamein); + MReadMarsFile read("Events", filenamein); + read.ls(); + read.AddFile(filenamein1); + read.DisableAutoScheme(); + + //--------------------- + MF selectCrabNebula + ("(MRawRunHeader.fRunNumber>3126 && MRawRunHeader.fRunNumber<3216)||(MRawRunHeader.fRunNumber>3283 && MRawRunHeader.fRunNumber<3335)||(MRawRunHeader.fRunNumber>3339 && MRawRunHeader.fRunNumber<3417)"); + selectCrabNebula.SetInverted(); + + MF selectOffCrab1 + ("(MRawRunHeader.fRunNumber>3215 && MRawRunHeader.fRunNumber<3275)||(MRawRunHeader.fRunNumber>3416 && MRawRunHeader.fRunNumber<3444)"); + selectOffCrab1.SetInverted(); + + MF selectMrk421 + ("(MRawRunHeader.fRunNumber>3443 && MRawRunHeader.fRunNumber<3490)"); + selectMrk421.SetInverted(); + + MContinue contCrabNebula(&selectCrabNebula); + MContinue contOffCrab1(&selectOffCrab1); + MContinue contMrk421(&selectMrk421); + + + //.............................. + MGeomApply apply; + MMcPedestalCopy pcopy; + MMcPedestalNSBAdd pnsb; + + MPedestalWorkaround waround; + + // a way to find out whether one is dealing with MC : + MFDataMember fMC("MRawRunHeader.fRunType", '>', 255.5); // MC + fMC.SetName("Select MC"); + MFDataMember fDa("MRawRunHeader.fRunType", '<', 255.5); // data + fDa.SetName("Select Data"); + + MCerPhotCalc ncalc; + ncalc.SetFilter(&fMC); + MCerPhotAnal2 nanal; + nanal.SetFilter(&fDa); + //.............................. + + + //MPointingCorrCalc pointcorr(sourceName, "MPointingCorrCalc", + // "MPointingCorrCalc"); + MBlindPixelCalc blind; + blind.SetUseBlindPixels(); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + MFillH fillblind("BlindPixels[MHBlindPixels]", "MBlindPixels"); + fillblind.SetName("HBlind"); + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetName("HSigmaTheta"); + + MImgCleanStd clean; + + + // calculation of image parameters --------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + MHillasCalc hcalc; + hcalc.SetNameHillas(fHilName); + hcalc.SetNameHillasExt(fHilNameExt); + hcalc.SetNameNewImgPar(fImgParName); + + MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc); + hsrccalc.SetInput(fHilName); + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + // -------------------------------------------------- + + MFSelStandard selstandard(fHilNameSrc); + selstandard.SetHillasName(fHilName); + selstandard.SetImgParName(fImgParName); + selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); + MContinue contstandard(&selstandard); + contstandard.SetName("SelStandard"); + + //MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + MWriteRootFile write(outNameImage, "RECREATE"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + //MF daniel( "(MRawRunHeader.fRunNumber<13167)||(MRawRunHeader.fRunNumber>13167)" ); + //MContinue contdaniel(&daniel); + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&source); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + tliston.AddToList(&contCrabNebula); + + //...................... + tliston.AddToList(&fMC); + tliston.AddToList(&fDa); + tliston.AddToList(&apply); + tliston.AddToList(&pcopy); + //tliston.AddToList(&waround); + + tliston.AddToList(&pnsb); + tliston.AddToList(&ncalc); + tliston.AddToList(&nanal); + //...................... + + //tliston.AddToList(&pointcorr); + tliston.AddToList(&blind); + tliston.AddToList(&contbasic); + + tliston.AddToList(&fillblind); + tliston.AddToList(&sigbarcalc); + tliston.AddToList(&fillsigtheta); + tliston.AddToList(&clean); + + tliston.AddToList(&hcalc); + tliston.AddToList(&hsrccalc); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contstandard); + if (WON1) + tliston.AddToList(&write); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + //evtloop.ReadEnv(env, "", printEnv); + evtloop.SetProgressBar(&bar); + //if (WON1) + // evtloop.Write(); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + + pliston.FindObject("SigmaTheta", "MHSigmaTheta")->DrawClone(); + + pliston.FindObject("BlindPixels", "MHBlindPixels")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + + + //------------------------------------------- + // Write histograms onto a file + if (WHistON) + { + MHSigmaTheta *sigtheta = + (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); + + MHBlindPixels *blindpixels = + (MHBlindPixels*)pliston.FindObject("BlindPixels"); + if (!sigtheta || !blindpixels) + { + gLog << "Object 'SigmaTheta' or 'BlindPixels' not found" << endl; + return; + } + TH2D *fHSigTh = sigtheta->GetSigmaTheta(); + TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); + TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); + + TH2D *fHBlindId = blindpixels->GetBlindId(); + TH2D *fHBlindN = blindpixels->GetBlindN(); + + + TFile outfile(outNameSigTh, "RECREATE"); + fHSigTh->Write(); + fHSigPixTh->Write(); + fHDifPixTh->Write(); + + fHBlindId->Write(); + fHBlindN->Write(); + + gLog << "" << endl; + gLog << "File " << outNameSigTh << " was written out" << endl; + } + + + DeleteBinnings(&pliston); + + gLog << "Macro ONAnalysis : End of Job A_ON" << endl; + gLog << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job A_MC + //========= + + // read MC gamma data + + // - to pad them + // (using the 2D-histogram "sigmabar versus Theta" + // (SigmaTheta_ON.root) of the ON data) + + // - to write a file of padded MC gamma events (MC1.root) + // (after the standard cuts, before the g/h separation) + // (to be used together with the corresponding hadron file + // for the optimization of the g/h separation) + + + if (JobA_MC) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job A_MC" << endl; + + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobA_MC, WMC1 = " + << (JobA_MC ? "kTRUE" : "kFALSE") << ", " + << (WMC1 ? "kTRUE" : "kFALSE") << endl; + + + // name of input root file + TString filenamein(mcfile); + + // name of output root file + TString outNameImage = outPath; + outNameImage += "MC"; + outNameImage += "1.root"; + + //------------------------------------------------ + // use for padding sigmabar vs. Theta from ON data + TString typeHist = "ON"; + gLog << "typeHist = " << typeHist << endl; + + // name of file containing the histograms for the padding + TString outNameSigTh = outPath; + outNameSigTh += "SigmaTheta_"; + outNameSigTh += typeHist; + outNameSigTh += ".root"; + + + //------------------------------------ + // Get the histograms "2D-ThetaSigmabar" + // and "3D-ThetaPixSigma" + // and "3D-ThetaPixDiff" + // and "2D-IdBlindPixels" + // and "2D-NBlindPixels" + + + gLog << "Reading in file " << outNameSigTh << endl; + + TFile *infile = new TFile(outNameSigTh); + infile->ls(); + + TH2D *fHSigmaTheta = + (TH2D*) gROOT->FindObject("2D-ThetaSigmabar"); + if (!fHSigmaTheta) + { + gLog << "Object '2D-ThetaSigmabar' not found on root file" << endl; + return; + } + gLog << "Object '2D-ThetaSigmabar' was read in" << endl; + + TH3D *fHSigmaPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); + if (!fHSigmaPixTheta) + { + gLog << "Object '3D-ThetaPixSigma' not found on root file" << endl; + return; + } + gLog << "Object '3D-ThetaPixSigma' was read in" << endl; + + TH3D *fHDiffPixTheta = + (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); + if (!fHDiffPixTheta) + { + gLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; + return; + } + gLog << "Object '3D-ThetaPixDiff' was read in" << endl; + + + TH2D *fHIdBlindPixels = + (TH2D*) gROOT->FindObject("2D-IdBlindPixels"); + if (!fHIdBlindPixels) + { + gLog << "Object '2D-IdBlindPixels' not found on root file" << endl; + return; + } + gLog << "Object '2D-IdBlindPixels' was read in" << endl; + + TH2D *fHNBlindPixels = + (TH2D*) gROOT->FindObject("2D-NBlindPixels"); + if (!fHNBlindPixels) + { + gLog << "Object '2D-NBlindPixels' not found on root file" << endl; + return; + } + gLog << "Object '2D-NBlindPixels' was read in" << endl; + + //------------------------------------ + + MTaskList tlist; + MParList plist; + + char *sourceName = "MSrcPosCam"; + MSrcPosCam source(sourceName); + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)plist->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + //MCT1ReadPreProc read(filenamein); + MReadMarsFile read("Events", filenamein); + read.DisableAutoScheme(); + + //.............................. + MGeomApply apply; + MMcPedestalCopy pcopy; + MMcPedestalNSBAdd pnsb; + + MPedestalWorkaround waround; + + // a way to find out whether one is dealing with MC : + MFDataMember fMC("MRawRunHeader.fRunType", '>', 255.5); // MC + fMC.SetName("Select MC"); + MFDataMember fDa("MRawRunHeader.fRunType", '<', 255.5); // data + fDa.SetName("Select Data"); + + MCerPhotCalc ncalc; + ncalc.SetFilter(&fMC); + MCerPhotAnal2 nanal; + nanal.SetFilter(&fDa); + //.............................. + + MBlindPixelCalc blindbeforepad; + blindbeforepad.SetUseBlindPixels(); + blindbeforepad.SetName("BlindBeforePadding"); + + MBlindPixelCalc blind; + blind.SetUseBlindPixels(); + blind.SetName("BlindAfterPadding"); + + MFSelBasic selbasic; + MContinue contbasic(&selbasic); + contbasic.SetName("SelBasic"); + + + // There are 2 options for Thomas Schweizer's padding + // fPadFlag = 1 get Sigmabar from fHSigmaTheta + // and Sigma from fHDiffPixTheta + // fPadFlag = 2 get Sigma from fHSigmaPixTheta + + MPadON padthomas("MPadON","Task for the padding"); + padthomas.SetHistograms(fHSigmaTheta, fHSigmaPixTheta, fHDiffPixTheta, + fHIdBlindPixels, fHNBlindPixels); + padthomas.SetPadFlag(1); + + MFillH fillblind("MCBlindPixels[MHBlindPixels]", "MBlindPixels"); + fillblind.SetName("HBlind"); + + + //........................................... + + MSigmabarCalc sigbarcalc; + + MFillH fillsigtheta ("MCSigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetName("HSigmaTheta"); + + MImgCleanStd clean; + + // calculation of image parameters --------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + MHillasCalc hcalc; + hcalc.SetNameHillas(fHilName); + hcalc.SetNameHillasExt(fHilNameExt); + hcalc.SetNameNewImgPar(fImgParName); + + MHillasSrcCalc hsrccalc(sourceName, fHilNameSrc); + hsrccalc.SetInput(fHilName); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + // -------------------------------------------------- + + MFSelStandard selstandard(fHilNameSrc); + selstandard.SetHillasName(fHilName); + selstandard.SetImgParName(fImgParName); + selstandard.SetCuts(92, 4, 60, 0.4, 1.05, 0.0, 0.0); + MContinue contstandard(&selstandard); + contstandard.SetName("SelStandard"); + + + //MWriteRootFile &write = *(new MWriteRootFile(outNameImage)); + MWriteRootFile write(outNameImage, "RECREATE"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + InitBinnings(&plist); + + plist.AddToList(&source); + + + //***************************** + // entries in MTaskList + + tlist.AddToList(&read); + + //...................... + tlist.AddToList(&fMC); + tlist.AddToList(&fDa); + tlist.AddToList(&apply); + tlist.AddToList(&pcopy); + //tlist.AddToList(&waround); + + tlist.AddToList(&pnsb); + tlist.AddToList(&ncalc); + tlist.AddToList(&nanal); + //...................... + + tlist.AddToList(&blindbeforepad); + tlist.AddToList(&padthomas); + tlist.AddToList(&blind); + + tlist.AddToList(&contbasic); + tlist.AddToList(&fillblind); + tlist.AddToList(&sigbarcalc); + tlist.AddToList(&fillsigtheta); + tlist.AddToList(&clean); + + tlist.AddToList(&hcalc); + tlist.AddToList(&hsrccalc); + + tlist.AddToList(&hfill1); + tlist.AddToList(&hfill2); + tlist.AddToList(&hfill3); + tlist.AddToList(&hfill4); + tlist.AddToList(&hfill5); + + tlist.AddToList(&contstandard); + if (WMC1) + tlist.AddToList(&write); + + //***************************** + + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + //evtloop.ReadEnv(env, "", printEnv); + evtloop.SetProgressBar(&bar); + //if (WMC1) + // evtloop.Write(); + + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tlist.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + plist.FindObject("MCSigmaTheta", "MHSigmaTheta")->DrawClone(); + plist.FindObject("MCBlindPixels", "MHBlindPixels")->DrawClone(); + + plist.FindObject("MHHillas")->DrawClone(); + plist.FindObject("MHHillasExt")->DrawClone(); + plist.FindObject("MHHillasSrc")->DrawClone(); + plist.FindObject("MHNewImagePar")->DrawClone(); + plist.FindObject("MHStarMap")->DrawClone(); + + + + DeleteBinnings(&plist); + + gLog << "Macro ONAnalysis : End of Job A_MC" + << endl; + gLog << "=========================================================" + << endl; + } + + + + //--------------------------------------------------------------------- + // Job B_RF_UP + //============ + + + // - create (or read in) the matrices of training events for gammas + // and hadrons + // - create (or read in) the trees + // - then read ON1.root (or MC1.root) file + // - calculate the hadroness for the method of RANDOM FOREST + // - update input root file with the hadroness + + + if (JobB_RF_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job B_RF_UP" << endl; + + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = " + << (JobB_RF_UP ? "kTRUE" : "kFALSE") << ", " + << (RTrainRF ? "kTRUE" : "kFALSE") << ", " + << (CTrainRF ? "kTRUE" : "kFALSE") << ", " + << (RTree ? "kTRUE" : "kFALSE") << ", " + << (WRF ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // parameters for the random forest + Int_t NumTrees = 100; + Int_t NumTry = 3; + Int_t NdSize = 1; + + + TString hadRFName = "HadRF"; + Float_t maxhadronness = 0.23; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + + TString extin = "1.root"; + TString extout = "2.root"; + + //-------------------------------------------- + // for the analysis using ON data only set typeMatrixHadrons = "ON" + // ON and OFF data = "OFF" + TString typeMatrixHadrons = "ON"; + gLog << "typeMatrixHadrons = " << typeMatrixHadrons << endl; + + + // file to be updated (ON, OFF or MC) + + //TString typeInput = "ON"; + TString typeInput = "OFF"; + //TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + // name of input root file + TString NameData = outPath; + NameData += typeInput; + TString inNameData(NameData); + inNameData += extin; + gLog << "inNameData = " << inNameData << endl; + + // name of output root file + TString outNameData(NameData); + outNameData += extout; + gLog << "outNameData = " << outNameData << endl; + + //-------------------------------------------- + // files to be read for generating + // - the matrices of training events + // - and the root files of training and test events + + + // "hadrons" : + TString filenameHad = outPath; + filenameHad += typeMatrixHadrons; + filenameHad += extin; + Int_t howManyHadronsTrain = 12000; + Int_t howManyHadronsTest = 12000; + gLog << "filenameHad = " << filenameHad << ", howManyHadronsTrain = " + << howManyHadronsTrain << ", howManyHadronsTest = " + << howManyHadronsTest << endl; + + + // "gammas" : + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += extin; + Int_t howManyGammasTrain = 12000; + Int_t howManyGammasTest = 12000; + gLog << "filenameMC = " << filenameMC << ", howManyGammasTrain = " + << howManyGammasTrain << ", howManyGammasTest = " + << howManyGammasTest << endl; + + //-------------------------------------------- + // files for the matrices of training events + + TString NameGammas = outPath; + NameGammas += "RFmatrix_gammas_Train_"; + NameGammas += "MC"; + NameGammas += extin; + + TString NameHadrons = outPath; + NameHadrons += "RFmatrix_hadrons_Train_"; + NameHadrons += typeMatrixHadrons; + NameHadrons += extin; + + + //-------------------------------------------- + // root files for the training events + + TString NameGammasTrain = outPath; + NameGammasTrain += "RF_gammas_Train_"; + NameGammasTrain += "MC"; + TString inNameGammasTrain(NameGammasTrain); + inNameGammasTrain += extin; + TString outNameGammasTrain(NameGammasTrain); + outNameGammasTrain += extout; + + + TString NameHadronsTrain = outPath; + NameHadronsTrain += "RF_hadrons_Train_"; + NameHadronsTrain += typeMatrixHadrons; + TString inNameHadronsTrain(NameHadronsTrain); + inNameHadronsTrain += extin; + TString outNameHadronsTrain(NameHadronsTrain); + outNameHadronsTrain += extout; + + + //-------------------------------------------- + // root files for the test events + + TString NameGammasTest = outPath; + NameGammasTest += "RF_gammas_Test_"; + NameGammasTest += "MC"; + TString inNameGammasTest(NameGammasTest); + inNameGammasTest += extin; + TString outNameGammasTest(NameGammasTest); + outNameGammasTest += extout; + + TString NameHadronsTest = outPath; + NameHadronsTest += "RF_hadrons_Test_"; + NameHadronsTest += typeMatrixHadrons; + TString inNameHadronsTest(NameHadronsTest); + inNameHadronsTest += extin; + TString outNameHadronsTest(NameHadronsTest); + outNameHadronsTest += extout; + + //-------------------------------------------------------------------- + + + MHMatrix matrixg("MatrixGammas"); + matrixg.EnableGraphicalOutput(); + + matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrixg.AddColumn("MSigmabar.fSigmabar"); + matrixg.AddColumn("log10(MHillas.fSize)"); + matrixg.AddColumn("MHillasSrc.fDist"); + matrixg.AddColumn("MHillas.fWidth"); + matrixg.AddColumn("MHillas.fLength"); + matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))"); + matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + matrixg.AddColumn("MNewImagePar.fConc"); + matrixg.AddColumn("MNewImagePar.fLeakage1"); + + MHMatrix matrixh("MatrixHadrons"); + matrixh.EnableGraphicalOutput(); + + matrixh.AddColumns(matrixg.GetColumns()); + + //-------------------------------------------- + // file of trees of the random forest + + TString outRF = outPath; + outRF += "RF.root"; + + + //************************************************************************* + // read in matrices of training events +if (RTrainRF) + { + const char* mtxName = "MatrixGammas"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Get matrix for (gammas)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << NameGammas << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'NameGammas' + // + TFile fileg(NameGammas); + + matrixg.Read(mtxName); + matrixg.Print("SizeCols"); + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Get matrix for (hadrons)" << endl; + gLog << "matrix name = " << mtxName << endl; + gLog << "name of root file = " << NameHadrons << endl; + gLog << "" << endl; + + + // read in the object with the name 'mtxName' from file 'NameHadrons' + // + TFile fileh(NameHadrons); + + matrixh.Read(mtxName); + matrixh.Print("SizeCols"); + } + + + //************************************************************************* + // create matrices of training events + // and root files of training and test events + +if (CTrainRF) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << " Create matrices of training events and root files of training and test events" + << endl; + gLog << " Gammas :" << endl; + gLog << "---------" << endl; + + MParList plistg; + MTaskList tlistg; + + MReadMarsFile readg("Events", filenameMC); + readg.DisableAutoScheme(); + + TString mgname("costhg"); + MBinning bing("Binning"+mgname); + bing.SetEdges(10, 0., 1.0); + + MH3 gref("cos(MMcEvt.fTelescopeTheta)"); + gref.SetName(mgname); + MH::SetBinning(&gref.GetHist(), &bing); + for (Int_t i=1; i<=gref.GetNbins(); i++) + gref.GetHist().SetBinContent(i, 1.0); + + MFEventSelector2 selectorg(gref); + selectorg.SetNumMax(howManyGammasTrain+howManyGammasTest); + selectorg.SetName("selectGammasTrainTest"); + selectorg.SetInverted(); + selectorg.SetUseOrigDistribution(kTRUE); + + MContinue contg(&selectorg); + contg.SetName("ContGammas"); + + Double_t probg = ( (Double_t) howManyGammasTrain ) + / ( (Double_t)(howManyGammasTrain+howManyGammasTest) ); + MFRandomSplit splitg(probg); + + MFillH fillmatg("MatrixGammas"); + fillmatg.SetFilter(&splitg); + fillmatg.SetName("fillGammas"); + + //----------------------- + // for writing the root files of training and test events + // for gammas + + MWriteRootFile writetraing(inNameGammasTrain, "RECREATE"); + writetraing.SetName("WriteGammasTrain"); + writetraing.SetFilter(&splitg); + + writetraing.AddContainer("MRawRunHeader", "RunHeaders"); + writetraing.AddContainer("MTime", "Events"); + writetraing.AddContainer("MMcEvt", "Events"); + writetraing.AddContainer("ThetaOrig", "Events"); + writetraing.AddContainer("MSrcPosCam", "Events"); + writetraing.AddContainer("MSigmabar", "Events"); + writetraing.AddContainer("MHillas", "Events"); + writetraing.AddContainer("MHillasExt", "Events"); + writetraing.AddContainer("MHillasSrc", "Events"); + writetraing.AddContainer("MNewImagePar", "Events"); + + MContinue contgtrain(&splitg); + contgtrain.SetName("ContGammaTrain"); + + MWriteRootFile writetestg(inNameGammasTest, "RECREATE"); + writetestg.SetName("WriteGammasTest"); + + writetestg.AddContainer("MRawRunHeader", "RunHeaders"); + writetestg.AddContainer("MTime", "Events"); + writetestg.AddContainer("MMcEvt", "Events"); + writetestg.AddContainer("ThetaOrig", "Events"); + writetestg.AddContainer("MSrcPosCam", "Events"); + writetestg.AddContainer("MSigmabar", "Events"); + writetestg.AddContainer("MHillas", "Events"); + writetestg.AddContainer("MHillasExt", "Events"); + writetestg.AddContainer("MHillasSrc", "Events"); + writetestg.AddContainer("MNewImagePar", "Events"); + + //----------------------- + + //***************************** fill gammas *** + // entries in MParList + + plistg.AddToList(&tlistg); + InitBinnings(&plistg); + + plistg.AddToList(&matrixg); + + //***************************** + // entries in MTaskList + + tlistg.AddToList(&readg); + tlistg.AddToList(&contg); + + tlistg.AddToList(&splitg); + tlistg.AddToList(&fillmatg); + tlistg.AddToList(&writetraing); + tlistg.AddToList(&contgtrain); + + tlistg.AddToList(&writetestg); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtloopg; + evtloopg.SetName("FillGammaMatrix"); + evtloopg.SetParList(&plistg); + //evtloopg.ReadEnv(env, "", printEnv); + evtloopg.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtloopg.Eventloop(maxevents)) + return; + + tlistg.PrintStatistics(0, kTRUE); + + matrixg.Print("SizeCols"); + Int_t generatedgTrain = matrixg.GetM().GetNrows(); + if ( fabs(generatedgTrain-howManyGammasTrain) > + 3.0*sqrt(howManyGammasTrain) ) + { + gLog << "Macro ONAnalysis.C : no.of generated gamma training events (" + << generatedgTrain << ") is incompatible with the no.of requested events (" + << howManyGammasTrain << ")" << endl; + } + + + Int_t generatedgTest = writetestg.GetNumExecutions(); + if ( fabs(generatedgTest-howManyGammasTest) > + 3.0*sqrt(howManyGammasTest) ) + { + gLog << "Macro ONAnalysis.C : no.of generated gamma test events (" + << generatedgTest << ") is incompatible with the no.of requested events (" + << howManyGammasTest << ")" << endl; + } + + //***************************** fill hadrons *** + gLog << "---------------------------------------------------------------" + << endl; + gLog << " Hadrons :" << endl; + gLog << "----------" << endl; + + MParList plisth; + MTaskList tlisth; + + MReadMarsFile readh("Events", filenameHad); + readh.DisableAutoScheme(); + + TString mhname("costhh"); + MBinning binh("Binning"+mhname); + binh.SetEdges(10, 0., 1.0); + + //MH3 href("cos(MMcEvt.fTelescopeTheta)"); + //href.SetName(mhname); + //MH::SetBinning(&href.GetHist(), &binh); + //for (Int_t i=1; i<=href.GetNbins(); i++) + // href.GetHist().SetBinContent(i, 1.0); + + //use the original distribution from the gammas + MH3 &href = *(selectorg.GetHistOrig()); + + MFEventSelector2 selectorh(href); + selectorh.SetNumMax(howManyHadronsTrain+howManyHadronsTest); + selectorh.SetName("selectHadronsTrainTest"); + selectorh.SetInverted(); + + MContinue conth(&selectorh); + conth.SetName("ContHadrons"); + + Double_t probh = ( (Double_t) howManyHadronsTrain ) + / ( (Double_t)(howManyHadronsTrain+howManyHadronsTest) ); + MFRandomSplit splith(probh); + + MFillH fillmath("MatrixHadrons"); + fillmath.SetFilter(&splith); + fillmath.SetName("fillHadrons"); + + //----------------------- + // for writing the root files of training and test events + // for hadrons + + MWriteRootFile writetrainh(inNameHadronsTrain, "RECREATE"); + writetrainh.SetName("WriteHadronsTrain"); + writetrainh.SetFilter(&splith); + + writetrainh.AddContainer("MRawRunHeader", "RunHeaders"); + writetrainh.AddContainer("MTime", "Events"); + writetrainh.AddContainer("MMcEvt", "Events"); + writetrainh.AddContainer("ThetaOrig", "Events"); + writetrainh.AddContainer("MSrcPosCam", "Events"); + writetrainh.AddContainer("MSigmabar", "Events"); + writetrainh.AddContainer("MHillas", "Events"); + writetrainh.AddContainer("MHillasExt", "Events"); + writetrainh.AddContainer("MHillasSrc", "Events"); + writetrainh.AddContainer("MNewImagePar", "Events"); + + MContinue conthtrain(&splith); + + MWriteRootFile writetesth(inNameHadronsTest, "RECREATE"); + writetesth.SetName("WriteHadronsTest"); + + writetesth.AddContainer("MRawRunHeader", "RunHeaders"); + writetesth.AddContainer("MTime", "Events"); + writetesth.AddContainer("MMcEvt", "Events"); + writetesth.AddContainer("ThetaOrig", "Events"); + writetesth.AddContainer("MSrcPosCam", "Events"); + writetesth.AddContainer("MSigmabar", "Events"); + writetesth.AddContainer("MHillas", "Events"); + writetesth.AddContainer("MHillasExt", "Events"); + writetesth.AddContainer("MHillasSrc", "Events"); + writetesth.AddContainer("MNewImagePar", "Events"); + + + //***************************** + // entries in MParList + + plisth.AddToList(&tlisth); + InitBinnings(&plisth); + + plisth.AddToList(&matrixh); + + //***************************** + // entries in MTaskList + + tlisth.AddToList(&readh); + tlisth.AddToList(&conth); + + tlisth.AddToList(&splith); + tlisth.AddToList(&fillmath); + tlisth.AddToList(&writetrainh); + tlisth.AddToList(&conthtrain); + + tlisth.AddToList(&writetesth); + + //***************************** + + MProgressBar matrixbar; + MEvtLoop evtlooph; + evtlooph.SetName("FillHadronMatrix"); + evtlooph.SetParList(&plisth); + //evtlooph.ReadEnv(env, "", printEnv); + evtlooph.SetProgressBar(&matrixbar); + + Int_t maxevents = -1; + if (!evtlooph.Eventloop(maxevents)) + return; + + tlisth.PrintStatistics(0, kTRUE); + + matrixh.Print("SizeCols"); + Int_t generatedhTrain = matrixh.GetM().GetNrows(); + if ( fabs(generatedhTrain-howManyHadronsTrain) > + 3.0*sqrt(howManyHadronsTrain) ) + { + gLog << "Macro ONAnalysis.C : no.of generated hadron training events (" + << generatedhTrain << ") is incompatible with the no.of requested events (" + << howManyHadronsTrain << ")" << endl; + } + + + Int_t generatedhTest = writetesth.GetNumExecutions(); + if ( fabs(generatedhTest-howManyHadronsTest) > + 3.0*sqrt(howManyHadronsTest) ) + { + gLog << "Macro ONAnalysis.C : no.of generated gamma test events (" + << generatedhTest << ") is incompatible with the no.of requested events (" + << howManyHadronsTest << ")" << endl; + } + + + //***************************************************** + + + // write out matrices of training events + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Write out matrices of training events" << endl; + + + //------------------------------------------- + // "gammas" + gLog << "Gammas :" << endl; + matrixg.Print("SizeCols"); + + TFile writeg(NameGammas, "RECREATE", ""); + matrixg.Write(); + + gLog << "" << endl; + gLog << "Macro ONAnalysis : matrix of training events for gammas written onto file " + << NameGammas << endl; + + //------------------------------------------- + // "hadrons" + gLog << "Hadrons :" << endl; + matrixh.Print("SizeCols"); + + TFile writeh(NameHadrons, "RECREATE", ""); + matrixh.Write(); + + gLog << "" << endl; + gLog << "Macro ONAnalysis : matrix of training events for hadrons written onto file " + << NameHadrons << endl; + + } + //********** end of creating matrices of training events *********** + + + MRanForest *fRanForest; + MRanTree *fRanTree; + //----------------------------------------------------------------- + // read in the trees of the random forest + if (RTree) + { + MParList plisttr; + MTaskList tlisttr; + plisttr.AddToList(&tlisttr); + + MReadTree readtr("TREE", outRF); + readtr.DisableAutoScheme(); + + MRanForestFill rffill; + rffill.SetNumTrees(NumTrees); + + // list of tasks for the loop over the trees + + tlisttr.AddToList(&readtr); + tlisttr.AddToList(&rffill); + + //------------------- + // Execute tree loop + // + MEvtLoop evtlooptr; + evtlooptr.SetName("ReadRFTrees"); + evtlooptr.SetParList(&plisttr); + if (!evtlooptr.Eventloop()) + return; + + tlisttr.PrintStatistics(0, kTRUE); + + gLog << "ONAnalysis : RF trees were read in from file " + << outRF << endl; + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plisttr->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plisttr->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + + //----------------------------------------------------------------- + // grow the trees of the random forest (event loop = tree loop) + + if (!RTree) + { + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Macro ONAnalysis : start growing trees" << endl; + + MTaskList tlist2; + MParList plist2; + plist2.AddToList(&tlist2); + + plist2.AddToList(&matrixg); + plist2.AddToList(&matrixh); + + MRanForestGrow rfgrow2; + rfgrow2.SetNumTrees(NumTrees); + rfgrow2.SetNumTry(NumTry); + rfgrow2.SetNdSize(NdSize); + + MWriteRootFile rfwrite2(outRF); + rfwrite2.AddContainer("MRanTree", "TREE"); + + MFillH fillh2("MHRanForestGini"); + + // list of tasks for the loop over the trees + + tlist2.AddToList(&rfgrow2); + tlist2.AddToList(&rfwrite2); + tlist2.AddToList(&fillh2); + + //------------------- + // Execute tree loop + // + MEvtLoop treeloop; + treeloop.SetName("GrowRFTrees"); + treeloop.SetParList(&plist2); + + if ( !treeloop.Eventloop() ) + return; + + tlist2.PrintStatistics(0, kTRUE); + + plist2.FindObject("MHRanForestGini")->DrawClone(); + + + // get adresses of objects which are used in the next eventloop + fRanForest = (MRanForest*)plist2->FindObject("MRanForest"); + if (!fRanForest) + { + gLog << err << dbginf << "MRanForest not found... aborting." << endl; + return kFALSE; + } + + fRanTree = (MRanTree*)plist2->FindObject("MRanTree"); + if (!fRanTree) + { + gLog << err << dbginf << "MRanTree not found... aborting." << endl; + return kFALSE; + } + + } + // end of growing the trees of the random forest + //----------------------------------------------------------------- + + + //----------------------------------------------------------------- + // Update the root files with the RF hadronness + // + + if (WRF) + { + //TString fileName(inNameHadronsTrain); + //TString outName(outNameHadronsTrain); + + //TString fileName(inNameHadronsTest); + //TString outName(outNameHadronsTest); + + //TString fileName(inNameGammasTrain); + //TString outName(outNameGammasTrain); + + //TString fileName(inNameGammasTest); + //TString outName(outNameGammasTest); + + TString fileName(inNameData); + TString outName(outNameData); + + + + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update root file '" << fileName + << "' with the RF hadronness; ==> " << outName << endl; + + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", fileName); + read.DisableAutoScheme(); + + + //....................................................................... + // calculate hadronnes for method of RANDOM FOREST + + + MRanForestCalc rfcalc; + rfcalc.SetHadronnessName(hadRFName); + + + //....................................................................... + + //MWriteRootFile write(outName, "UPDATE"); + MWriteRootFile write(outName, "RECREATE"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + write.AddContainer(hadRFName, "Events"); + + //----------------------------------------------------------------- + + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadRFName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillranfor("MHRanForest"); + fillranfor.SetName("HRanForest"); + + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadRFName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(fRanForest); + pliston.AddToList(fRanTree); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&rfcalc); + tliston.AddToList(&fillranfor); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetName("UpdateRootFile"); + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + pliston.FindObject("MHRanForest")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->Print(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); + TH1 &alphaHist = absalpha->GetHist(); + alphaHist.SetXTitle("|alpha| [\\circ]"); + alphaHist.SetName("alpha-macro"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << fileName << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + + DeleteBinnings(&pliston); + } + + gLog << "Macro ONAnalysis : End of Job B_RF_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + //--------------------------------------------------------------------- + // Job B_SC_UP + //============ + + // - create (or read in) optimum supercuts parameter values + // + // - calculate the hadroness for the supercuts + // + // - update input root file, including the hadroness + + + if (JobB_SC_UP) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job B_SC_UP" << endl; + + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobB_SC_UP, CMatrix, RMatrix, WOptimize, RTest, WSC = " + << (JobB_SC_UP ? "kTRUE" : "kFALSE") << ", " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " + << (RMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "kTRUE" : "kFALSE") << ", " + << (WSC ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // file which contains the initial parameter values for the supercuts + // if parSCinit ="" the initial values are taken from the constructor of + // MSupercuts + + TString parSCinit = outPath; + //parSCinit += "parSC_1709d"; + parSCinit = ""; + + gLog << "parSCinit = " << parSCinit << endl; + + //--------------- + // file onto which the optimal parameter values for the supercuts + // are written + + TString parSCfile = outPath; + parSCfile += "parSC_2310a"; + + gLog << "parSCfile = " << parSCfile << endl; + + //-------------------------------------------- + // file to be updated (either ON or MC) + + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; + gLog << "typeInput = " << typeInput << endl; + + // name of input root file + TString filenameData = outPath; + filenameData += typeInput; + filenameData += "2.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of output root file + TString outNameImage = outPath; + outNameImage += typeInput; + outNameImage += "3.root"; + + + //TString outNameImage = filenameData; + + gLog << "outNameImage = " << outNameImage << endl; + + //-------------------------------------------- + // files to be read for optimizing the supercuts + // + // for the training + TString filenameTrain = outPath; + filenameTrain += "ON"; + filenameTrain += "1.root"; + Int_t howManyTrain = 800000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = outPath; + filenameTest += "ON"; + filenameTest += "1.root"; + Int_t howManyTest = 800000; + + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // files to contain the matrices (generated from filenameTrain and + // filenameTest) + // + // for the training + TString fileMatrixTrain = outPath; + fileMatrixTrain += "MatrixTrainSC"; + fileMatrixTrain += ".root"; + gLog << "fileMatrixTrain = " << fileMatrixTrain << endl; + + // for testing + TString fileMatrixTest = outPath; + fileMatrixTest += "MatrixTestSC"; + fileMatrixTest += ".root"; + gLog << "fileMatrixTest = " << fileMatrixTest << endl; + + + + //--------------------------------------------------------------------- + // Training and test matrices : + // - either create them and write them onto a file + // - or read them from a file + + + MFindSupercuts findsuper; + findsuper.SetFilenameParam(parSCfile); + findsuper.SetHadronnessName("HadSC"); + findsuper.SetUseOrigDistribution(kTRUE); + + //-------------------------- + // create matrices and write them onto files + if (CMatrix) + { + TString mname("costheta"); + MBinning bin("Binning"+mname); + bin.SetEdges(10, 0., 1.0); + + MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mname); + MH::SetBinning(&mh3.GetHist(), &bin); + for (Int_t i=1; i<=mh3.GetNbins(); i++) + mh3.GetHist().SetBinContent(i, 1.0); + + + if (filenameTrain == filenameTest) + { + if ( !findsuper.DefineTrainTestMatrix( + filenameTrain, mh3, + howManyTrain, howManyTest, + fileMatrixTrain, fileMatrixTest) ) + { + *fLog << "Macro ONAnalysis.C : DefineTrainTestMatrix failed" << endl; + return; + } + + } + else + { + if ( !findsuper.DefineTrainMatrix(filenameTrain, mh3, + howManyTrain, fileMatrixTrain) ) + { + *fLog << "Macro ONAnalysis.C : DefineTrainMatrix failed" << endl; + return; + } + + if ( !findsuper.DefineTestMatrix( filenameTest, mh3, + howManyTest, fileMatrixTest) ) + { + *fLog << "Macro ONAnalysis.C : DefineTestMatrix failed" << endl; + return; + } + } + } + + //-------------------------- + // read matrices from files + // + + if (RMatrix) + findsuper.ReadMatrix(fileMatrixTrain, fileMatrixTest); + //-------------------------- + + + + //--------------------------------------------------------------------- + // optimize supercuts using the training sample + // + // the initial values are taken + // - from the file parSCinit (if != "") + // - or from the arrays params and steps (if their sizes are != 0) + // - or from the MSupercuts constructor + + +if (WOptimize) + { + gLog << "Macro ONAnalysis.C : optimize the supercuts using the training matrix" + << endl; + + TArrayD params(0); + TArrayD steps(0); + + if (parSCinit == "") + { + Double_t vparams[104] = { + // LengthUp + 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362, + 0.007388, -0.013463, + // LengthLo + 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993, + 0.000080, -0.007157, + // WidthUp + 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353, + 0.020711, -0.016703, + // WidthLo + 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101, + 0.006126, -0.002849, + // DistUp + 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909, + 0.189697, 0.0, + // DistLo + 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374, + -0.001750, 0.0, + // AsymUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + Double_t vsteps[104] = { + // LengthUp + 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002, + 0.0008, 0.002, + // LengthLo + 0.02, 0.003, 0.05, 0.006, 0.002, 0.3, + 0.0001, 0.0008, + // WidthUp + 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006, + 0.002, 0.002, + // WidthLo + 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002, + 0.0007, 0.003, + // DistUp + 0.2, 0.0, 0.3, 0.02, 0.0, 0.03, + 0.02, 0.0 + // DistLo + 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005, + 0.0002, 0.0 + // AsymUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + params.Set(104, vparams); + steps.Set (104, vsteps ); + } + + Bool_t rf; + rf = findsuper.FindParams(parSCinit, params, steps); + + if (!rf) + { + gLog << "Macro ONAnalysis.C : optimization of supercuts failed" << endl; + return; + } + } + + //-------------------------------------- + // test the supercuts on the test sample + // + + if (RTest) + { + gLog << "Macro ONAnalysis.C : test the supercuts on the test matrix" << endl; + Bool_t rt = findsuper.TestParams(); + if (!rt) + { + gLog << "Macro ONAnalysis.C : test of supercuts on the test matrix failed" + << endl; + } + + } + + + //----------------------------------------------------------------- + // Update the input files with the SC hadronness + // + + if (WSC) + { + gLog << "" << endl; + gLog << "========================================================" << endl; + gLog << "Update input file '" << filenameData + << "' with the SC hadronness" << endl; + + + //---------------------------------------------------- + // read in optimum parameter values for the supercuts + + TFile inparam(parSCfile); + MSupercuts scin; + scin.Read("MSupercuts"); + inparam.Close(); + + gLog << "Parameter values for supercuts were read in from file '" + << parSCfile << "'" << endl; + + TArrayD supercutsPar; + supercutsPar = scin.GetParameters(); + + TArrayD supercutsStep; + supercutsStep = scin.GetStepsizes(); + + gLog << "Parameter values for supercuts : " << endl; + for (Int_t i=0; iFindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + + //....................................................................... + // calculation of hadroness for the supercuts + // (=0.25 if fullfilled, =0.75 otherwise) + + TString hadSCName = "HadSC"; + MSupercutsCalc sccalc(fHilName, fHilNameSrc); + sccalc.SetHadronnessName(hadSCName); + + + //....................................................................... + + + //MWriteRootFile write(outNameImage, "UPDATE"); + //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE"); + + + MWriteRootFile write(outNameImage, "RECREATE"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + write.AddContainer("HadRF", "Events"); + write.AddContainer(hadSCName, "Events"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadSCName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadSCName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + + TH1 &alphahist = alphaabs->GetHist(); + + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + pliston.AddToList(&supercuts); + + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&sccalc); + tliston.AddToList(&fillhadsc); + + tliston.AddToList(&write); + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + //------------------------------------------- + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); + TH1 &alphaHist = absalpha->GetHist(); + alphaHist.SetXTitle("|alpha| [\\circ]"); + alphaHist.SetName("alpha-macro"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << filenameData << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + DeleteBinnings(&pliston); + } + + + gLog << "Macro ONAnalysis : End of Job B_SC_UP" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + //--------------------------------------------------------------------- + // Job C + //====== + + // - read ON1 and MC1 data files + // which should have been updated to contain the hadronnesses + // for the method of Random Forest and for the SUPERCUTS + // - produce Neyman-Pearson plots + + if (JobC) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job C" << endl; + + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobC = " + << (JobC ? "kTRUE" : "kFALSE") << endl; + + + // name of input data file + TString filenameData = outPath; + filenameData += "ON"; + filenameData += "3.root"; + gLog << "filenameData = " << filenameData << endl; + + // name of input MC file + TString filenameMC = outPath; + filenameMC += "MC"; + filenameMC += "3.root"; + gLog << "filenameMC = " << filenameMC << endl; + + + //----------------------------------------------------------------- + + MTaskList tliston; + MParList pliston; + + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + + //....................................................................... + // names of hadronness containers + + //TString hadNNName = "HadNN"; + TString hadSCName = "HadSC"; + TString hadRFName = "HadRF"; + + //....................................................................... + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + Float_t maxhadronness = 0.40; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(hadSCName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + //MFillH fillhadnn("hadNN[MHHadronness]", hadNNName); + //fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", hadSCName); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", hadRFName); + fillhadrf.SetName("HhadRF"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(hadSCName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + //tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 35000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + DeleteBinnings(&pliston); + + gLog << "Macro ONAnalysis : End of Job C" << endl; + gLog << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job D + //====== + + // - select g/h separation method XX + // - read ON2 (or MC2) root file + // - apply cuts in hadronness + // - make plots + + + if (JobD) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job D" << endl; + + gLog << "" << endl; + gLog << "Macro ONAnalysis : JobD = " + << (JobD ? "kTRUE" : "kFALSE") << endl; + + + // type of data to be analysed + TString typeData = "ON"; + //TString typeData = "OFF"; + //TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString ext = "3.root"; + + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.233; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + + + //************************************************************************* + // + // Analyse the data + // + + MTaskList tliston; + MParList pliston; + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + TString mh3name = "abs(Alpha)"; + MBinning binsalphaabs("Binning"+mh3name); + binsalphaabs.SetEdges(50, -2.0, 98.0); + + MH3 alphaabs("abs(MHillasSrc.fAlpha)"); + alphaabs.SetName(mh3name); + + TH1 &alphahist = alphaabs->GetHist(); + + MFillH alpha(&alphaabs); + alpha.SetName("FillAlphaAbs"); + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + pliston.AddToList(&binsalphaabs); + pliston.AddToList(&alphaabs); + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + tliston.AddToList(&contfinalgh); + + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&alpha); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + //Int_t maxevents = 10000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + pliston.FindObject("MHHillas")->DrawClone(); + pliston.FindObject("MHHillasExt")->DrawClone(); + pliston.FindObject("MHHillasSrc")->DrawClone(); + pliston.FindObject("MHNewImagePar")->DrawClone(); + pliston.FindObject("MHStarMap")->DrawClone(); + + + //------------------------------------------- + + // fit alpha distribution to get the number of excess events and + // calculate significance of gamma signal in the alpha plot + + MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3")); + TH1 &alphaHist = absalpha->GetHist(); + alphaHist.SetXTitle("|alpha| [\\circ]"); + alphaHist.SetName("alpha-JobD"); + + Double_t alphasig = 13.1; + Double_t alphamin = 30.0; + Double_t alphamax = 90.0; + Int_t degree = 2; + Double_t significance = -99.0; + Bool_t drawpoly = kTRUE; + Bool_t fitgauss = kTRUE; + Bool_t print = kTRUE; + + MHFindSignificance findsig; + findsig.SetRebin(kTRUE); + findsig.SetReduceDegree(kFALSE); + + findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, + alphasig, drawpoly, fitgauss, print); + significance = findsig.GetSignificance(); + Float_t alphasi = findsig.GetAlphasi(); + + gLog << "For file '" << filenameData << "' : " << endl; + gLog << "Significance of gamma signal after supercuts : " + << significance << " (for |alpha| < " << alphasi << " degrees)" + << endl; + + findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print); + + //------------------------------------------- + + + DeleteBinnings(&pliston); + + gLog << "Macro ONAnalysis : End of Job D" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + + + + + //--------------------------------------------------------------------- + // Job E_XX + //========= + + // - select g/h separation method XX + // - read MC_XX2.root file +   // - calculate eff. collection area + // - read ON_XX2.root file + // - apply final cuts + // - calculate flux + // - write root file for ON data after final cuts (ON_XX3.root)) + + + if (JobE_XX) + { + gLog << "=====================================================" << endl; + gLog << "Macro ONAnalysis : Start of Job E_XX" << endl; + + gLog << "" << endl; + gLog << "Macro +Analysis : JobE_XX, CCollArea, OEEst, WEX = " + << (JobE_XX ? "kTRUE" : "kFALSE") << ", " + << (CCollArea?"kTRUE" : "kFALSE") << ", " + << (OEEst ? "kTRUE" : "kFALSE") << ", " + << (WEX ? "kTRUE" : "kFALSE") << endl; + + + // type of data to be analysed + //TString typeData = "ON"; + //TString typeData = "OFF"; + TString typeData = "MC"; + gLog << "typeData = " << typeData << endl; + + TString typeMC = "MC"; + TString ext = "3.root"; + TString extout = "4.root"; + + //------------------------------ + // selection of g/h separation method + // and definition of final selections + + //TString XX("SC"); + TString XX("RF"); + TString fhadronnessName("Had"); + fhadronnessName += XX; + gLog << "fhadronnessName = " << fhadronnessName << endl; + + // maximum values of the hadronness, |ALPHA| and DIST + Float_t maxhadronness = 0.23; + Float_t maxalpha = 20.0; + Float_t maxdist = 10.0; + gLog << "Maximum values of hadronness, |ALPHA| and DIST = " + << maxhadronness << ", " << maxalpha << ", " + << maxdist << endl; + + //------------------------------ + // name of MC file to be used for optimizing the energy estimator + TString filenameOpt(outPath); + filenameOpt += typeMC; + filenameOpt += ext; + gLog << "filenameOpt = " << filenameOpt << endl; + + //------------------------------ + // name of file containing the parameters of the energy estimator + TString energyParName(outPath); + energyParName += "energyest_"; + energyParName += XX; + energyParName += ".root"; + gLog << "energyParName = " << energyParName << endl; + + //------------------------------ + // name of MC file to be used for calculating the eff. collection areas + TString filenameArea(outPath); + filenameArea += typeMC; + filenameArea += ext; + gLog << "filenameArea = " << filenameArea << endl; + + //------------------------------ + // name of file containing the eff. collection areas + TString collareaName(outPath); + collareaName += "area_"; + collareaName += XX; + collareaName += ".root"; + gLog << "collareaName = " << collareaName << endl; + + //------------------------------ + // name of data file to be analysed + TString filenameData(outPath); + filenameData += typeData; + filenameData += ext; + gLog << "filenameData = " << filenameData << endl; + + //------------------------------ + // name of output data file (after the final cuts) + TString filenameDataout(outPath); + filenameDataout += typeData; + filenameDataout += "_"; + filenameDataout += XX; + filenameDataout += extout; + gLog << "filenameDataout = " << filenameDataout << endl; + + //------------------------------ + // name of file containing histograms for flux calculastion + TString filenameResults(outPath); + filenameResults += typeData; + filenameResults += "Results_"; + filenameResults += XX; + filenameResults += extout; + gLog << "filenameResults = " << filenameResults << endl; + + + //==================================================================== + + MHMcCT1CollectionArea collarea; + collarea.SetEaxis(MHMcCT1CollectionArea::kLinear); + + MParList parlist; + InitBinnings(&parlist); + + if (CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Start calculation of effective collection areas" << endl; + + + MTaskList tasklist; + + //--------------------------------------- + // Setup the tasks to be executed + // + MReadMarsFile reader("Events", filenameArea); + reader.DisableAutoScheme(); + + MFSelFinal cuthadrons; + cuthadrons.SetHadronnessName(fhadronnessName); + cuthadrons.SetCuts(maxhadronness, maxalpha, maxdist); + + MContinue conthadrons(&cuthadrons); + + + MFillH filler("MHMcCT1CollectionArea", "MMcEvt"); + filler.SetName("CollectionArea"); + + //******************************** + // entries in MParList + + parlist.AddToList(&tasklist); + + parlist.AddToList(&collarea); + + //******************************** + // entries in MTaskList + + tasklist.AddToList(&reader); + tasklist.AddToList(&conthadrons); + tasklist.AddToList(&filler); + + //******************************** + + //----------------------------------------- + // Execute event loop + // + MEvtLoop magic; + magic.SetParList(&parlist); + + MProgressBar bar; + magic.SetProgressBar(&bar); + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(0, kTRUE); + + // Calculate effective collection areas + // and display the histograms + // + //MHMcCT1CollectionArea *collarea = + // (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea"); + collarea.CalcEfficiency(); + collarea.DrawClone(); + + + + //--------------------------------------------- + // Write histograms to a file + // + + TFile f(collareaName, "RECREATE"); + //collarea.GetHist()->Write(); + //collarea.GetHAll()->Write(); + //collarea.GetHSel()->Write(); + collarea.Write(); + + f.Close(); + + gLog << "Collection area plots written onto file " << collareaName << endl; + + gLog << "Calculation of effective collection areas done" << endl; + gLog << "-----------------------------------------------" << endl; + //------------------------------------------------------------------ + } + + if (!CCollArea) + { + gLog << "-----------------------------------------------" << endl; + gLog << "Read in effective collection areas from file " + << collareaName << endl; + + TFile collfile(collareaName); + collfile.ls(); + collarea.Read("MHMcCT1CollectionArea"); + collarea.DrawClone(); + + gLog << "Effective collection areas were read in from file " + << collareaName << endl; + gLog << "-----------------------------------------------" << endl; + } + + + // save binnings for call to CT1EEst + MBinning *binsE = (MBinning*)parlist.FindObject("BinningE"); + if (!binsE) + { + gLog << "Object 'BinningE' not found in MParList" << endl; + return; + } + MBinning *binsTheta = (MBinning*)parlist.FindObject("BinningTheta"); + if (!binsTheta) + { + gLog << "Object 'BinningTheta' not found in MParList" << endl; + return; + } + + //------------------------------------- + TString fHilName = "MHillas"; + TString fHilNameExt = "MHillasExt"; + TString fHilNameSrc = "MHillasSrc"; + TString fImgParName = "MNewImagePar"; + + + if (OEEst) + { + //=========================================================== + // + // Optimization of energy estimator + // + gLog << "Macro ONAnalysis.C : calling CT1EEst" << endl; + + TString inpath(""); + TString outpath(""); + Int_t howMany = 2000; + CT1EEst(inpath, filenameOpt, outpath, energyParName, + fHilName, fHilNameSrc, fhadronnessName, + howMany, maxhadronness, maxalpha, maxdist, + binsE, binsTheta); + gLog << "Macro ONAnalysis.C : returning from CT1EEst" << endl; + } + + if (WEX) + { + //----------------------------------------------------------- + // + // Read in parameters of energy estimator ("MMcEnergyEst") + // and migration matrix ("MHMcEnergyMigration") + // + gLog << "================================================================" + << endl; + gLog << "Macro ONAnalysis.C : read parameters of energy estimator and migration matrix from file '" + << energyParName << "'" << endl; + TFile enparam(energyParName); + enparam.ls(); + MMcEnergyEst mcest("MMcEnergyEst"); + mcest.Read("MMcEnergyEst"); + + //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst")); + gLog << "Parameters of energy estimator were read in" << endl; + + + gLog << "Read in Migration matrix" << endl; + + MHMcEnergyMigration mighiston("MHMcEnergyMigration"); + mighiston.Read("MHMcEnergyMigration"); + //MHMcEnergyMigration &mighiston = + // *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration")); + + gLog << "Migration matrix was read in" << endl; + + + TArrayD parA(mcest.GetNumCoeffA()); + TArrayD parB(mcest.GetNumCoeffB()); + for (Int_t i=0; iFindCreateObj("MGeomCamMagic", "MGeomCam"); + + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + //....................................................................... + + gLog << "Macro ONAnalysis.C : write root file '" << filenameDataout + << "'" << endl; + + //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout)); + + + MWriteRootFile write(filenameDataout, "RECREATE"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MSrcPosCam", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MHillasSrc", "Events"); + write.AddContainer("MNewImagePar", "Events"); + + //write.AddContainer("HadNN", "Events"); + write.AddContainer("HadSC", "Events"); + write.AddContainer("HadRF", "Events"); + + write.AddContainer("MEnergyEst", "Events"); + + + //----------------------------------------------------------------- + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + MFSelFinal selfinalgh(fHilNameSrc); + selfinalgh.SetCuts(maxhadronness, 100.0, maxdist); + selfinalgh.SetHadronnessName(fhadronnessName); + selfinalgh.SetName("SelFinalgh"); + MContinue contfinalgh(&selfinalgh); + contfinalgh.SetName("ContSelFinalgh"); + + //MFillH fillhadnn("hadNN[MHHadronness]", "HadNN"); + //fillhadnn.SetName("HhadNN"); + MFillH fillhadsc("hadSC[MHHadronness]", "HadSC"); + fillhadsc.SetName("HhadSC"); + MFillH fillhadrf("hadRF[MHHadronness]", "HadRF"); + fillhadrf.SetName("HhadRF"); + + //--------------------------- + // calculate estimated energy + + MEnergyEstParam eeston(fHilName); + eeston.Add(fHilNameSrc); + + eeston.SetCoeffA(parA); + eeston.SetCoeffB(parB); + + //--------------------------- + // calculate estimated energy using Daniel's parameters + + //MEnergyEstParamDanielMkn421 eeston(fHilName); + //eeston.Add(fHilNameSrc); + //eeston.SetCoeffA(parA); + //eeston.SetCoeffB(parB); + + + //--------------------------- + + + MFillH hfill1("MHHillas", fHilName); + hfill1.SetName("HHillas"); + + MFillH hfill2("MHStarMap", fHilName); + hfill2.SetName("HStarMap"); + + MFillH hfill3("MHHillasExt", fHilNameSrc); + hfill3.SetName("HHillasExt"); + + MFillH hfill4("MHHillasSrc", fHilNameSrc); + hfill4.SetName("HHillasSrc"); + + MFillH hfill5("MHNewImagePar", fImgParName); + hfill5.SetName("HNewImagePar"); + + //--------------------------- + // new from Robert + + MFillH hfill6("MHTimeDiffTheta", "MMcEvt"); + hfill6.SetName("HTimeDiffTheta"); + + MFillH hfill6a("MHTimeDiffTime", "MMcEvt"); + hfill6a.SetName("HTimeDiffTime"); + + MFillH hfill7("MHAlphaEnergyTheta", fHilNameSrc); + hfill7.SetName("HAlphaEnergyTheta"); + + MFillH hfill7a("MHAlphaEnergyTime", fHilNameSrc); + hfill7a.SetName("HAlphaEnergyTime"); + + MFillH hfill7b("MHThetabarTime", fHilNameSrc); + hfill7b.SetName("HThetabarTime"); + + MFillH hfill7c("MHEnergyTime", "MMcEvt"); + hfill7c.SetName("HEnergyTime"); + + + //--------------------------- + + MFSelFinal selfinal(fHilNameSrc); + selfinal.SetCuts(maxhadronness, maxalpha, maxdist); + selfinal.SetHadronnessName(fhadronnessName); + selfinal.SetName("SelFinal"); + MContinue contfinal(&selfinal); + contfinal.SetName("ContSelFinal"); + + + //***************************** + // entries in MParList + + pliston.AddToList(&tliston); + InitBinnings(&pliston); + + + //***************************** + // entries in MTaskList + + tliston.AddToList(&read); + + // robert + tliston.AddToList(&hfill6); //timediff + tliston.AddToList(&hfill6a); //timediff + + tliston.AddToList(&contfinalgh); + tliston.AddToList(&eeston); + + tliston.AddToList(&write); + + //tliston.AddToList(&fillhadnn); + tliston.AddToList(&fillhadsc); + tliston.AddToList(&fillhadrf); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + tliston.AddToList(&hfill5); + + //robert + tliston.AddToList(&hfill7); + tliston.AddToList(&hfill7a); + tliston.AddToList(&hfill7b); + tliston.AddToList(&hfill7c); + + tliston.AddToList(&contfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = -1; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + + //pliston.FindObject("hadNN", "MHHadronness")->DrawClone(); + + gLog << "before hadRF" << endl; + pliston.FindObject("hadRF", "MHHadronness")->DrawClone(); + + gLog << "before hadSC" << endl; + pliston.FindObject("hadSC", "MHHadronness")->DrawClone(); + + gLog << "before MHHillas" << endl; + pliston.FindObject("MHHillas")->DrawClone(); + + gLog << "before MHHillasExt" << endl; + pliston.FindObject("MHHillasExt")->DrawClone(); + + gLog << "before MHHillasSrc" << endl; + pliston.FindObject("MHHillasSrc")->DrawClone(); + + gLog << "before MHNewImagePar" << endl; + pliston.FindObject("MHNewImagePar")->DrawClone(); + + gLog << "before MHStarMap" << endl; + pliston.FindObject("MHStarMap")->DrawClone(); + + gLog << "before DeleteBinnings" << endl; + + DeleteBinnings(&pliston); + + gLog << "before Robert's code" << endl; + + +//rwagner write all relevant histograms onto a file + + if (WRobert) + { + gLog << "=======================================================" << endl; + gLog << "Write results onto file '" << filenameResults << "'" << endl; + + TFile outfile(filenameResults,"recreate"); + + MHHillasSrc* hillasSrc = + (MHHillasSrc*)(pliston->FindObject("MHHillasSrc")); + TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha()); + alphaHist->Write(); + gLog << "Alpha plot has been written out" << endl; + + + MHAlphaEnergyTheta* aetH = + (MHAlphaEnergyTheta*)(pliston->FindObject("MHAlphaEnergyTheta")); + TH3D* aetHist = (TH3D*)(aetH->GetHist()); + aetHist->SetName("aetHist"); + aetHist->Write(); + gLog << "AlphaEnergyTheta plot has been written out" << endl; + + MHAlphaEnergyTime* aetH2 = + (MHAlphaEnergyTime*)(pliston->FindObject("MHAlphaEnergyTime")); + TH3D* aetHist2 = (TH3D*)(aetH2->GetHist()); + aetHist2->SetName("aetimeHist"); +// aetHist2->DrawClone(); + aetHist2->Write(); + gLog << "AlphaEnergyTime plot has been written out" << endl; + + MHThetabarTime* tbt = + (MHThetabarTime*)(pliston->FindObject("MHThetabarTime")); + TProfile* tbtHist = (TProfile*)(tbt->GetHist()); + tbtHist->SetName("tbtHist"); + tbtHist->Write(); + gLog << "ThetabarTime plot has been written out" << endl; + + MHEnergyTime* ent = + (MHEnergyTime*)(pliston->FindObject("MHEnergyTime")); + TH2D* entHist = (TH2D*)(ent->GetHist()); + entHist->SetName("entHist"); + entHist->Write(); + gLog << "EnergyTime plot has been written out" << endl; + + MHTimeDiffTheta *time = (MHTimeDiffTheta*)pliston.FindObject("MHTimeDiffTheta"); + TH2D* timeHist = (TH2D*)(time->GetHist()); + timeHist->SetName("MHTimeDiffTheta"); + timeHist->SetTitle("Time diffs"); + timeHist->Write(); + gLog << "TimeDiffTheta plot has been written out" << endl; + + + MHTimeDiffTime *time2 = (MHTimeDiffTime*)pliston.FindObject("MHTimeDiffTime"); + TH2D* timeHist2 = (TH2D*)(time2->GetHist()); + timeHist2->SetName("MHTimeDiffTime"); + timeHist2->SetTitle("Time diffs"); + timeHist2->Write(); + gLog << "TimeDiffTime plot has been written out" << endl; + +//rwagner write also collareas to same file + collarea->GetHist()->Write(); + collarea->GetHAll()->Write(); + collarea->GetHSel()->Write(); + gLog << "Effective collection areas have been written out" << endl; + +//rwagner todo: write alpha_cut, type of g/h sep (RF, SC, NN), type of data +//rwagner (ON/OFF/MC), MJDmin, MJDmax to this file + + gLog << "before closing outfile" << endl; + + //outfile.Close(); + gLog << "Results were written onto file '" << filenameResults + << "'" << endl; + gLog << "=======================================================" << endl; + } + + } + + gLog << "Macro ONAnalysis : End of Job E_XX" << endl; + gLog << "=======================================================" << endl; + } + //--------------------------------------------------------------------- + +} + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2745) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2746) @@ -48,4 +48,5 @@ #pragma link C++ class MCT1PadSchweizer+; #pragma link C++ class MCT1PadONOFF+; +#pragma link C++ class MPadON+; #pragma link C++ class MPadONOFF+; Index: trunk/MagicSoft/Mars/manalysis/MPadON.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadON.cc (revision 2746) +++ trunk/MagicSoft/Mars/manalysis/MPadON.cc (revision 2746) @@ -0,0 +1,940 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek, 12/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPadON +// +// This task applies padding such that for a given pixel and for a given +// Theta bin the resulting distribution of the pedestal sigma is identical +// to the distributions given by fHSigmaPixTheta and fHDiffPixTheta. +// +// The number of photons, its error and the pedestal sigmas are altered. +// On average, the number of photons added is zero. +// +// The formulas used can be found in Thomas Schweizer's Thesis, +// Section 2.2.1 +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// THE FOLLOWING OPTION HAS NOT BEEN TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +// +// The padding has to be done before the image cleaning because the +// image cleaning depends on the pedestal sigmas. +// +// For random numbers gRandom is used. +// +///////////////////////////////////////////////////////////////////////////// +#include "MPadON.h" + +#include + +#include +#include +#include +#include +#include + +#include "MBinning.h" +#include "MSigmabar.h" +#include "MMcEvt.hxx" +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" +#include "MGeomCam.h" + +#include "MCerPhotPix.h" +#include "MCerPhotEvt.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MBlindPixels.h" + +ClassImp(MPadON); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPadON::MPadON(const char *name, const char *title) +{ + fName = name ? name : "MPadON"; + fTitle = title ? title : "Task for the padding (Schweizer)"; + + fHSigmaTheta = NULL; + fHSigmaPixTheta = NULL; + fHDiffPixTheta = NULL; + fHBlindPixIdTheta = NULL; + fHBlindPixNTheta = NULL; + + fHSigmaPedestal = NULL; + fHPhotons = NULL; + fHNSB = NULL; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MPadON::~MPadON() +{ + if (fHSigmaPedestal != NULL) delete fHSigmaPedestal; + if (fHPhotons != NULL) delete fHPhotons; + if (fHNSB != NULL) delete fHNSB; +} + +// -------------------------------------------------------------------------- +// +// Set the references to the histograms to be used in the padding +// +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// fHSigmaPixTheta 3D-hiostogram (Theta, pixel, sigma) +// fHDiffPixTheta 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) +// fHBlindPixIdTheta 2D-histogram (Theta, blind pixel Id) +// fHBlindPixNTheta 2D-histogram (Theta, no.of blind pixels ) +// +// +void MPadON::SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, + TH2D *hist2Pix, TH2D *hist2PixN) +{ + fHSigmaTheta = hist2; + fHSigmaPixTheta = hist3; + fHDiffPixTheta = hist3Diff; + fHBlindPixIdTheta = hist2Pix; + fHBlindPixNTheta = hist2PixN; + + fHSigmaTheta->SetDirectory(NULL); + fHSigmaPixTheta->SetDirectory(NULL); + fHDiffPixTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetDirectory(NULL); + + Print(); +} + +// -------------------------------------------------------------------------- +// +// Set the option for the padding +// +// There are 2 options for the padding : +// +// 1) fPadFlag = 1 : +// Generate first a Sigmabar using the 2D-histogram Sigmabar vs. Theta +// (fHSigmaTheta). Then generate a pedestal sigma for each pixel using +// the 3D-histogram Theta, pixel no., Sigma^2-Sigmabar^2 +// (fHDiffPixTheta). +// +// This is the preferred option as it takes into account the +// correlations between the Sigma of a pixel and Sigmabar. +// +// NOT YET TESTED : +// 2) fPadFlag = 2 : +// Generate a pedestal sigma for each pixel using the 3D-histogram +// Theta, pixel no., Sigma (fHSigmaPixTheta). +// +void MPadON::SetPadFlag(Int_t padflag) +{ + fPadFlag = padflag; + *fLog << inf << "MPadON::SetPadFlag(); choose option " << fPadFlag << endl; +} + +// -------------------------------------------------------------------------- +// +// Set the pointers and prepare the histograms +// +Int_t MPadON::PreProcess(MParList *pList) +{ + if ( !fHSigmaTheta || !fHSigmaPixTheta || !fHDiffPixTheta || + !fHBlindPixIdTheta || !fHBlindPixNTheta) + { + *fLog << err + << "MPadON : At least one of the histograms needed for the padding is not defined ... aborting." + << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MPadON : MMcEvt not found... aborting." + << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << err << "MPadON : MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << err + << "MPadON : MGeomCam not found (no geometry information available)... aborting." + << endl; + return kFALSE; + } + + fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); + if (!fEvt) + { + *fLog << err << "MPadON : MCerPhotEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSigmabar) + { + *fLog << err << "MPadON : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlinds) + { + *fLog << err << "MPadON : MBlindPixels not found... aborting." << endl; + return kFALSE; + } + + + //-------------------------------------------------------------------- + // histograms for checking the padding + // + // plot pedestal sigmas + fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", + 100, 0.0, 5.0, 100, 0.0, 5.0); + fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); + fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); + + // plot no.of photons (before vs. after padding) + fHPhotons = new TH2D("Photons","Photons: after vs.before padding", + 100, -10.0, 90.0, 100, -10, 90); + fHPhotons->SetXTitle("no.of photons before padding"); + fHPhotons->SetYTitle("no.of photons after padding"); + + // plot of added NSB + fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0); + fHNSB->SetXTitle("no.of added NSB photons"); + fHNSB->SetYTitle("no.of pixels"); + + + //-------------------------------------------------------------------- + + fIter = 20; + + memset(fErrors, 0, sizeof(fErrors)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Do the Padding +// idealy the events to be padded should have been generated without NSB +// +Int_t MPadON::Process() +{ + *fLog << all << "Entry MPadON::Process();" << endl; + + // this is the conversion factor from the number of photons + // to the number of photo electrons + // later to be taken from MCalibration + // fPEperPhoton = PW * LG * QE * 1D + Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + + //------------------------------------------------ + // add pixels to MCerPhotEvt which are not yet in; + // set their numnber of photons equal to zero + + UInt_t imaxnumpix = fCam->GetNumPixels(); + + for (UInt_t i=0; iGetNumPixels(); + for (UInt_t j=0; jAddPixel(i, 0.0, (*fPed)[i].GetPedestalRms()); + } + + + + //----------------------------------------- + Int_t rc=0; + + const UInt_t npix = fEvt->GetNumPixels(); + + //fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << all << "before padding : " << endl; + //fSigmabar->Print(""); + + + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + // to simulate the situation that before the padding the NSB and + // electronic noise are zero : set Sigma = 0 for all pixels + //for (UInt_t i=0; ioperator[](i); + // Int_t j = pix.GetPixId(); + + // MPedestalPix &ppix = fPed->operator[](j); + // ppix.SetMeanRms(0.0); + //} + //$$$$$$$$$$$$$$$$$$$$$$$$$$ + + //------------------------------------------- + // Calculate average sigma of the event + // + Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "before padding : " << endl; + fSigmabar->Print(""); + Double_t sigbarold = sigbarold_phot * fPEperPhoton; + + Double_t sigbarold2 = sigbarold*sigbarold; + + + // for MC data : expect sigmabar to be zero before the padding + if (sigbarold_phot > 0) + { + *fLog << err + << "MPadON::Process(); sigmabar of event to be padded is > 0 : " + << sigbarold_phot << ". Stop event loop " << endl; + // input data should have sigmabar = 0; stop event loop + + rc = 1; + fErrors[rc]++; + return kERROR; + } + + const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + *fLog << all << "theta = " << theta << endl; + + + //------------------------------------------- + // for the current theta, + // generate blind pixels according to the histograms + // fHBlindPixNTheta and fHBlindPixIDTheta + // + + + Int_t binPix = fHBlindPixNTheta->GetXaxis()->FindBin(theta); + + if ( binPix < 1 || binPix > fHBlindPixNTheta->GetNbinsX() ) + { + *fLog << warn + << "MPadON::Process(); binNumber out of range : theta, binPix = " + << theta << ", " << binPix << "; Skip event " << endl; + // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; + return kCONTINUE; + } + + // numBlind is the number of blind pixels in this event + TH1D *nblind; + UInt_t numBlind; + + nblind = fHBlindPixNTheta->ProjectionY("", binPix, binPix, ""); + if ( nblind->GetEntries() == 0.0 ) + { + *fLog << warn << "MPadON::Process(); projection for Theta bin " + << binPix << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete nblind; + + rc = 7; + fErrors[rc]++; + return kCONTINUE; + } + else + { + numBlind = (Int_t) (nblind->GetRandom()+0.5); + } + delete nblind; + + + // throw the Id of numBlind different pixels in this event + TH1D *hblind; + UInt_t idBlind; + UInt_t listId[npix]; + UInt_t nlist = 0; + Bool_t equal; + + hblind = fHBlindPixIdTheta->ProjectionY("", binPix, binPix, ""); + if ( hblind->GetEntries() > 0.0 ) + for (UInt_t i=0; iGetRandom()+0.5); + equal = kFALSE; + for (UInt_t j=0; jSetPixelBlind(idBlind); + *fLog << "idBlind = " << idBlind << endl; + } + fBlinds->SetReadyToSave(); + + delete hblind; + + + //------------------------------------------- + // for the current theta, + // generate a sigmabar according to the histogram fHSigmaTheta + // + Double_t sigmabar_phot = 0; + Double_t sigmabar = 0; + Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(theta); + + if (binPix != binNumber) + { + *fLog << err + << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + TH1D *hsigma; + + hsigma = fHSigmaTheta->ProjectionY("", binNumber, binNumber, ""); + if ( hsigma->GetEntries() == 0.0 ) + { + *fLog << warn << "MPadON::Process(); projection for Theta bin " + << binNumber << " has no entries; Skip event " << endl; + // event cannot be padded; skip event + delete hsigma; + + rc = 3; + fErrors[rc]++; + return kCONTINUE; + } + else + { + sigmabar_phot = hsigma->GetRandom(); + sigmabar = sigmabar_phot * fPEperPhoton; + *fLog << all << "Theta, bin number = " << theta << ", " << binNumber + << ", sigmabar_phot = " << sigmabar_phot << endl; + } + delete hsigma; + + const Double_t sigmabar2 = sigmabar*sigmabar; + + //------------------------------------------- + + *fLog << all << "MPadON::Process(); sigbarold, sigmabar = " + << sigbarold << ", "<< sigmabar << endl; + + // Skip event if target sigmabar is <= sigbarold + if (sigmabar <= sigbarold) + { + *fLog << all + << "MPadON::Process(); target sigmabar is less than sigbarold : " + << sigmabar << ", " << sigbarold << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; + return kCONTINUE; + } + + + //------------------------------------------- + // + // Calculate average number of NSB photons to be added (lambdabar) + // from the value of sigmabar, + // - making assumptions about the average electronic noise (elNoise2) and + // - using a fixed value (F2excess) for the excess noise factor + + Double_t elNoise2; // [photo electrons] + Double_t F2excess = 1.3; + Double_t lambdabar; // [photo electrons] + + + + Int_t binTheta = fHDiffPixTheta->GetXaxis()->FindBin(theta); + if (binTheta != binNumber) + { + *fLog << err + << "MPadON::Process(); The binnings of the 2 histograms are not identical; aborting" + << endl; + return kERROR; + } + + // Get RMS of (Sigma^2-sigmabar^2) in this Theta bin. + // The average electronic noise (to be added) has to be well above this RMS, + // otherwise the electronic noise of an individual pixel (elNoise2Pix) + // may become negative + + TH1D *hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + 0, 9999, ""); + Double_t RMS_phot = hnoise->GetRMS(1); + Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton; + delete hnoise; + + elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); + *fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; + + // This value of lambdabar is the same for all pixels; + // note that lambdabar is normalized to the area of pixel 0 + + //---------- start loop over pixels --------------------------------- + // do the padding for each pixel + // + // pad only pixels - which are used (before image cleaning) + // + Double_t sig_phot = 0; + Double_t sig = 0; + + Double_t sigma2 = 0; + + Double_t diff_phot = 0; + Double_t diff = 0; + + Double_t addSig2_phot= 0; + Double_t addSig2 = 0; + + Double_t elNoise2Pix = 0; + + + for (UInt_t i=0; iGetPixRatio(j); + + MPedestalPix &ppix = (*fPed)[j]; + Double_t oldsigma_phot = ppix.GetPedestalRms(); + Double_t oldsigma = oldsigma_phot * fPEperPhoton; + Double_t oldsigma2 = oldsigma*oldsigma; + + //--------------------------------- + // throw the Sigma for this pixel + // + Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j ); + + Int_t count; + Bool_t ok; + + TH1D *hdiff; + TH1D *hsig; + + switch (fPadFlag) + { + case 1 : + // throw the Sigma for this pixel from the distribution fHDiffPixTheta + + hdiff = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + if ( hdiff->GetEntries() == 0 ) + { + *fLog << warn << "MPadON::Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hdiff; + + rc = 5; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + diff = diff_phot * fPEperPhoton * fPEperPhoton; + + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2/ratioArea + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + *fLog << all << "theta, j, count, sigmabar, diff = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << diff << endl; + diff = lambdabar*F2excess + oldsigma2/ratioArea - sigmabar2; + } + delete hdiff; + sigma2 = diff + sigmabar2; + break; + + case 2 : + // throw the Sigma for this pixel from the distribution fHSigmaPixTheta + + hsig = fHSigmaPixTheta->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + if ( hsig->GetEntries() == 0 ) + { + *fLog << err << "MPadON::Process(); projection for Theta bin " + << binTheta << " and pixel bin " << binPixel + << " has no entries; aborting " << endl; + delete hsig; + + rc = 6; + fErrors[rc]++; + return kCONTINUE; + } + + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + sig = sig_phot * fPEperPhoton; + + sigma2 = sig*sig/ratioArea; + // the following condition ensures that elNoise2Pix > 0.0 + + if ( (sigma2-oldsigma2/ratioArea-lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + if (!ok) + { + + *fLog << "theta, j, count, sigmabar, sig = " << theta << ", " + << j << ", " << count << ", " << sigmabar << ", " + << sig << endl; + sigma2 = lambdabar*F2excess + oldsigma2/ratioArea; + } + delete hsig; + break; + } + + //--------------------------------- + // get the additional sigma^2 for this pixel (due to the padding) + + addSig2 = sigma2*ratioArea - oldsigma2; + addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); + + //--------------------------------- + // get the additional electronic noise for this pixel + + elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea; + + + //--------------------------------- + // throw actual number of additional NSB photons (NSB) + // and its RMS (sigmaNSB) + + Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea); + Double_t arg = NSB0*(F2excess-1.0) + elNoise2Pix; + Double_t sigmaNSB0; + + if (arg >= 0) + { + sigmaNSB0 = sqrt( arg ); + } + else + { + sigmaNSB0 = 0.0000001; + if (arg < -1.e-10) + { + *fLog << warn << "MPadON::Process(); argument of sqrt < 0 : " + << arg << endl; + } + } + + + //--------------------------------- + // smear NSB0 according to sigmaNSB0 + // and subtract lambdabar because of AC coupling + + Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea; + Double_t NSB_phot = NSB / fPEperPhoton; + + //--------------------------------- + + // add additional NSB to the number of photons + Double_t oldphotons_phot = pix.GetNumPhotons(); + Double_t oldphotons = oldphotons_phot * fPEperPhoton; + Double_t newphotons = oldphotons + NSB; + Double_t newphotons_phot = newphotons / fPEperPhoton; + pix.SetNumPhotons(newphotons_phot); + + + fHNSB->Fill( NSB_phot/sqrt(ratioArea) ); + fHPhotons->Fill( oldphotons_phot/ratioArea, + newphotons_phot/ratioArea ); + + + // error: add sigma of padded noise quadratically + Double_t olderror_phot = pix.GetErrorPhot(); + Double_t newerror_phot = + sqrt(olderror_phot*olderror_phot + addSig2_phot); + pix.SetErrorPhot(newerror_phot); + + + Double_t newsigma = sqrt(oldsigma2 + addSig2); + Double_t newsigma_phot = newsigma / fPEperPhoton; + ppix.SetPedestalRms(newsigma_phot); + + fHSigmaPedestal->Fill(oldsigma_phot, newsigma_phot); + } + //---------- end of loop over pixels --------------------------------- + + // Calculate sigmabar again and crosscheck + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + *fLog << all << "MPadON::Process(); after padding : " << endl; + fSigmabar->Print(""); + + + *fLog << all << "Exit MPadON::Process();" << endl; + + rc = 0; + fErrors[rc]++; + + return kTRUE; + +} + +// -------------------------------------------------------------------------- +// +// +Int_t MPadON::PostProcess() +{ + if (GetNumExecutions() != 0) + { + + *fLog << inf << endl; + *fLog << GetDescriptor() << " execution statistics:" << endl; + *fLog << dec << setfill(' '); + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) + << (int)(fErrors[2]*100/GetNumExecutions()) + << "%) Evts skipped due to: Zenith angle out of range" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) + << (int)(fErrors[3]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigmabar" << endl; + + *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) + << (int)(fErrors[4]*100/GetNumExecutions()) + << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl; + + *fLog << " " << setw(7) << fErrors[5] << " (" << setw(3) + << (int)(fErrors[5]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma^2-Sigmabar^2" << endl; + + *fLog << " " << setw(7) << fErrors[6] << " (" << setw(3) + << (int)(fErrors[6]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Sigma" << endl; + + *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) + << (int)(fErrors[7]*100/GetNumExecutions()) + << "%) Evts skipped due to: No data for generating Blind pixels" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts survived the padding!" << endl; + *fLog << endl; + + } + + //--------------------------------------------------------------- + TCanvas &c = *(MH::MakeDefCanvas("PadON", "", 900, 1200)); + c.Divide(3, 4); + gROOT->SetSelectedPad(NULL); + + c.cd(1); + fHNSB->SetDirectory(NULL); + fHNSB->DrawCopy(); + fHNSB->SetBit(kCanDelete); + + c.cd(2); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->DrawCopy(); + fHSigmaPedestal->SetBit(kCanDelete); + + c.cd(3); + fHPhotons->SetDirectory(NULL); + fHPhotons->DrawCopy(); + fHPhotons->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(4); + fHSigmaTheta->SetDirectory(NULL); + fHSigmaTheta->SetTitle("(Input) 2D : Sigmabar, \\Theta"); + fHSigmaTheta->DrawCopy(); + fHSigmaTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(7); + fHBlindPixNTheta->SetDirectory(NULL); + fHBlindPixNTheta->SetTitle("(Input) 2D : no.of blind pixels, \\Theta"); + fHBlindPixNTheta->DrawCopy(); + fHBlindPixNTheta->SetBit(kCanDelete); + + //-------------------------------------------------------------------- + + + c.cd(10); + fHBlindPixIdTheta->SetDirectory(NULL); + fHBlindPixIdTheta->SetTitle("(Input) 2D : blind pixel Id, \\Theta"); + fHBlindPixIdTheta->DrawCopy(); + fHBlindPixIdTheta->SetBit(kCanDelete); + + + //-------------------------------------------------------------------- + // draw the 3D histogram (input): Theta, pixel, Sigma^2-Sigmabar^2 + + c.cd(5); + TH2D *l1; + l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx"); + l1->SetDirectory(NULL); + l1->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); + l1->SetXTitle("\\Theta [\\circ]"); + l1->SetYTitle("Sigma^2-Sigmabar^2"); + + l1->DrawCopy("box"); + l1->SetBit(kCanDelete);; + + c.cd(8); + TH2D *l2; + l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy"); + l2->SetDirectory(NULL); + l2->SetTitle("(Input) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); + l2->SetXTitle("pixel"); + l2->SetYTitle("Sigma^2-Sigmabar^2"); + + l2->DrawCopy("box"); + l2->SetBit(kCanDelete);; + + //-------------------------------------------------------------------- + // draw the 3D histogram (input): Theta, pixel, Sigma + + c.cd(6); + TH2D *k1; + k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx"); + k1->SetDirectory(NULL); + k1->SetTitle("(Input) Sigma vs. \\Theta (all pixels)"); + k1->SetXTitle("\\Theta [\\circ]"); + k1->SetYTitle("Sigma"); + + k1->DrawCopy("box"); + k1->SetBit(kCanDelete); + + c.cd(9); + TH2D *k2; + k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy"); + k2->SetDirectory(NULL); + k2->SetTitle("(Input) Sigma vs. pixel number (all \\Theta)"); + k2->SetXTitle("pixel"); + k2->SetYTitle("Sigma"); + + k2->DrawCopy("box"); + k2->SetBit(kCanDelete);; + + + //-------------------------------------------------------------------- + + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadON.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadON.h (revision 2746) +++ trunk/MagicSoft/Mars/manalysis/MPadON.h (revision 2746) @@ -0,0 +1,82 @@ +#ifndef MARS_MPadON +#define MARS_MPadON + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class TH3D; + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MSigmabar; +class MParList; +class MBlindPixels; + +class MPadON : public MTask +{ +private: + MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + MBlindPixels *fBlinds; + + Int_t fPadFlag; + Int_t fIter; + + Int_t fRunType; + Int_t fGroup; + + Int_t fErrors[8]; + + // plots used for the padding + TH2D *fHBlindPixIdTheta; // 2D-histogram (blind pixel Id vs. Theta) + TH2D *fHBlindPixNTheta; // 2D-histogram (no.of blind pixels vs. Theta) + TH2D *fHSigmaTheta; // 2D-histogram (sigmabar vs. Theta) + TH3D *fHSigmaPixTheta; // 3D-histogram (Theta, pixel, sigma) + TH3D *fHDiffPixTheta; // 3D-histogram (Theta, pixel, sigma^2-sigmabar^2) + + // plots for checking the padding + TH2D *fHSigmaPedestal; // 2D-histogram : pedestal sigma after + // versus before padding + TH2D *fHPhotons; // 2D-histogram : no.of photons after + // versus before padding + TH1D *fHNSB; // 1D-histogram : additional NSB + + +public: + MPadON(const char *name=NULL, const char *title=NULL); + ~MPadON(); + + void SetHistograms(TH2D *hist2, TH3D *hist3, TH3D *hist3Diff, + TH2D *hist2Pix, TH2D *hist2PixN); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + void SetPadFlag(Int_t padflag); + + ClassDef(MPadON, 0) // task for the padding (Schweizer) +}; + +#endif + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2745) +++ trunk/MagicSoft/Mars/manalysis/MPadONOFF.cc (revision 2746) @@ -110,5 +110,5 @@ fPadFlag = 1; - *fLog << all << "fPadFlag = " << fPadFlag << endl; + *fLog << all << "MadONOFF : fPadFlag = " << fPadFlag << endl; fType = ""; @@ -175,6 +175,6 @@ TH2D *blindnthon, TH2D *blindnthoff) { - *fLog << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "Merge the ON and OFF histograms to obtain the target distributions for the padding" + *fLog << all << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "MPadONOFF::MergeHistograms(); Merge the ON and OFF histograms to obtain the target distributions for the padding" << endl; @@ -242,5 +242,5 @@ // at the end, it should be zero - *fLog << all << "MergeHistograms; bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + *fLog << all << "MPadONOFF::MergeHistograms(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; for (Int_t l=1; l<=nbinstheta; l++) { @@ -417,5 +417,5 @@ if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) - *fLog << err << "Read; inconsistency in results; a, b, c = " + *fLog << err << "MPadONOFF::MergeHistograms(); inconsistency in results; a, b, c = " << a << ", " << b << ", " << c << endl; } @@ -642,5 +642,6 @@ - *fLog << all << "The target distributions for the padding have been created" + *fLog << all + << "MPadONOFF::MergeHistograms(); The target distributions for the padding have been created" << endl; *fLog << all << "----------------------------------------------------------" @@ -679,5 +680,6 @@ Bool_t MPadONOFF::ReadTargetDist(const char* namefilein) { - *fLog << all << "Read padding histograms from file " << namefilein << endl; + *fLog << all << "MPadONOFF : Read padding histograms from file " + << namefilein << endl; fInfile = new TFile(namefilein); @@ -690,9 +692,11 @@ if (!fHSigmaTheta) { - *fLog << all << "Object '2D-ThetaSigmabar' not found on root file" + *fLog << all + << "MPadONOFF : Object '2D-ThetaSigmabar' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '2D-ThetaSigmabar' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '2D-ThetaSigmabar' was read in" << endl; fHSigmaThetaON = @@ -700,9 +704,11 @@ if (!fHSigmaThetaON) { - *fLog << all << "Object '2D-ThetaSigmabarON' not found on root file" + *fLog << all + << "MPadONOFF : Object '2D-ThetaSigmabarON' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '2D-ThetaSigmabarON' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '2D-ThetaSigmabarON' was read in" << endl; fHSigmaThetaOFF = @@ -710,9 +716,11 @@ if (!fHSigmaThetaOFF) { - *fLog << "Object '2D-ThetaSigmabarOFF' not found on root file" + *fLog << all + << "MPadONOFF : Object '2D-ThetaSigmabarOFF' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '2D-ThetaSigmabarOFF' was read in" << endl; + *fLog << all + << "MPadONOFF :Object '2D-ThetaSigmabarOFF' was read in" << endl; fHSigmaPixTheta = @@ -720,9 +728,11 @@ if (!fHSigmaPixTheta) { - *fLog << all << "Object '3D-ThetaPixSigma' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigma' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-ThetaPixSigma' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigma' was read in" << endl; fHSigmaPixThetaON = @@ -730,9 +740,11 @@ if (!fHSigmaPixThetaON) { - *fLog << all << "Object '3D-ThetaPixSigmaON' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigmaON' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-ThetaPixSigmaON' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigmaON' was read in" << endl; fHSigmaPixThetaOFF = @@ -740,9 +752,11 @@ if (!fHSigmaPixThetaOFF) { - *fLog << all << "Object '3D-ThetaPixSigmaOFF' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' not found on root file" << endl; return kFALSE; } - *fLog << "Object '3D-ThetaPixSigmaOFF' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixSigmaOFF' was read in" << endl; fHDiffPixTheta = @@ -750,9 +764,11 @@ if (!fHDiffPixTheta) { - *fLog << all << "Object '3D-ThetaPixDiff' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiff' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-ThetaPixDiff' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiff' was read in" << endl; fHDiffPixThetaON = @@ -760,9 +776,11 @@ if (!fHDiffPixThetaON) { - *fLog << all << "Object '3D-ThetaPixDiffON' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiffON' not found on root file" << endl; return kFALSE; } - *fLog << "Object '3D-ThetaPixDiffON' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiffON' was read in" << endl; fHDiffPixThetaOFF = @@ -770,9 +788,11 @@ if (!fHDiffPixThetaOFF) { - *fLog << all << "Object '3D-ThetaPixDiffOFF' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiffOFF' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-ThetaPixDiffOFF' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-ThetaPixDiffOFF' was read in" << endl; fHgON = @@ -780,9 +800,11 @@ if (!fHgON) { - *fLog << all << "Object '3D-PaddingMatrixON' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-PaddingMatrixON' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-PaddingMatrixON' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-PaddingMatrixON' was read in" << endl; fHgOFF = @@ -790,9 +812,11 @@ if (!fHgOFF) { - *fLog << all << "Object '3D-PaddingMatrixOFF' not found on root file" + *fLog << all + << "MPadONOFF : Object '3D-PaddingMatrixOFF' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '3D-PaddingMatrixOFF' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '3D-PaddingMatrixOFF' was read in" << endl; @@ -801,9 +825,11 @@ if (!fHBlindPixIdTheta) { - *fLog << all << "Object '2D-ThetaBlindId' not found on root file" + *fLog << all + << "MPadONOFF : Object '2D-ThetaBlindId' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '2D-ThetaBlindId' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '2D-ThetaBlindId' was read in" << endl; fHBlindPixNTheta = @@ -811,9 +837,11 @@ if (!fHBlindPixNTheta) { - *fLog << all << "Object '2D-ThetaBlindN' not found on root file" + *fLog << all + << "MPadONOFF : Object '2D-ThetaBlindN' not found on root file" << endl; return kFALSE; } - *fLog << all << "Object '2D-ThetaBlindN' was read in" << endl; + *fLog << all + << "MPadONOFF : Object '2D-ThetaBlindN' was read in" << endl; @@ -831,5 +859,6 @@ Bool_t MPadONOFF::WriteTargetDist(const char* namefileout) { - *fLog << all << "Write padding histograms onto file " << namefileout << endl; + *fLog << all << "MPadONOFF : Write padding histograms onto file " + << namefileout << endl; TFile outfile(namefileout, "RECREATE"); @@ -853,5 +882,6 @@ fHgOFF->Write(); - *fLog << all << "WriteTargetDist; target histograms written onto file " + *fLog << all + << "MPadONOFF::WriteTargetDist(); target histograms written onto file " << namefileout << endl; @@ -870,5 +900,5 @@ { fType = type; - *fLog << all << "SetDataType(); type of data to be padded : " + *fLog << all << "MPadONOFF::SetDataType(); type of data to be padded : " << fType << endl; @@ -900,5 +930,6 @@ { fPadFlag = padflag; - *fLog << all << "SetPadFlag(); choose option " << fPadFlag << endl; + *fLog << all << "MPadONOFF::SetPadFlag(); choose option " << fPadFlag + << endl; } @@ -915,5 +946,7 @@ !fHgON || !fHgOFF) { - *fLog << err << "At least one of the histograms needed for the padding is not defined ... aborting." << endl; + *fLog << err + << "MPadONOFF : At least one of the histograms needed for the padding is not defined ... aborting." + << endl; return kFALSE; } @@ -922,5 +955,6 @@ if (!fMcEvt) { - *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; + *fLog << err << dbginf << "MPadONOFF : MMcEvt not found... aborting." + << endl; return kFALSE; } @@ -929,5 +963,6 @@ if (!fPed) { - *fLog << err << "MPedestalCam not found... aborting." << endl; + *fLog << err << "MPadONOFF : MPedestalCam not found... aborting." + << endl; return kFALSE; } @@ -936,5 +971,7 @@ if (!fCam) { - *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl; + *fLog << err + << "MPadONOFF : MGeomCam not found (no geometry information available)... aborting." + << endl; return kFALSE; } @@ -943,5 +980,6 @@ if (!fEvt) { - *fLog << err << "MCerPhotEvt not found... aborting." << endl; + *fLog << err << "MPadONOFF : MCerPhotEvt not found... aborting." + << endl; return kFALSE; } @@ -950,5 +988,5 @@ if (!fSigmabar) { - *fLog << err << "MSigmabar not found... aborting." << endl; + *fLog << err << "MPadONOFF : MSigmabar not found... aborting." << endl; return kFALSE; } @@ -957,5 +995,6 @@ if (!fBlinds) { - *fLog << err << "MBlindPixels not found... aborting." << endl; + *fLog << err << "MPadONOFF : MBlindPixels not found... aborting." + << endl; return kFALSE; } @@ -963,5 +1002,7 @@ if (fType !="ON" && fType !="OFF" && fType !="MC") { - *fLog << err << "Type of data to be padded not defined... aborting." << endl; + *fLog << err + << "MPadONOFF : Type of data to be padded not defined... aborting." + << endl; return kFALSE; } @@ -1005,10 +1046,10 @@ Int_t MPadONOFF::Process() { - *fLog << all << "Entry Process();" << endl; + *fLog << all << "Entry MPadONOFF::Process();" << endl; // this is the conversion factor from the number of photons // to the number of photo electrons // later to be taken from MCalibration - // fPEpePhoton = PW * LG * QE * 1D + // fPEperPhoton = PW * LG * QE * 1D Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; @@ -1069,5 +1110,5 @@ // Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "before padding : " << endl; + *fLog << all << "MPadONOFF::Process(); before padding : " << endl; fSigmabar->Print(""); Double_t sigbarold = sigbarold_phot * fPEperPhoton; @@ -1082,5 +1123,6 @@ if (sigbarold_phot > 0) { - *fLog << warn << "Process(); sigmabar of event to be padded is > 0 : " + *fLog << err + << "MPadONOFF::Process(); sigmabar of event to be padded is > 0 : " << sigbarold_phot << ". Stop event loop " << endl; // input data should have sigmabar = 0; stop event loop @@ -1088,5 +1130,5 @@ rc = 1; fErrors[rc]++; - return kCONTINUE; + return kERROR; } } @@ -1098,5 +1140,6 @@ if ( binTheta < 1 || binTheta > fHBlindPixNTheta->GetNbinsX() ) { - *fLog << warn << "Process(); binNumber out of range : theta, binTheta = " + *fLog << warn + << "MPadONOFF::Process(); binNumber out of range : theta, binTheta = " << theta << ", " << binTheta << "; Skip event " << endl; // event cannot be padded; skip event @@ -1170,5 +1213,5 @@ if ( nblind->Integral() == 0.0 ) { - *fLog << warn << "Process(); projection for Theta bin " + *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " << binTheta << " has no entries; Skip event " << endl; // event cannot be padded; skip event @@ -1319,7 +1362,7 @@ hsigma = fHSigmaTheta->ProjectionY("", binTheta, binTheta, ""); - if ( hsigma->Integral() == 0.0 ) + if ( hsigma->GetEntries() == 0.0 ) { - *fLog << warn << "Process(); projection for Theta bin " + *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " << binTheta << " has no entries; Skip event " << endl; // event cannot be padded; skip event @@ -1345,5 +1388,5 @@ //------------------------------------------- - *fLog << all << "Process(); sigbarold, sigmabar = " + *fLog << all << "MPadONOFF::Process(); sigbarold, sigmabar = " << sigbarold << ", "<< sigmabar << endl; @@ -1351,5 +1394,5 @@ if (sigmabar <= sigbarold) { - *fLog << all << "Process(); target sigmabar is less than sigbarold : " + *fLog << all << "MPadONOFF::Process(); target sigmabar is less than sigbarold : " << sigmabar << ", " << sigbarold << ", Skip event" << endl; @@ -1367,7 +1410,7 @@ // - using a fixed value (F2excess) for the excess noise factor - Double_t elNoise2; // [photons] + Double_t elNoise2; // [photo electrons] Double_t F2excess = 1.3; - Double_t lambdabar; // [photons] + Double_t lambdabar; // [photo electrons] @@ -1376,6 +1419,7 @@ if (binTheta != bincheck) { - cout << "The binnings of the 2 histograms are not identical; aborting" - << endl; + *fLog << err + << "MPadONOFF::Process(); The binnings of the 2 histograms are not identical; aborting" + << endl; return kERROR; } @@ -1395,5 +1439,6 @@ else { - *fLog << err << "Process; illegal data type... aborting" << endl; + *fLog << err << "MPadONOFF::Process(); illegal data type... aborting" + << endl; return kERROR; } @@ -1406,5 +1451,5 @@ *fLog << all << "elNoise2 = " << elNoise2 << endl; - lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // [photo electrons] + lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; // This value of lambdabar is the same for all pixels; @@ -1438,5 +1483,6 @@ //if ( pix.GetNumPhotons() == 0.0) //{ - // *fLog << warn << "Process(); no.of photons is 0 for used pixel" + // *fLog << warn + // << "MPadONOFF::Process(); no.of photons is 0 for used pixel" // << endl; // continue; @@ -1479,7 +1525,7 @@ binPixel, binPixel, ""); - if ( hdiff->Integral() == 0 ) + if ( hdiff->GetEntries() == 0 ) { - *fLog << warn << "Process(); projection for Theta bin " + *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " << binTheta << " and pixel bin " << binPixel << " has no entries; aborting " << endl; @@ -1541,7 +1587,7 @@ binPixel, binPixel, ""); - if ( hsig->Integral() == 0 ) + if ( hsig->GetEntries() == 0 ) { - *fLog << warn << "Process(); projection for Theta bin " + *fLog << warn << "MPadONOFF::Process(); projection for Theta bin " << binTheta << " and pixel bin " << binPixel << " has no entries; aborting " << endl; @@ -1621,5 +1667,5 @@ if (arg < -1.e-10) { - *fLog << err << "Process(); argument of sqrt < 0 : " + *fLog << warn << "MPadONOFF::Process(); argument of sqrt < 0 : " << arg << endl; } @@ -1640,5 +1686,5 @@ Double_t oldphotons = oldphotons_phot * fPEperPhoton; Double_t newphotons = oldphotons + NSB; - Double_t newphotons_phot = newphotons / fPEperPhoton;; + Double_t newphotons_phot = newphotons / fPEperPhoton; pix.SetNumPhotons( newphotons_phot ); @@ -1657,5 +1703,5 @@ Double_t newsigma = sqrt( oldsigma2 + addSig2 ); - Double_t newsigma_phot = newsigma / fPEperPhoton;; + Double_t newsigma_phot = newsigma / fPEperPhoton; ppix.SetPedestalRms( newsigma_phot ); @@ -1667,6 +1713,8 @@ fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "after padding : " << endl; + *fLog << all << "MPadONOFF::Process(); after padding : " << endl; fSigmabar->Print(""); + + *fLog << all << "Exit MPadONOFF::Process();" << endl; rc = 0; Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2745) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.cc (revision 2746) @@ -16,5 +16,6 @@ ! ! -! Author(s): Sabrina Stark 12/2003 +! Author(s): Sabrina Stark 12/2003 +! Wolfgang Wittek 12/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 @@ -27,6 +28,17 @@ // MPedestalWorkaround // // // -// This class contains the function to store the pedestal value and the // -// RMS of the pedestal in units of photons. // +// copies the pedestal value and the pedestal RMS // +// from the container MPedPhotCam // +// into the contaier MPedestalCam // +// // +// put the zenith angle into MMcEvt; // +// the zenith angle is taken from the runbooks // +// (1 fixed zenith angle for a given run) // +// // +// this workaround is necessary // +// - as long as the analysis classes // +// take the pedestals from MPedestalCam; in the long run they have to // +// be taken from MPedPhotCam // +// - and as long as the container for the zenith angle is not defined // // // ///////////////////////////////////////////////////////////////////////////// @@ -42,4 +54,6 @@ #include "MPedestalPix.h" #include "MPedPhotCam.h" +#include "MMcEvt.hxx" +#include "MRawRunHeader.h" using namespace std; @@ -65,4 +79,5 @@ return kFALSE; } + fPedPhot = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); if (!fPedPhot) @@ -71,4 +86,5 @@ return kFALSE; } + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fCam) @@ -78,5 +94,19 @@ } - return kTRUE; + fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRun) + { + *fLog << err << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; } @@ -85,6 +115,9 @@ Int_t MPedestalWorkaround::Process() { - - UInt_t imaxnumpix = fCam->GetNumPixels(); + //------------------------------------------------------------------- + // copy the pedestal and the pedestal RMS + // from MPedPhotCam into MPedestalCam + + UInt_t imaxnumpix = fCam->GetNumPixels(); for (UInt_t i=0; iGetRunNumber(); + Double_t thetadeg; + Double_t thetarad; + + if (run == 3127) thetadeg = 27.2; // Crab + else if (run == 3128) thetadeg = 25.6; + else if (run == 3129) thetadeg = 24.3; + else if (run == 3130) thetadeg = 23.9; + + + else if (run == 3132) thetadeg = 22.9; + else if (run == 3133) thetadeg = 22.8; + else if (run == 3134) thetadeg = 22.3; + else if (run == 3135) thetadeg = 21.9; + else if (run == 3136) thetadeg = 21.5; + else if (run == 3137) thetadeg = 21.1; + else if (run == 3138) thetadeg = 20.8; + else if (run == 3139) thetadeg = 20.4; + + else if (run == 3140) thetadeg = 19.5; + else if (run == 3141) thetadeg = 19.4; + else if (run == 3142) thetadeg = 19.0; + else if (run == 3143) thetadeg = 18.6; + else if (run == 3144) thetadeg = 13.0; + else if (run == 3145) thetadeg = 12.4; + else if (run == 3146) thetadeg = 12.1; + else if (run == 3147) thetadeg = 11.7; + else if (run == 3148) thetadeg = 11.3; + else if (run == 3149) thetadeg = 11.9; + + else if (run == 3150) thetadeg = 10.6; + else if (run == 3151) thetadeg = 10.3; + else if (run == 3152) thetadeg = 10.0; + else if (run == 3153) thetadeg = 9.6; + else if (run == 3154) thetadeg = 9.3; + else if (run == 3155) thetadeg = 9.0; + else if (run == 3156) thetadeg = 8.7; + else if (run == 3157) thetadeg = 8.4; + else if (run == 3158) thetadeg = 8.1; + else if (run == 3159) thetadeg = 7.9; + + else if (run == 3160) thetadeg = 7.7; + else if (run == 3161) thetadeg = 7.3; + else if (run == 3162) thetadeg = 7.2; + else if (run == 3163) thetadeg = 7.0; + else if (run == 3164) thetadeg = 6.8; + else if (run == 3165) thetadeg = 6.7; + else if (run == 3166) thetadeg = 6.6; + else if (run == 3167) thetadeg = 6.5; + else if (run == 3168) thetadeg = 6.4; + else if (run == 3169) thetadeg = 6.4; + + else if (run == 3170) thetadeg = 6.4; + else if (run == 3171) thetadeg = 6.4; + else if (run == 3172) thetadeg = 6.5; + else if (run == 3173) thetadeg = 6.6; + else if (run == 3174) thetadeg = 6.7; + + else if (run == 3176) thetadeg = 7.1; + else if (run == 3177) thetadeg = 7.4; + else if (run == 3178) thetadeg = 7.6; + else if (run == 3179) thetadeg = 7.9; + + + else if (run == 3182) thetadeg = 8.4; + else if (run == 3183) thetadeg = 8.9; + else if (run == 3184) thetadeg = 9.2; + else if (run == 3185) thetadeg = 9.5; + else if (run == 3186) thetadeg = 9.8; + else if (run == 3187) thetadeg = 10.5; + else if (run == 3188) thetadeg = 10.9; + else if (run == 3189) thetadeg = 11.2; + + else if (run == 3190) thetadeg = 11.6; + else if (run == 3191) thetadeg = 11.6; + else if (run == 3192) thetadeg = 12.4; + else if (run == 3193) thetadeg = 12.7; + else if (run == 3194) thetadeg = 13.1; + else if (run == 3195) thetadeg = 13.5; + else if (run == 3196) thetadeg = 13.9; + else if (run == 3197) thetadeg = 14.3; + else if (run == 3198) thetadeg = 14.7; + else if (run == 3199) thetadeg = 15.1; + + else if (run == 3200) thetadeg = 15.6; + else if (run == 3201) thetadeg = 16.0; + else if (run == 3202) thetadeg = 16.5; + else if (run == 3203) thetadeg = 16.9; + else if (run == 3204) thetadeg = 17.3; + else if (run == 3205) thetadeg = 17.7; + else if (run == 3206) thetadeg = 18.2; + else if (run == 3207) thetadeg = 18.6; + else if (run == 3208) thetadeg = 19.0; + else if (run == 3209) thetadeg = 19.4; + + else if (run == 3210) thetadeg = 19.9; + else if (run == 3211) thetadeg = 20.4; + else if (run == 3212) thetadeg = 20.8; + else if (run == 3213) thetadeg = 21.2; + else if (run == 3214) thetadeg = 21.7; + else if (run == 3215) thetadeg = 22.2; + else if (run == 3216) thetadeg = 25.6; // Off Crab1 + else if (run == 3217) thetadeg = 25.0; + else if (run == 3218) thetadeg = 24.5; + else if (run == 3219) thetadeg = 24.0; + + else if (run == 3220) thetadeg = 23.5; + else if (run == 3221) thetadeg = 22.5; + else if (run == 3222) thetadeg = 22.1; + + else if (run == 3225) thetadeg = 15.1; + else if (run == 3226) thetadeg = 15.0; + else if (run == 3227) thetadeg = 14.5; + else if (run == 3228) thetadeg = 14.1; + else if (run == 3229) thetadeg = 13.8; + + else if (run == 3230) thetadeg = 13.3; + else if (run == 3231) thetadeg = 13.0; + else if (run == 3232) thetadeg = 12.6; + else if (run == 3233) thetadeg = 12.3; + else if (run == 3234) thetadeg = 12.0; + else if (run == 3235) thetadeg = 11.6; + else if (run == 3236) thetadeg = 11.3; + else if (run == 3237) thetadeg = 11.0; + else if (run == 3238) thetadeg = 10.8; + else if (run == 3239) thetadeg = 10.4; + + else if (run == 3240) thetadeg = 10.1; + else if (run == 3241) thetadeg = 9.9; + else if (run == 3242) thetadeg = 9.6; + else if (run == 3243) thetadeg = 9.4; + else if (run == 3244) thetadeg = 9.2; + else if (run == 3245) thetadeg = 9.0; + else if (run == 3246) thetadeg = 8.9; + else if (run == 3247) thetadeg = 8.8; + else if (run == 3248) thetadeg = 8.7; + else if (run == 3249) thetadeg = 8.6; + + else if (run == 3250) thetadeg = 8.6; + else if (run == 3251) thetadeg = 8.6; + else if (run == 3252) thetadeg = 8.6; + else if (run == 3253) thetadeg = 8.7; + else if (run == 3254) thetadeg = 8.8; + else if (run == 3255) thetadeg = 8.9; + else if (run == 3256) thetadeg = 9.1; + else if (run == 3257) thetadeg = 9.3; + else if (run == 3258) thetadeg = 9.5; + else if (run == 3259) thetadeg = 9.7; + else if (run == 3260) thetadeg = 9.9; + + else if (run == 3261) thetadeg = 10.2; + else if (run == 3262) thetadeg = 10.5; + else if (run == 3263) thetadeg = 10.8; + else if (run == 3264) thetadeg = 11.1; + else if (run == 3265) thetadeg = 11.4; + else if (run == 3266) thetadeg = 11.8; + else if (run == 3267) thetadeg = 12.1; + else if (run == 3268) thetadeg = 12.5; + else if (run == 3269) thetadeg = 12.8; + + else if (run == 3270) thetadeg = 13.2; + else if (run == 3271) thetadeg = 13.5; + else if (run == 3272) thetadeg = 13.9; + else if (run == 3273) thetadeg = 14.0; + else if (run == 3274) thetadeg = 14.4; + + else if (run == 3284) thetadeg = 7.0; // Crab + else if (run == 3285) thetadeg = 7.2; + else if (run == 3286) thetadeg = 7.4; + else if (run == 3287) thetadeg = 7.5; + else if (run == 3288) thetadeg = 8.4; + else if (run == 3289) thetadeg = 9.0; + + else if (run == 3290) thetadeg = 9.4; + else if (run == 3291) thetadeg = 9.8; + else if (run == 3292) thetadeg = 10.2; + else if (run == 3293) thetadeg = 10.5; + else if (run == 3294) thetadeg = 10.9; + else if (run == 3295) thetadeg = 11.3; + else if (run == 3296) thetadeg = 11.8; + else if (run == 3297) thetadeg = 12.2; + else if (run == 3298) thetadeg = 12.6; + else if (run == 3299) thetadeg = 13.0; + + else if (run == 3300) thetadeg = 13.5; + else if (run == 3301) thetadeg = 13.9; + else if (run == 3302) thetadeg = 14.3; + else if (run == 3303) thetadeg = 14.8; + else if (run == 3304) thetadeg = 15.2; + else if (run == 3305) thetadeg = 15.7; + else if (run == 3306) thetadeg = 16.2; + else if (run == 3307) thetadeg = 16.6; + else if (run == 3308) thetadeg = 17.1; + else if (run == 3309) thetadeg = 17.6; + + else if (run == 3310) thetadeg = 17.9; + else if (run == 3311) thetadeg = 18.4; + else if (run == 3312) thetadeg = 18.9; + else if (run == 3313) thetadeg = 19.3; + else if (run == 3314) thetadeg = 19.8; + else if (run == 3315) thetadeg = 20.1; + else if (run == 3316) thetadeg = 20.7; + else if (run == 3317) thetadeg = 21.2; + else if (run == 3318) thetadeg = 21.7; + else if (run == 3319) thetadeg = 22.1; + + else if (run == 3320) thetadeg = 22.5; + else if (run == 3321) thetadeg = 23.1; + else if (run == 3322) thetadeg = 23.6; + else if (run == 3323) thetadeg = 24.1; + else if (run == 3324) thetadeg = 24.6; + else if (run == 3325) thetadeg = 24.9; + else if (run == 3326) thetadeg = 25.5; + else if (run == 3327) thetadeg = 26.0; + else if (run == 3328) thetadeg = 26.0; + else if (run == 3329) thetadeg = 26.6; + + else if (run == 3330) thetadeg = 26.6; + else if (run == 3331) thetadeg = 27.1; + else if (run == 3332) thetadeg = 27.7; + else if (run == 3333) thetadeg = 28.2; + else if (run == 3334) thetadeg = 28.5; + + else if (run == 3340) thetadeg = 10.5; + else if (run == 3341) thetadeg = 10.3; + else if (run == 3342) thetadeg = 9.6; + else if (run == 3343) thetadeg = 9.2; + else if (run == 3344) thetadeg = 8.9; + else if (run == 3345) thetadeg = 8.6; + else if (run == 3346) thetadeg = 8.3; + else if (run == 3347) thetadeg = 8.0; + else if (run == 3348) thetadeg = 7.7; + else if (run == 3349) thetadeg = 7.5; + + else if (run == 3350) thetadeg = 7.2; + else if (run == 3351) thetadeg = 7.0; + else if (run == 3352) thetadeg = 6.8; + else if (run == 3353) thetadeg = 6.7; + else if (run == 3354) thetadeg = 6.6; + else if (run == 3355) thetadeg = 6.5; + else if (run == 3356) thetadeg = 6.4; + else if (run == 3357) thetadeg = 6.4; + else if (run == 3358) thetadeg = 6.4; + else if (run == 3359) thetadeg = 6.5; + + else if (run == 3360) thetadeg = 6.6; + + else if (run == 3362) thetadeg = 6.7; + else if (run == 3363) thetadeg = 6.8; + else if (run == 3364) thetadeg = 7.0; + else if (run == 3365) thetadeg = 7.2; + else if (run == 3366) thetadeg = 7.5; + else if (run == 3367) thetadeg = 7.7; + else if (run == 3368) thetadeg = 8.0; + else if (run == 3369) thetadeg = 8.3; + + else if (run == 3370) thetadeg = 8.6; + else if (run == 3371) thetadeg = 9.0; + else if (run == 3372) thetadeg = 9.3; + else if (run == 3373) thetadeg = 9.6; + else if (run == 3374) thetadeg = 10.0; + else if (run == 3375) thetadeg = 10.4; + else if (run == 3376) thetadeg = 10.7; + else if (run == 3377) thetadeg = 11.1; + else if (run == 3378) thetadeg = 11.5; + else if (run == 3379) thetadeg = 11.9; + + else if (run == 3380) thetadeg = 12.3; + else if (run == 3381) thetadeg = 12.7; + else if (run == 3382) thetadeg = 13.1; + else if (run == 3383) thetadeg = 13.5; + else if (run == 3384) thetadeg = 13.9; + else if (run == 3385) thetadeg = 14.3; + else if (run == 3386) thetadeg = 14.7; + else if (run == 3387) thetadeg = 15.2; + else if (run == 3388) thetadeg = 15.6; + else if (run == 3389) thetadeg = 16.0; + + else if (run == 3390) thetadeg = 16.4; + else if (run == 3391) thetadeg = 16.7; + else if (run == 3392) thetadeg = 17.9; + else if (run == 3393) thetadeg = 18.3; + else if (run == 3394) thetadeg = 18.7; + else if (run == 3395) thetadeg = 19.2; + else if (run == 3396) thetadeg = 19.6; + else if (run == 3397) thetadeg = 20.0; + else if (run == 3398) thetadeg = 20.5; + else if (run == 3399) thetadeg = 20.9; + + else if (run == 3400) thetadeg = 21.4; + else if (run == 3401) thetadeg = 21.8; + else if (run == 3402) thetadeg = 22.1; + else if (run == 3403) thetadeg = 22.6; + else if (run == 3404) thetadeg = 23.1; + else if (run == 3405) thetadeg = 23.4; + else if (run == 3406) thetadeg = 23.9; + else if (run == 3407) thetadeg = 24.3; + else if (run == 3408) thetadeg = 24.7; + else if (run == 3409) thetadeg = 26.9; + + else if (run == 3410) thetadeg = 27.3; + else if (run == 3411) thetadeg = 27.7; + else if (run == 3412) thetadeg = 28.2; + else if (run == 3413) thetadeg = 28.7; + else if (run == 3414) thetadeg = 29.1; + else if (run == 3415) thetadeg = 29.2; + else if (run == 3416) thetadeg = 30.0; + else if (run == 3417) thetadeg = 18.5; // Off Crab1 + else if (run == 3418) thetadeg = 18.4; + else if (run == 3419) thetadeg = 17.5; + + else if (run == 3420) thetadeg = 17.2; + else if (run == 3421) thetadeg = 16.8; + else if (run == 3422) thetadeg = 16.4; + else if (run == 3423) thetadeg = 16.0; + else if (run == 3424) thetadeg = 15.6; + else if (run == 3425) thetadeg = 15.3; + else if (run == 3426) thetadeg = 14.9; + else if (run == 3427) thetadeg = 14.5; + else if (run == 3428) thetadeg = 14.1; + else if (run == 3429) thetadeg = 13.7; + + else if (run == 3430) thetadeg = 13.4; + else if (run == 3431) thetadeg = 13.0; + else if (run == 3432) thetadeg = 12.7; + else if (run == 3433) thetadeg = 12.3; + else if (run == 3434) thetadeg = 12.0; + else if (run == 3435) thetadeg = 12.0; + else if (run == 3436) thetadeg = 11.6; + else if (run == 3437) thetadeg = 11.3; + else if (run == 3438) thetadeg = 11.0; + else if (run == 3439) thetadeg = 10.7; + + else if (run == 3440) thetadeg = 10.4; + else if (run == 3441) thetadeg = 10.1; + else if (run == 3442) thetadeg = 9.9; + else if (run == 3443) thetadeg = 9.8; + else if (run == 3444) thetadeg = 30.8; // Mkn 421 + else if (run == 3445) thetadeg = 30.6; + else if (run == 3446) thetadeg = 29.7; + else if (run == 3447) thetadeg = 29.3; + else if (run == 3448) thetadeg = 28.9; + else if (run == 3449) thetadeg = 28.5; + + else if (run == 3450) thetadeg = 28.1; + else if (run == 3451) thetadeg = 27.7; + else if (run == 3452) thetadeg = 27.3; + else if (run == 3453) thetadeg = 26.9; + else if (run == 3454) thetadeg = 26.5; + else if (run == 3455) thetadeg = 26.1; + else if (run == 3456) thetadeg = 25.7; + else if (run == 3457) thetadeg = 25.3; + else if (run == 3458) thetadeg = 24.9; + else if (run == 3459) thetadeg = 24.5; + + else if (run == 3460) thetadeg = 24.1; + + else if (run == 3463) thetadeg = 23.2; + else if (run == 3464) thetadeg = 22.8; + else if (run == 3465) thetadeg = 22.4; + else if (run == 3466) thetadeg = 22.0; + else if (run == 3467) thetadeg = 21.6; + else if (run == 3468) thetadeg = 21.2; + else if (run == 3469) thetadeg = 20.8; + + else if (run == 3470) thetadeg = 20.4; + else if (run == 3471) thetadeg = 20.1; + else if (run == 3472) thetadeg = 19.7; + else if (run == 3473) thetadeg = 19.3; + else if (run == 3474) thetadeg = 18.9; + else if (run == 3475) thetadeg = 18.5; + else if (run == 3476) thetadeg = 18.2; + else if (run == 3477) thetadeg = 18.1; + else if (run == 3478) thetadeg = 17.7; + + + else if (run == 3480) thetadeg = 17.5; + else if (run == 3481) thetadeg = 17.1; + else if (run == 3482) thetadeg = 16.7; + else if (run == 3483) thetadeg = 16.3; + else if (run == 3484) thetadeg = 16.0; + else if (run == 3485) thetadeg = 15.6; + else if (run == 3486) thetadeg = 15.3; + else if (run == 3487) thetadeg = 15.0; + else if (run == 3488) thetadeg = 14.8; + else if (run == 3489) thetadeg = 14.8; + + else + { + *fLog << warn << "MPedestalWorkaround : no zenith angle for run number = " + << run << endl; + thetadeg = 90.0; + } + + thetarad = thetadeg / kRad2Deg; + fMcEvt->SetTelescopeTheta(thetarad); + + return kTRUE; } @@ -113,2 +547,8 @@ } + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2745) +++ trunk/MagicSoft/Mars/manalysis/MPedestalWorkaround.h (revision 2746) @@ -19,5 +19,6 @@ class MGeomCam; class MParList; - +class MMcEvt; +class MRawRunHeader; class MPedestalWorkaround : public MTask @@ -25,7 +26,9 @@ private: - MPedestalCam *fPed; // - MPedPhotCam *fPedPhot; // - MGeomCam *fCam; // + MPedestalCam *fPed; + MPedPhotCam *fPedPhot; + MGeomCam *fCam; + MMcEvt *fMcEvt; + MRawRunHeader *fRun; Int_t PreProcess(MParList *pList); Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2745) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2746) @@ -74,4 +74,5 @@ MFiltercutsCalc.cc \ MCT1PadONOFF.cc \ + MPadON.cc \ MPadONOFF.cc \ MPedestalWorkaround.cc \ Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 2745) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 2746) @@ -142,4 +142,8 @@ void SetImpact(Float_t Impact) { fImpact=Impact;} //Set impact parameter + + void SetTelescopeTheta(Float_t Theta) { fTelescopeTheta=Theta; } + void SetTelescopePhi(Float_t Phi) { fTelescopePhi=Phi; } + /* void SetPhotIni(Short_t PhotIni) */