Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 1808) +++ trunk/MagicSoft/Mars/Changelog (revision 1809) @@ -1,3 +1,49 @@ -*-*- END OF LINE -*-*- + + + 2003/03/08: Wolfgang Wittek + + * macros/AnalyseCT1.C + - for the CT1 data analysis + + * mhist/MHMatrix.[h,cc] + - let refcolumn start at 1 (not at 0) + + * mhist/MHSigmaTheta.[h,cc] + - Draw replaced by DrawCopy + - add SetDirectory(NULL) + + * manalysis/ MSelBasic.[h,cc] + MSelStandard.[h,cc] + MSelFinal.[h,cc] + - more detailed output for errors + - bugs removed + + * manalysis/ MPadSchweizer.[h,cc] + - add SetDirectory(NULL) + - add fErrors + + * mfilter/ MFEventSelector.[h,cc] + - add fErrors + + * manalysis/MMultiDimDistCalc.[h,cc] + - check division by zero + + * mhist/MHHadronness.[h,cc] + - check division by zero + - normalize distributions of hadronness + + * mfileio/MCT1ReadPreProc.[h,cc] + - add event number (event.isecs_since_midday) + - change definition of "fIsMcFile", + because outpars.bmontecarlo is set wrongly sometimes + - copy pedestalRMS for each event from the header information + - check for the presence of a footer record even after reading + a run header + + * mmc/ MMcEvt.[hxx,cxx]] + - add GetEvtNumber() + + 2003/02/27: Abelardo Moralejo @@ -63,5 +109,4 @@ * mmc/MMcTrigHeader.hxx - Added GetMultiplicity() and GetMeanThreshold() - Index: trunk/MagicSoft/Mars/macros/AnalyseCT1.C =================================================================== --- trunk/MagicSoft/Mars/macros/AnalyseCT1.C (revision 1808) +++ trunk/MagicSoft/Mars/macros/AnalyseCT1.C (revision 1809) @@ -1,51 +1,112 @@ void AnalyseCT1() { - const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; - const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; + //----------------------------------------------- + //const char *offfile = "/hd43/rwagner/CT1/PreprocData/offdata.preproc"; + //const char *onfile = "/hd43/rwagner/CT1/PreprocData/mkn421_on.preproc"; //const char *mcfile = "/hd43/rwagner/CT1/PreprocData/mc_Gammas.preproc.0.6"; //const char *mcfile = "/hd43/rwagner/mkn421_mc.preproc"; - const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc"; + //const char *mcfile = "/hd43/rwagner/mkn421_mc_pedrms_0.001.preproc"; + + //----------------------------------------------- + const char *offfile = "~/Mars200203/CT1data/offdata.preproc"; + + const char *onfile = "~/Mars200203/CT1data/mkn421_on.preproc"; + //const char *onfile = "~/Mars200203/CT1data/Crab_preproc_daniel_0.6"; + + //const char *mcfile = "~/Mars200203/CT1data/mc_gammas.preproc"; //Version 0.5 + //const char *mcfile = "~/Mars200203/CT1data/mc_Gammas.preproc.0.6"; + const char *mcfile = "~/Mars200203/CT1data/mkn421_mc_pedrms_0.001.preproc"; + + //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6"; + //const char *mcfile = "~/Mars200203/CT1data/mc_preproc_daniel_0.6_040303"; + + //----------------------------------------------- + //const char *offfile = "~magican/home/ct1test/PreprocData/offdata.preproc"; + //const char *onfile = "~magican/home/ct1test/PreprocData/mkn421_on.preproc"; + //const char *mcfile = "~magican/home/ct1test/PreprocData/mc_gammas.preproc"; + + const char *filename; - Bool_t Data = kFALSE; // loop over data ? - Bool_t WLook = kFALSE; // write out look-alike events for hadrons ? - Bool_t WHist = kFALSE; // write out histograms for padding ? - - Bool_t MC = kTRUE; // loop over MC gamma data ? - Bool_t RHist = kTRUE; // read in histograms for padding ? - Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? - - - cout << "" << endl; - cout << "Macro AnalyseCT1 : Data, WHist = " << Data << ", " - << WHist << endl; - - cout << " RHist, MC = " << RHist << ", " - << MC << endl; - cout << "" << endl; + //************************************************************************ + // Job 1 : read OFF data + // (generate sigmabar vs. Theta plot; write root file for OFF data (OFF1); + // generate hadron matrix for g/h separation) + + Bool_t Job1 = kFALSE; + Bool_t WHistOFF = kTRUE; // write out histogram sigmabar vs. Theta ? + Bool_t WLookOFF = kTRUE; // write out look-alike events for hadrons ? + Bool_t WOFF1 = kTRUE; // write out root file OFF1 ? + + + // Job 2 : read ON data + // (generate sigmabar vs. Theta plot; write root file for ON data (ON1)) + + Bool_t Job2 = kFALSE; + Bool_t WHistON = kFALSE; // write out histogram sigmabar vs. Theta ? + Bool_t WLookON = kFALSE; // write out look-alike events for hadrons ? + Bool_t WON1 = kFALSE; // write out root file ON1 ? + + // Job 3 : read MC gamma data, read sigmabar vs. Theta plot from OFF data + // (do padding; write root file for MC gammas (MC1); + // generate gamma matrix for g/h separation) + + Bool_t Job3 = kFALSE; + Bool_t WLookMC = kTRUE; // write out look-alike events for gammas ? + Bool_t WMC1 = kTRUE; // write out root file MC1 ? + + + + // Job 4 : read OFF1 and MC1 files, read hadron and gamma matrix + // (produce the Neyman-Pearson plot) + Bool_t Job4 = kTRUE; + + + // Job 5 : read MC1 file, read hadron and gamma matrix + // (do g/h separation; write root file for MC gammas after final cuts (MC2)) + Bool_t Job5 = kFALSE; + Bool_t WMC2 = kFALSE; // write out root file MC2 ? + + + // Job 6 : read ON data, read hadron and gamma matrix + // (do g/h separation; write root file for ON data after final cuts (ON2)) + Bool_t Job6 = kFALSE; + Bool_t WON2 = kTRUE; // write out root file ON2 ? + //************************************************************************ + + + //--------------------------------------------------------------------- - // start of section for ON data - // (in the CT1 analysis the ON data are also used as a sample representing - // the hadrons for the g/h separation) - - // read ON data file - - // - to produce the 2D-histogram "sigmabar versus Theta" + // Job 1 + //====== + + // read OFF data file + + // - produce the 2D-histogram "sigmabar versus Theta" for OFF data // (to be used for the padding of the MC gamma data) - // - to write a file of look-alike events (for hadrons) + // - write a file of look-alike events (hadron matrix) // (to be used later together with the corresponding MC gamma file // for the g/h separation in the analysis of the ON data) - // - to write a file of ON events (after the standard cuts, - // before the g/h separation) + // - write a file of OFF events (OFF1) + // (after the standard cuts, before the g/h separation) // (to be used together with the corresponding MC gamma file // for the optimization of the g/h separation) - if (Data) + if (Job1) { cout << "=====================================================" << endl; - cout << "Macro AnalyseCT1 : Start of section for ON data" << endl; + cout << "Macro AnalyseCT1 : Start of Job 1" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job1, WHistOFF, WLookOFF, WOFF1 = " + << Job1 << ", " << WHistOFF << ", " << WLookOFF << ", " + << WOFF1 << endl; + + + TString typeData = "OFF"; + cout << "typeData = " << typeData << endl; MTaskList tliston; @@ -112,7 +173,7 @@ // - filename = onfile; - readname = "ReadCT1data"; - MCT1ReadPreProc read(filename, readname); + filename = offfile; + RName = "ReadCT1data"; + MCT1ReadPreProc read(filename, RName); MSelBasic selbasic; @@ -150,5 +211,5 @@ MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); - MSelStandard selstand(fHillas, fHillasSrc); + MSelStandard selstand(fHilName, fHilSrcName); MFillH hfill1("MHHillas", fHilName); @@ -168,6 +229,8 @@ const char* mtxName = "MatrixHadrons"; - Int_t howMany = 30000; - char* outName = "matrix_hadrons.root"; + Int_t howMany = 50000; + TString outName = "matrix_hadrons_"; + outName += typeData; + outName += ".root"; cout << "" << endl; @@ -182,5 +245,5 @@ // matrix limitation for look-alike events (approximate number) - MFEventSelector selector("", "", readname); + MFEventSelector selector("", "", RName); selector.SetNumSelectEvts(howMany); MFilterList flist; @@ -192,18 +255,16 @@ MHMatrix *pmatrix = new MHMatrix(mtxName); MHMatrix& matrix = *pmatrix; + + matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrix.AddColumn("MSigmabar.fSigmabar"); + matrix.AddColumn("log10(Hillas.fSize)"); + matrix.AddColumn("HillasSrc.fDist"); matrix.AddColumn("Hillas.fWidth"); matrix.AddColumn("Hillas.fLength"); - matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); - matrix.AddColumn("abs(Hillas.fAsym)"); + matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); matrix.AddColumn("abs(Hillas.fM3Long)"); - matrix.AddColumn("abs(Hillas.fM3Trans)"); - matrix.AddColumn("abs(HillasSrc.fHeadTail)"); matrix.AddColumn("Hillas.fConc"); matrix.AddColumn("Hillas.fConc1"); - matrix.AddColumn("HillasSrc.fDist"); - matrix.AddColumn("log10(Hillas.fSize)"); - matrix.AddColumn("HillasSrc.fAlpha"); - matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); - //matrix.AddColumn("MMcEvt.fTheta"); + pliston.AddToList(pmatrix); @@ -211,7 +272,21 @@ fillmatg.SetFilter(&flist); + + if (WOFF1) + { + TString outNameImage = typeData; + outNameImage += "1.root"; + MWriteRootFile write(outNameImage); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("Hillas", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("HillasSrc", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MMcRunHeader","RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + } + //+++++++++++++++++++++++++++++++++++++++++++++++++++ - MSelFinal selfinal(fHillas, fHillasSrc); //***************************** @@ -229,4 +304,10 @@ tliston.AddToList(&selstand); + if (WOFF1) + tliston.AddToList(&write); + + tliston.AddToList(&flist); + tliston.AddToList(&fillmatg); + tliston.AddToList(&hfill1); tliston.AddToList(&hfill2); @@ -234,8 +315,4 @@ tliston.AddToList(&hfill4); - tliston.AddToList(&flist); - tliston.AddToList(&fillmatg); - - //tliston.AddToList(&selfinal); //***************************** @@ -243,8 +320,11 @@ // Execute event loop // + MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); Int_t maxevents = 1000000000; + //Int_t maxevents = 1000; if ( !evtloop.Eventloop(maxevents) ) return; @@ -255,11 +335,13 @@ // Display the histograms // - pliston.FindObject("MHHillas")->DrawClone(); - pliston.FindObject("MHHillasSrc")->DrawClone(); + TObject *c; + c = pliston.FindObject("MHHillas")->DrawClone(); + + c = pliston.FindObject("MHHillasSrc")->DrawClone(); //pliston.FindObject("MHHillasExt")->DrawClone(); - pliston.FindObject(nxt)->DrawClone(); - - pliston.FindObject("MHStarMap")->DrawClone(); + c = pliston.FindObject(nxt)->DrawClone(); + + c = pliston.FindObject("MHStarMap")->DrawClone(); @@ -277,7 +359,7 @@ Int_t nmaxevts = 2000; - // target distribution for the variable in column refcolumn; + // target distribution for the variable in column refcolumn (start at 1); // set refcolumn negative if distribution is not to be changed - Int_t refcolumn = 12; + Int_t refcolumn = 1; Int_t nbins = 10; Float_t frombin = 0.5; @@ -311,5 +393,5 @@ // write out look-alike events (for hadrons) // - if (WLook) + if (WLookOFF) { TFile *writejens = new TFile(outName, "RECREATE", ""); @@ -324,5 +406,5 @@ //------------------------------------------- // Write histograms onto a file - if (WHist) + if (WHistOFF) { MHSigmaTheta *sigtheta = @@ -330,5 +412,5 @@ if (!sigtheta) { - *fLog << "Object 'SigmaTheta' not found" << endl; + cout << "Object 'SigmaTheta' not found" << endl; return; } @@ -337,5 +419,9 @@ TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); - TFile outfile("SigmaTheta.root", "RECREATE"); + TString outNameSigTh = "SigmaTheta_"; + outNameSigTh += typeData; + outNameSigTh += ".root"; + + TFile outfile(outNameSigTh, "RECREATE"); fHSigTh->Write(); fHSigPixTh->Write(); @@ -343,9 +429,9 @@ outfile.Close(); - cout << "File 'SigmaTheta.root' was written out" << endl; + cout << "File " << outNameSigTh << " was written out" << endl; } - cout << "Macro AnalyseCT1 : End of section for ON data" << endl; + cout << "Macro AnalyseCT1 : End of Job 1" << endl; cout << "===================================================" << endl; } @@ -353,5 +439,366 @@ //--------------------------------------------------------------------- - // start of section for MC gamma data + // Job 2 + //====== + // read ON data file + + // - produce the 2D-histogram "sigmabar versus Theta" for ON data + // (to be used for the padding of the MC gamma data) + + // - write a file of look-alike events (hadron matrix) + // (to be used later together with the corresponding MC gamma file + // for the g/h separation in the analysis of the ON data) + + // - write a file of ON events (ON1) + // (after the standard cuts, before the g/h separation) + // (to be used together with the corresponding MC gamma file + // for the optimization of the g/h separation) + + if (Job2) + { + cout << "=====================================================" << endl; + cout << "Macro AnalyseCT1 : Start of Job 2" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job2, WHistON, WLookON, WON1 = " + << Job2 << ", " << WHistON << ", " << WLookON << ", " + << WON1 << endl; + + TString typeData = "ON"; + cout << "typeData = " << typeData << endl; + + MTaskList tliston; + MParList pliston; + pliston.AddToList(&tliston); + + + //:::::::::::::::::::::::::::::::::::::::::: + + MBinning binssize("BinningSize"); + binssize.SetEdgesLog(50, 10, 1.0e5); + pliston.AddToList(&binssize); + + MBinning binsdistc("BinningDist"); + binsdistc.SetEdges(50, 0, 1.4); + pliston.AddToList(&binsdistc); + + MBinning binswidth("BinningWidth"); + binswidth.SetEdges(50, 0, 1.0); + pliston.AddToList(&binswidth); + + MBinning binslength("BinningLength"); + binslength.SetEdges(50, 0, 1.0); + pliston.AddToList(&binslength); + + MBinning binsalpha("BinningAlpha"); + binsalpha.SetEdges(100, -100, 100); + pliston.AddToList(&binsalpha); + + MBinning binsasym("BinningAsym"); + binsasym.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsasym); + + MBinning binsht("BinningHeadTail"); + binsht.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsht); + + MBinning binsm3l("BinningM3Long"); + binsm3l.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3l); + + MBinning binsm3t("BinningM3Trans"); + binsm3t.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3t); + + + //..... + MBinning binsb("BinningSigmabar"); + binsb.SetEdges( 100, 0.0, 5.0); + pliston.AddToList(&binsb); + + MBinning binth("BinningTheta"); + binth.SetEdges( 70, -0.5, 69.5); + pliston.AddToList(&binth); + + MBinning binsdiff("BinningDiffsigma2"); + binsdiff.SetEdges(100, -5.0, 20.0); + pliston.AddToList(&binsdiff); + //:::::::::::::::::::::::::::::::::::::::::: + + + //------------------------------------------- + // create the tasks which should be executed + // + + filename = onfile; + RName = "ReadCT1data"; + MCT1ReadPreProc read(filename, RName); + + MSelBasic selbasic; + + MBlindPixelCalc blind; + blind.SetUseInterpolation(); + blind.SetUseBlindPixels(); + // blind.SetUseCetralPixel(); + + // create an object of class MSigmabar, + // because class MHSigmaTheta will look for it + MSigmabar sigbar; + pliston.AddToList(&sigbar); + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + fillsigtheta.SetTitle("Task to make 2D plot Sigmabar vs Theta"); + + MImgCleanStd clean; + + // calculate also extended image parameters + TString fHilName = "Hillas"; + MHillasExt hext(fHilName); + pliston.AddToList(&hext); + MHillasExt *fHillas = &hext; + + // name of output container for MHillasCalc + // = name of MHillas object + MHillasCalc hcalc(fHilName); + + // name of output container for MHillasSrcCalc + // = name of MHillasSrc object + TString fHilSrcName = "HillasSrc"; + MHillasSrc hilsrc(fHilSrcName); + MHillasSrc *fHillasSrc = &hilsrc; + pliston.AddToList(&hilsrc); + MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); + + MSelStandard selstand(fHilName, fHilSrcName); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + TString nxt = "HillasExt"; + MHHillasExt hhilext(nxt, "", fHilName); + pliston.AddToList(&hhilext); + TString namext = nxt; + namext += "[MHHillasExt]"; + MFillH hfill3(namext, fHilSrcName); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + // prepare writing of look-alike events (for hadrons) + + const char* mtxName = "MatrixHadrons"; + Int_t howMany = 50000; + TString outName = "matrix_hadrons_"; + outName += typeData; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Matrix for (hadrons)" << endl; + cout << "input file = " << filename<< endl; + cout << "matrix name = " << mtxName << endl; + cout << "max. no.of look-alike events = " << howMany << endl; + cout << "name of output root file = " << outName << endl; + cout << "" << endl; + + + // matrix limitation for look-alike events (approximate number) + MFEventSelector selector("", "", RName); + selector.SetNumSelectEvts(howMany); + MFilterList flist; + flist.AddToList(&selector); + + // + // --- setup the matrix and add it to the parlist + // + MHMatrix *pmatrix = new MHMatrix(mtxName); + MHMatrix& matrix = *pmatrix; + + matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrix.AddColumn("MSigmabar.fSigmabar"); + matrix.AddColumn("log10(Hillas.fSize)"); + matrix.AddColumn("HillasSrc.fDist"); + matrix.AddColumn("Hillas.fWidth"); + matrix.AddColumn("Hillas.fLength"); + matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); + matrix.AddColumn("abs(Hillas.fM3Long)"); + matrix.AddColumn("Hillas.fConc"); + matrix.AddColumn("Hillas.fConc1"); + + pliston.AddToList(pmatrix); + + MFillH fillmatg(mtxName); + fillmatg.SetFilter(&flist); + + + if (WON1) + { + TString outNameImage = typeData; + outNameImage += "1.root"; + MWriteRootFile write(outNameImage); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("Hillas", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("HillasSrc", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MMcRunHeader","RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + } + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ + + + //***************************** + // tasks to be executed + + tliston.AddToList(&read); + + tliston.AddToList(&selbasic); + tliston.AddToList(&blind); + tliston.AddToList(&fillsigtheta); + tliston.AddToList(&clean); + + tliston.AddToList(&hcalc); + tliston.AddToList(&csrc1); + + tliston.AddToList(&selstand); + if (WON1) + tliston.AddToList(&write); + + tliston.AddToList(&flist); + tliston.AddToList(&fillmatg); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + //------------------------------------------- + // Display the histograms + // + TObject *c; + c = pliston.FindObject("MHHillas")->DrawClone(); + + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + + //pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject(nxt)->DrawClone(); + + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + // + // ---------------------------------------------------------- + // Definition of the reference sample of look-alike events + // (this is a subsample of the original sample) + // ---------------------------------------------------------- + // + cout << "" << endl; + cout << "========================================================" << endl; + // Select a maximum of nmaxevts events from the sample of look-alike + // events. They will form the reference sample. + Int_t nmaxevts = 2000; + + // target distribution for the variable in column refcolumn (start at 1); + // set refcolumn negative if distribution is not to be changed + Int_t refcolumn = 1; + Int_t nbins = 10; + Float_t frombin = 0.5; + Float_t tobin = 1.0; + TH1F *thsh = new TH1F("thsh","target distribution", + nbins, frombin, tobin); + thsh->SetXTitle("cos( \\Theta)"); + thsh->SetYTitle("Counts"); + Float_t dbin = (tobin-frombin)/nbins; + Float_t lbin = frombin +dbin*0.5; + for (Int_t j=1; j<=nbins; j++) + { + thsh->Fill(lbin,1.0); + lbin += dbin; + } + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Macro AnalyseCT1 : define reference sample for hadrons" << endl; + cout << "Macro AnalyseCT1 : Parameters for reference sample : refcolumn, nmaxevts = " + << refcolumn << ", " << nmaxevts << endl; + + if ( !matrix.DefRefMatrix(refcolumn, thsh, nmaxevts) ) + { + cout << "Macro AnalyseCT1 : Reference matrix for hadrons cannot be defined" << endl; + return; + } + //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + //------------------------------------------- + // write out look-alike events (for hadrons) + // + if (WLookON) + { + TFile *writejens = new TFile(outName, "RECREATE", ""); + matrix.Write(); + cout << "Macro AnalyseCT1 : matrix of look-alike events for hadrons written onto file " + << outName << endl; + + writejens->Close(); + delete writejens; + } + + //------------------------------------------- + // Write histograms onto a file + if (WHistON) + { + MHSigmaTheta *sigtheta = + (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); + if (!sigtheta) + { + cout << "Object 'SigmaTheta' not found" << endl; + return; + } + TH2D *fHSigTh = sigtheta->GetSigmaTheta(); + TH3D *fHSigPixTh = sigtheta->GetSigmaPixTheta(); + TH3D *fHDifPixTh = sigtheta->GetDiffPixTheta(); + + TString outNameSigTh = "SigmaTheta_"; + outNameSigTh += typeData; + outNameSigTh += ".root"; + + TFile outfile(outNameSigTh, "RECREATE"); + fHSigTh->Write(); + fHSigPixTh->Write(); + fHDifPixTh->Write(); + outfile.Close(); + + cout << "File " << outNameSigTh << " was written out" << endl; + } + + + cout << "Macro AnalyseCT1 : End of Job 2" << endl; + cout << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job 3 + //====== // read MC gamma data @@ -360,28 +807,31 @@ // (using the 2D-histogram "sigmabar versus Theta" of the OFF data) - // - to write a file of look-alike events (for gammas) + // - to write a file of look-alike events (gammas matrix) // (to be used later together with the corresponding hadron file // for the g/h separation in the analysis of the ON data) - // - to write a file of padded MC gamma events + // - to write a file of padded MC gamma events (MC1) // (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) - // - to write a file of padded MC gamma events (after the final cuts) - // (to be used for the optimization of the energy estimator - // and for calculating the effective collection areas) - - if (MC) + + if (Job3) { - cout << "=============================================================" - << endl; - cout << "Macro : AnalyseCT1 : Start of section for MC gamma data" - << endl; - - MTaskList tlist; - MParList plist; - - plist.AddToList(&tlist); + cout << "=====================================================" << endl; + cout << "Macro AnalyseCT1 : Start of Job 3" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job3, WLookMC, WMC1 = " + << Job3 << ", " << WLookMC << ", " << WMC1 << endl; + + + TString typeMC = "MC"; + cout << "typeMC = " << typeMC << endl; + + // use for padding sigmabar vs. Theta from ON or OFF data + TString typeData = "ON"; + //TString typeData = "OFF"; + cout << "typeData = " << typeData << endl; //------------------------------------ @@ -389,9 +839,14 @@ // and "3D-ThetaPixSigma" // and "3D-ThetaPixDiff" - if (RHist) - { - cout << "Reading in file 'SigmaTheta.root' " << endl; - - TFile *infile = new TFile("SigmaTheta.root"); + + + TString outNameSigTh = "SigmaTheta_"; + outNameSigTh += typeData; + outNameSigTh += ".root"; + + + cout << "Reading in file " << outNameSigTh << endl; + + TFile *infile = new TFile(outNameSigTh); infile->ls(); @@ -400,5 +855,5 @@ if (!fHSigmaTheta) { - *fLog << "Object '2D-ThetaSigmabar' not found on root file" << endl; + cout << "Object '2D-ThetaSigmabar' not found on root file" << endl; return; } @@ -409,5 +864,5 @@ if (!fHSigmaPixTheta) { - *fLog << "Object '3D-ThetaPixSigma' not found on root file" << endl; + cout << "Object '3D-ThetaPixSigma' not found on root file" << endl; return; } @@ -418,24 +873,15 @@ if (!fHSigmaPixTheta) { - *fLog << "Object '3D-ThetaPixDiff' not found on root file" << endl; + cout << "Object '3D-ThetaPixDiff' not found on root file" << endl; return; } cout << "Object '3D-ThetaPixDiff' was read in" << endl; - } - else - { - MHSigmaTheta *sigtheta = (MHSigmaTheta*)pliston.FindObject("SigmaTheta"); - if (!sigtheta) - { - cout << "Object with name 'SigmaTheta' not found" << endl; - return; - } - - TH2D *fHSigmaTheta = sigtheta->GetSigmaTheta(); - TH3D *fHSigmaPixTheta = sigtheta->GetSigmaPixTheta(); - TH3D *fHDiffPixTheta = sigtheta->GetDiffPixTheta(); - } + //------------------------------------ + MTaskList tlist; + MParList plist; + + plist.AddToList(&tlist); @@ -539,5 +985,5 @@ MHillasSrcCalc csrc1("MSrcPosCam", fHilSrcName, "MHillas", "", fHilName); - MSelStandard selstand(fHillas, fHillasSrc); + MSelStandard selstand(fHilName, fHilSrcName); MFillH hfill1("MHHillas", fHilName); @@ -557,6 +1003,12 @@ const char* mtxName = "MatrixGammas"; - Int_t howMany = 30000; - char* outName = "matrix_gammas.root"; + Int_t howMany = 50000; + + TString outName = "matrix_gammas_"; + outName += typeMC; + outName += "_"; + outName += typeData; + outName += ".root"; + cout << "" << endl; @@ -581,18 +1033,16 @@ MHMatrix *pmatrix = new MHMatrix(mtxName); MHMatrix& matrix = *pmatrix; + + matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrix.AddColumn("MSigmabar.fSigmabar"); + matrix.AddColumn("log10(Hillas.fSize)"); + matrix.AddColumn("HillasSrc.fDist"); matrix.AddColumn("Hillas.fWidth"); matrix.AddColumn("Hillas.fLength"); - matrix.AddColumn("Hillas.fWidth*Hillas.fLength/Hillas.fSize"); - matrix.AddColumn("abs(Hillas.fAsym)"); + matrix.AddColumn("(log10(Hillas.fSize))/(Hillas.fWidth*Hillas.fLength)"); matrix.AddColumn("abs(Hillas.fM3Long)"); - matrix.AddColumn("abs(Hillas.fM3Trans)"); - matrix.AddColumn("abs(HillasSrc.fHeadTail)"); matrix.AddColumn("Hillas.fConc"); matrix.AddColumn("Hillas.fConc1"); - matrix.AddColumn("HillasSrc.fDist"); - matrix.AddColumn("log10(Hillas.fSize)"); - matrix.AddColumn("HillasSrc.fAlpha"); - matrix.AddColumn("cos(MMcEvt.fTelescopeTheta)"); - //matrix.AddColumn("MMcEvt.fTheta"); + plist.AddToList(pmatrix); @@ -600,8 +1050,24 @@ fillmatg.SetFilter(&flist); + + if (WMC1) + { + TString outNameImage = typeMC; + outNameImage += "_"; + outNameImage += typeData; + outNameImage += "1.root"; + MWriteRootFile write(outNameImage); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("Hillas", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("HillasSrc", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MMcRunHeader","RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + } + + //+++++++++++++++++++++++++++++++++++++++++++++++++++ - - MSelFinal selfinal(fHillas, fHillasSrc); //***************************** @@ -619,4 +1085,10 @@ tlist.AddToList(&selstand); + if (WMC1) + tlist.AddToList(&write); + + tlist.AddToList(&flist); + tlist.AddToList(&fillmatg); + tlist.AddToList(&hfill1); tlist.AddToList(&hfill2); @@ -624,8 +1096,5 @@ tlist.AddToList(&hfill4); - tlist.AddToList(&flist); - tlist.AddToList(&fillmatg); - - //tlist.AddToList(&selfinal); + //***************************** @@ -634,8 +1103,11 @@ // Execute event loop // + MProgressBar bar; MEvtLoop evtloop; evtloop.SetParList(&plist); + evtloop.SetProgressBar(&bar); Int_t maxevents = 1000000000; + //Int_t maxevents = 100; if ( !evtloop.Eventloop(maxevents) ) return; @@ -670,5 +1142,5 @@ // target distribution for the variable in column refcolumn; // set refcolumn negative if distribution is not to be changed - Int_t refcolumn = 12; + Int_t refcolumn = 1; Int_t nbins = 10; Float_t frombin = 0.5; @@ -711,9 +1183,859 @@ } - cout << "Macro AnalyseCT1 : End of section for MC gamma data" + cout << "Macro AnalyseCT1 : End of Job 3" << endl; cout << "=========================================================" << endl; } + + + + //--------------------------------------------------------------------- + // Job 4 + //====== + + // read OFF1 and MC1 data files + // + // - produce Neyman-Pearson plot + + if (Job4) + { + cout << "=====================================================" << endl; + cout << "Macro AnalyseCT1 : Start of Job 4" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job4 = " << Job4 << endl; + + + TString typeMC = "MC"; + cout << "typeMC = " << typeMC << endl; + + // use hadron matrix from ON or OFF data + TString typeData = "ON"; + //TString typeData = "OFF"; + cout << "typeData = " << typeData << endl; + + //************************************************************************* + // read in matrices of look-alike events + + const char* mtxName = "MatrixGammas"; + + + TString outName = "matrix_gammas_"; + outName += typeMC; + outName += "_"; + outName += typeData; + outName += ".root"; + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Get matrix for (gammas)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outName' + // + TFile *fileg = new TFile(outName); + + MHMatrix matrixg; + matrixg.Read(mtxName); + pmatrixg = &matrixg; + + delete fileg; + + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + TString outName = "matrix_hadrons_"; + outName += typeData; + outName += ".root"; + + cout << "" << endl; + cout << "========================================================" << endl; + cout << " Get matrix for (hadrons)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name mtxName + // + TFile *fileh = new TFile(outName); + + MHMatrix matrixh; + matrixh.Read(mtxName); + pmatrixh = &matrixh; + + delete fileh; + + //***************************************************************** + + MTaskList tliston; + MParList pliston; + pliston.AddToList(&tliston); + + pliston.AddToList(pmatrixg); + pliston.AddToList(pmatrixh); + + //:::::::::::::::::::::::::::::::::::::::::: + + MBinning binssize("BinningSize"); + binssize.SetEdgesLog(50, 10, 1.0e5); + pliston.AddToList(&binssize); + + MBinning binsdistc("BinningDist"); + binsdistc.SetEdges(50, 0, 1.4); + pliston.AddToList(&binsdistc); + + MBinning binswidth("BinningWidth"); + binswidth.SetEdges(50, 0, 1.0); + pliston.AddToList(&binswidth); + + MBinning binslength("BinningLength"); + binslength.SetEdges(50, 0, 1.0); + pliston.AddToList(&binslength); + + MBinning binsalpha("BinningAlpha"); + binsalpha.SetEdges(100, -100, 100); + pliston.AddToList(&binsalpha); + + MBinning binsasym("BinningAsym"); + binsasym.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsasym); + + MBinning binsht("BinningHeadTail"); + binsht.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsht); + + MBinning binsm3l("BinningM3Long"); + binsm3l.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3l); + + MBinning binsm3t("BinningM3Trans"); + binsm3t.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3t); + + + //..... + MBinning binsb("BinningSigmabar"); + binsb.SetEdges( 100, 0.0, 5.0); + pliston.AddToList(&binsb); + + MBinning binth("BinningTheta"); + binth.SetEdges( 70, -0.5, 69.5); + pliston.AddToList(&binth); + + MBinning binsdiff("BinningDiffsigma2"); + binsdiff.SetEdges(100, -5.0, 20.0); + pliston.AddToList(&binsdiff); + //:::::::::::::::::::::::::::::::::::::::::: + + + //------------------------------------------- + // create the tasks which should be executed + // + + TString filenameData = typeData; + filenameData += "1.root"; + + TString filenameMC = typeMC; + filenameMC += "_"; + filenameMC += typeData; + filenameMC += "1.root"; + + + cout << "filenameData = " << filenameData << endl; + cout << "filenameMC = " << filenameMC << endl; + + MReadMarsFile read("Events", filenameMC); + read.AddFile(filenameData); + read.DisableAutoScheme(); + + //....................................................................... + // g/h separation + + MMultiDimDistCalc ghsep; + ghsep.SetUseNumRows(25); + ghsep.SetUseKernelMethod(kFALSE); + + //....................................................................... + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "Hillas"; + TString fHilSrcName = "HillasSrc"; + + Float_t hadcut = 0.2; + MSelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetHadronnessCut(hadcut); + selfinalgh.SetAlphaCut(100.0); + + MFillH fillh("MHHadronness"); + + MSelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetHadronnessCut(hadcut); + selfinal.SetAlphaCut(20.0); + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + TString nxt = "HillasExt"; + MHHillasExt hhilext(nxt, "", fHilName); + pliston.AddToList(&hhilext); + TString namext = nxt; + namext += "[MHHillasExt]"; + MFillH hfill3(namext, fHilSrcName); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + + + //***************************** + // tasks to be executed + + tliston.AddToList(&read); + + tliston.AddToList(&ghsep); + tliston.AddToList(&fillh); + + tliston.AddToList(&selfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + + tliston.AddToList(&selfinal); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + TObject *c; + + c = pliston.FindObject("MHHadronness")->DrawClone(); + c->Print(); + + //c = pliston.FindObject("MHHillas")->DrawClone(); + + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + + //pliston.FindObject("MHHillasExt")->DrawClone(); + //c = pliston.FindObject(nxt)->DrawClone(); + + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + + cout << "Macro AnalyseCT1 : End of Job 4" << endl; + cout << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job 5 + //====== + + // - read in the matrices of look-alike events for gammas and hadrons + + // then read MC1 data file + // + // - perform the g/h separation + // - apply the final cuts + // + // - write a file of MC gamma events (MC2) + // (after the g/h separation and after the final cuts) + + if (Job5) + { + cout << "=====================================================" << endl; + cout << "Macro AnalyseCT1 : Start of Job 5" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job5, WMC2 = " + << Job5 << ", " << WMC2 << endl; + + TString typeInput = "MC"; + cout << "typeInput = " << typeInput << endl; + + TString typeMC = "MC"; + cout << "typeMC = " << typeMC << endl; + + // use hadron matrix from ON or OFF data + TString typeData = "ON"; + //TString typeData = "OFF"; + cout << "typeData = " << typeData << endl; + + //************************************************************************* + // read in matrices of look-alike events + + const char* mtxName = "MatrixGammas"; + + TString outName = "matrix_gammas_"; + outName += typeMC; + outName += "_"; + outName += typeData; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Get matrix for (gammas)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outName' + // + TFile *fileg = new TFile(outName); + + MHMatrix matrixg; + matrixg.Read(mtxName); + pmatrixg = &matrixg; + + delete fileg; + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + TString outName = "matrix_hadrons_"; + outName += typeData; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << " Get matrix for (hadrons)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name mtxName + // + TFile *fileh = new TFile(outName); + + MHMatrix matrixh; + matrixh.Read(mtxName); + pmatrixh = &matrixh; + + delete fileh; + + //************************************************************************* + + + MTaskList tliston; + MParList pliston; + pliston.AddToList(&tliston); + + pliston.AddToList(pmatrixg); + pliston.AddToList(pmatrixh); + + //:::::::::::::::::::::::::::::::::::::::::: + + MBinning binssize("BinningSize"); + binssize.SetEdgesLog(50, 10, 1.0e5); + pliston.AddToList(&binssize); + + MBinning binsdistc("BinningDist"); + binsdistc.SetEdges(50, 0, 1.4); + pliston.AddToList(&binsdistc); + + MBinning binswidth("BinningWidth"); + binswidth.SetEdges(50, 0, 1.0); + pliston.AddToList(&binswidth); + + MBinning binslength("BinningLength"); + binslength.SetEdges(50, 0, 1.0); + pliston.AddToList(&binslength); + + MBinning binsalpha("BinningAlpha"); + binsalpha.SetEdges(100, -100, 100); + pliston.AddToList(&binsalpha); + + MBinning binsasym("BinningAsym"); + binsasym.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsasym); + + MBinning binsht("BinningHeadTail"); + binsht.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsht); + + MBinning binsm3l("BinningM3Long"); + binsm3l.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3l); + + MBinning binsm3t("BinningM3Trans"); + binsm3t.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3t); + + + //..... + MBinning binsb("BinningSigmabar"); + binsb.SetEdges( 100, 0.0, 5.0); + pliston.AddToList(&binsb); + + MBinning binth("BinningTheta"); + binth.SetEdges( 70, -0.5, 69.5); + pliston.AddToList(&binth); + + MBinning binsdiff("BinningDiffsigma2"); + binsdiff.SetEdges(100, -5.0, 20.0); + pliston.AddToList(&binsdiff); + //:::::::::::::::::::::::::::::::::::::::::: + + + //------------------------------------------- + // create the tasks which should be executed + // + + TString filenameMC = typeInput; + filenameMC += "_"; + filenameMC += typeData; + filenameMC += "1.root"; + + readname = "ReadCT1MCdata"; + MReadMarsFile read("Events", filenameMC); + read.DisableAutoScheme(); + + + //....................................................................... + // g/h separation + + MMultiDimDistCalc ghsep; + ghsep.SetUseNumRows(25); + ghsep.SetUseKernelMethod(kFALSE); + + + + //....................................................................... + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "Hillas"; + TString fHilSrcName = "HillasSrc"; + + Float_t hadcut = 0.2; + MSelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetHadronnessCut(hadcut); + selfinalgh.SetAlphaCut(100.0); + + MFillH fillh("MHHadronness"); + + MSelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetHadronnessCut(hadcut); + selfinal.SetAlphaCut(20.0); + + if (WMC2) + { + TString outNameImage = typeInput; + outNameImage += "_"; + outNameImage += typeData; + outNameImage += "2.root"; + MWriteRootFile write(outNameImage); + write.AddContainer("MHadronness", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("Hillas", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("HillasSrc", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MMcRunHeader","RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + } + + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + TString nxt = "HillasExt"; + MHHillasExt hhilext(nxt, "", fHilName); + pliston.AddToList(&hhilext); + TString namext = nxt; + namext += "[MHHillasExt]"; + MFillH hfill3(namext, fHilSrcName); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + + + + //***************************** + // tasks to be executed + + tliston.AddToList(&read); + + tliston.AddToList(&ghsep); + tliston.AddToList(&fillh); + + tliston.AddToList(&selfinalgh); + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + + tliston.AddToList(&selfinal); + if (WMC2) + tliston.AddToList(&write); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + TObject *c; + + c = pliston.FindObject("MHHadronness")->DrawClone(); + c->Print(); + + c = pliston.FindObject("MHHillas")->DrawClone(); + + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + + //pliston.FindObject("MHHillasExt")->DrawClone(); + c = pliston.FindObject(nxt)->DrawClone(); + + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + + cout << "Macro AnalyseCT1 : End of Job 5" << endl; + cout << "===================================================" << endl; + } + + + //--------------------------------------------------------------------- + // Job 6 + //====== + + // - read in the matrices of look-alike events for gammas and hadrons + + // then read ON1 data file + // + // - perform the g/h separation + // - apply the final cuts + // + // - write a file of ON events (ON2) + // (after the g/h separation and after the final cuts) + // (to be used for the flux calculation) + // + // - do the flux calculation + + if (Job6) + { + cout << "=====================================================" << endl; + cout << "Macro AnalyseCT1 : Start of Job 6" << endl; + + cout << "" << endl; + cout << "Macro AnalyseCT1 : Job6, WON2 = " + << Job6 << ", " << WON2 << endl; + + TString typeInput = "ON"; + cout << "typeInput = " << typeInput << endl; + + TString typeMC = "MC"; + cout << "typeMC = " << typeMC << endl; + + // use hadron matrix from ON or OFF data + TString typeData = "ON"; + //TString typeData = "OFF"; + cout << "typeData = " << typeData << endl; + + + //************************************************************************* + // read in matrices of look-alike events + + const char* mtxName = "MatrixGammas"; + + TString outName = "matrix_gammas_"; + outName += typeMC; + outName += "_"; + outName += typeData; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << "Get matrix for (gammas)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name 'mtxName' from file 'outName' + // + TFile *fileg = new TFile(outName); + + MHMatrix matrixg; + matrixg.Read(mtxName); + pmatrixg = &matrixg; + + delete fileg; + + //***************************************************************** + + const char* mtxName = "MatrixHadrons"; + + TString outName = "matrix_hadrons_"; + outName += typeData; + outName += ".root"; + + + cout << "" << endl; + cout << "========================================================" << endl; + cout << " Get matrix for (hadrons)" << endl; + cout << "matrix name = " << mtxName << endl; + cout << "name of root file = " << outName << endl; + cout << "" << endl; + + + // read in the object with the name mtxName + // + TFile *fileh = new TFile(outName); + + MHMatrix matrixh; + matrixh.Read(mtxName); + pmatrixh = &matrixh; + + delete fileh; + + //************************************************************************* + + + MTaskList tliston; + MParList pliston; + pliston.AddToList(&tliston); + + pliston.AddToList(pmatrixg); + pliston.AddToList(pmatrixh); + + //:::::::::::::::::::::::::::::::::::::::::: + + MBinning binssize("BinningSize"); + binssize.SetEdgesLog(50, 10, 1.0e5); + pliston.AddToList(&binssize); + + MBinning binsdistc("BinningDist"); + binsdistc.SetEdges(50, 0, 1.4); + pliston.AddToList(&binsdistc); + + MBinning binswidth("BinningWidth"); + binswidth.SetEdges(50, 0, 1.0); + pliston.AddToList(&binswidth); + + MBinning binslength("BinningLength"); + binslength.SetEdges(50, 0, 1.0); + pliston.AddToList(&binslength); + + MBinning binsalpha("BinningAlpha"); + binsalpha.SetEdges(100, -100, 100); + pliston.AddToList(&binsalpha); + + MBinning binsasym("BinningAsym"); + binsasym.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsasym); + + MBinning binsht("BinningHeadTail"); + binsht.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsht); + + MBinning binsm3l("BinningM3Long"); + binsm3l.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3l); + + MBinning binsm3t("BinningM3Trans"); + binsm3t.SetEdges(50, -1.5, 1.5); + pliston.AddToList(&binsm3t); + + + //..... + MBinning binsb("BinningSigmabar"); + binsb.SetEdges( 100, 0.0, 5.0); + pliston.AddToList(&binsb); + + MBinning binth("BinningTheta"); + binth.SetEdges( 70, -0.5, 69.5); + pliston.AddToList(&binth); + + MBinning binsdiff("BinningDiffsigma2"); + binsdiff.SetEdges(100, -5.0, 20.0); + pliston.AddToList(&binsdiff); + //:::::::::::::::::::::::::::::::::::::::::: + + + //------------------------------------------- + // create the tasks which should be executed + // + + TString filenameData = typeInput; + filenameData += "1.root"; + + readname = "ReadCT1ONdata"; + MReadMarsFile read("Events", filenameData); + read.DisableAutoScheme(); + + + //....................................................................... + // g/h separation + + MMultiDimDistCalc ghsep; + ghsep.SetUseNumRows(25); + ghsep.SetUseKernelMethod(kFALSE); + + + + //....................................................................... + + // geometry is needed in MHHillas... classes + MGeomCam *fGeom = + (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam"); + + TString fHilName = "Hillas"; + TString fHilSrcName = "HillasSrc"; + + Float_t hadcut = 0.2; + MSelFinal selfinalgh(fHilName, fHilSrcName); + selfinalgh.SetHadronnessCut(hadcut); + selfinalgh.SetAlphaCut(100.0); + + MFillH fillh("MHHadronness"); + + MSelFinal selfinal(fHilName, fHilSrcName); + selfinal.SetHadronnessCut(hadcut); + selfinal.SetAlphaCut(20.0); + + if (WON2) + { + TString outNameImage = typeInput; + outNameImage += "_"; + outNameImage += typeData; + outNameImage += "2.root"; + MWriteRootFile write(outNameImage); + write.AddContainer("MHadronness", "Events"); + write.AddContainer("MSigmabar", "Events"); + write.AddContainer("Hillas", "Events"); + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("HillasSrc", "Events"); + write.AddContainer("MRawRunHeader", "RunHeaders"); + //write.AddContainer("MMcRunHeader","RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + } + + + MFillH hfill1("MHHillas", fHilName); + MFillH hfill2("MHStarMap", fHilName); + + TString nxt = "HillasExt"; + MHHillasExt hhilext(nxt, "", fHilName); + pliston.AddToList(&hhilext); + TString namext = nxt; + namext += "[MHHillasExt]"; + MFillH hfill3(namext, fHilSrcName); + + MFillH hfill4("MHHillasSrc", fHilSrcName); + + + + //***************************** + // tasks to be executed + + tliston.AddToList(&read); + + tliston.AddToList(&ghsep); + tliston.AddToList(&fillh); + + tliston.AddToList(&selfinalgh); + + tliston.AddToList(&hfill1); + tliston.AddToList(&hfill2); + tliston.AddToList(&hfill3); + tliston.AddToList(&hfill4); + + tliston.AddToList(&selfinal); + if (WON2) + tliston.AddToList(&write); + + + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&pliston); + evtloop.SetProgressBar(&bar); + + Int_t maxevents = 1000000000; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; + + tliston.PrintStatistics(0, kTRUE); + + + //------------------------------------------- + // Display the histograms + // + TObject *c; + + c = pliston.FindObject("MHHadronness")->DrawClone(); + c->Print(); + + //c = pliston.FindObject("MHHillas")->DrawClone(); + + c = pliston.FindObject("MHHillasSrc")->DrawClone(); + + ////pliston.FindObject("MHHillasExt")->DrawClone(); + //c = pliston.FindObject(nxt)->DrawClone(); + + c = pliston.FindObject("MHStarMap")->DrawClone(); + + + + cout << "Macro AnalyseCT1 : End of Job 6" << endl; + cout << "=======================================================" << endl; + } + //--------------------------------------------------------------------- } @@ -721,2 +2043,8 @@ + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MMultiDimDistCalc.cc (revision 1809) @@ -208,6 +208,18 @@ } - //fHadronness->SetHadronness(dg/(dg+dh)); - fHadronness->SetHadronness(exp(-dh/dg)); + Double_t arg; + + if (dg+dh != 0.0) + arg = dg / (dg+dh); + else + arg = 1.e10; + //fHadronness->SetHadronness(arg); + + if (dg != 0.0) + arg = exp(-dh/dg); + else + arg = 0.0; + fHadronness->SetHadronness(arg); + return kTRUE; Index: trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc (revision 1809) @@ -1,3 +1,2 @@ - /* ======================================================================== *\ ! @@ -102,4 +101,8 @@ fHDiffPixTheta = fHist3Diff; + fHSigmaTheta->SetDirectory(NULL); + fHSigmaPixTheta->SetDirectory(NULL); + fHDiffPixTheta->SetDirectory(NULL); + Print(); } @@ -239,4 +242,8 @@ fHDiffPixTh->SetZTitle("Sigma^2-Sigmabar^2"); + //-------------------------------------------------------------------- + + memset(fErrors, 0, sizeof(fErrors)); + return kTRUE; } @@ -251,17 +258,26 @@ //*fLog << "Entry MPadSchweizer::Process();" << endl; + Int_t rc; + const UInt_t npix = fEvt->GetNumPixels(); + + Double_t SigbarOld; + + //*fLog << "before padding : " << endl; + //SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); + //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); - } + //for (UInt_t i=0; ioperator[](i); + // Int_t j = pix.GetPixId(); + + // MPedestalPix &ppix = fPed->operator[](j); + // ppix.SetMeanRms(0.0); + //} //$$$$$$$$$$$$$$$$$$$$$$$$$$ @@ -269,14 +285,17 @@ // Calculate average sigma of the event // - Double_t SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); + SigbarOld = fSigmabar->Calc(*fCam, *fPed, *fEvt); //fSigmabar->Print(""); - //if (SigbarOld > 0.0) - //{ + if (SigbarOld > 0.0) + { //*fLog << "MPadSchweizer::Process(); Sigmabar of event to be padded is > 0 : " // << SigbarOld << ". Stop event loop " << endl; // input data should have Sigmabar = 0; stop event loop - // return kFALSE; - //} + + rc = 1; + fErrors[rc]++; + return kCONTINUE; + } Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); @@ -296,4 +315,7 @@ << Theta << ", " << binNumber << "; Skip event " << endl; // event cannot be padded; skip event + + rc = 2; + fErrors[rc]++; return kCONTINUE; } @@ -308,4 +330,7 @@ // event cannot be padded; skip event delete fHSigma; + + rc = 3; + fErrors[rc]++; return kCONTINUE; } @@ -331,4 +356,7 @@ *fLog << "MPadSchweizer::Process(); target Sigmabar is less than SigbarOld : " << Sigmabar << ", " << SigbarOld << ", Skip event" << endl; + + rc = 4; + fErrors[rc]++; return kCONTINUE; } @@ -426,5 +454,8 @@ << binTheta << " and pixel bin " << binPixel << " has no entries; aborting " << endl; - return kERROR; + + rc = 5; + fErrors[rc]++; + return kCONTINUE; } @@ -465,5 +496,8 @@ << binTheta << " and pixel bin " << binPixel << " has no entries; aborting " << endl; - return kERROR; + + rc = 6; + fErrors[rc]++; + return kCONTINUE; } @@ -553,5 +587,7 @@ // Calculate Sigmabar again and crosscheck + Double_t SigbarNew = fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << "after padding : " << endl; //fSigmabar->Print(""); @@ -585,4 +621,7 @@ //*fLog << "Exit MPadSchweizer::Process();" << endl; + rc = 0; + fErrors[rc]++; + return kTRUE; @@ -594,4 +633,42 @@ Bool_t MPadSchweizer::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 << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts survived the padding!" << endl; + *fLog << endl; + + } + + //--------------------------------------------------------------- TCanvas &c = *(MH::MakeDefCanvas("PadSchweizer", "", 900, 1200)); c.Divide(3, 4); @@ -600,19 +677,22 @@ c.cd(1); + fHSigmaTheta->SetDirectory(NULL); fHSigmaTheta->SetTitle("2D : Sigmabar, \\Theta (reference sample)"); - fHSigmaTheta->DrawClone(); + fHSigmaTheta->DrawCopy(); fHSigmaTheta->SetBit(kCanDelete); //c.cd(1); - //fHSigmaPixTheta->DrawClone(); + //fHSigmaPixTheta->DrawCopy(); //fHSigmaPixTheta->SetBit(kCanDelete); c.cd(4); - fHSigbarTheta->DrawClone(); + fHSigbarTheta->SetDirectory(NULL); + fHSigbarTheta->DrawCopy(); fHSigbarTheta->SetBit(kCanDelete); c.cd(7); - fHNSB->DrawClone(); + fHNSB->SetDirectory(NULL); + fHNSB->DrawCopy(); fHNSB->SetBit(kCanDelete); @@ -625,21 +705,20 @@ c.cd(2); l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zx"); - + l->SetDirectory(NULL); l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); l->SetXTitle("\\Theta [\\circ]"); l->SetYTitle("Sigma^2-Sigmabar^2"); - *fLog << "before box" << endl; - - l->Draw("box"); + l->DrawCopy("box"); l->SetBit(kCanDelete);; c.cd(5); l = (TH2D*) ((TH3*)fHDiffPixTh)->Project3D("zy"); + l->SetDirectory(NULL); l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); l->SetXTitle("pixel"); l->SetYTitle("Sigma^2-Sigmabar^2"); - l->Draw("box"); + l->DrawCopy("box"); l->SetBit(kCanDelete);; @@ -655,18 +734,20 @@ c.cd(3); k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zx"); + k->SetDirectory(NULL); k->SetTitle("Sigma vs. \\Theta (all pixels)"); k->SetXTitle("\\Theta [\\circ]"); k->SetYTitle("Sigma"); - k->Draw("box"); + k->DrawCopy("box"); k->SetBit(kCanDelete); c.cd(6); k = (TH2D*) ((TH3*)fHSigPixTh)->Project3D("zy"); + k->SetDirectory(NULL); k->SetTitle("Sigma vs. pixel number (all \\Theta)"); k->SetXTitle("pixel"); k->SetYTitle("Sigma"); - k->Draw("box"); + k->DrawCopy("box"); k->SetBit(kCanDelete);; @@ -677,9 +758,11 @@ c.cd(10); - fHSigmaPedestal->DrawClone(); + fHSigmaPedestal->SetDirectory(NULL); + fHSigmaPedestal->DrawCopy(); fHSigmaPedestal->SetBit(kCanDelete); c.cd(11); - fHPhotons->DrawClone(); + fHPhotons->SetDirectory(NULL); + fHPhotons->DrawCopy(); fHPhotons->SetBit(kCanDelete); Index: trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MPadSchweizer.h (revision 1809) @@ -39,4 +39,6 @@ Int_t fRunType; Int_t fGroup; + + Int_t fErrors[7]; // plots used for the padding Index: trunk/MagicSoft/Mars/manalysis/MSelBasic.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelBasic.cc (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelBasic.cc (revision 1809) @@ -44,4 +44,5 @@ #include "MGeomCam.h" #include "MPedestalCam.h" +#include "MPedestalPix.h" #include "MGeomPix.h" @@ -119,7 +120,83 @@ Bool_t MSelBasic::Process() { + /* + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + *fLog << "=========================================================" << endl; + *fLog << "" << endl; + *fLog << "MmcEvt data : " << endl; + *fLog << "" << endl; + fMcEvt->Print(); + *fLog << "" << endl; + *fLog << "PartId() = " << fMcEvt->GetPartId() << endl; + *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl; + *fLog << "Theta() = " << fMcEvt->GetTheta() << endl; + + *fLog << "Phi() = " << fMcEvt->GetPhi() << endl; + //*fLog << "CoreD() = " << fMcEvt->GetCoreD() << endl; + //*fLog << "CoreX() = " << fMcEvt->GetCoreX() << endl; + + //*fLog << "CoreY() = " << fMcEvt->GetCoreY() << endl; + *fLog << "Impact() = " << fMcEvt->GetImpact() << endl; + //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl; + + //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl; + //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl; + //*fLog << "PassPhotCone() = " << fMcEvt->GetPassPhotCone() << endl; + + //*fLog << "PhotElfromShower() = " << fMcEvt->GetPhotElfromShower() << endl; + //*fLog << "PhotElinCamera() = " << fMcEvt->GetPhotElinCamera() << endl; + *fLog << "TelescopePhi() = " << fMcEvt->GetTelescopePhi() << endl; + + *fLog << "TelescopeTheta() = " << fMcEvt->GetTelescopeTheta() << endl; + *fLog << "Impact() = " << fMcEvt->GetImpact() << endl; + //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl; + *fLog << "" << endl; + *fLog << "=========================================================" << endl; + + *fLog << "=========================================================" << endl; + *fLog << "" << endl; + *fLog << "MPedestalPix data : " << endl; + *fLog << "" << endl; + + Int_t ntot; + ntot = fPed->GetSize(); + *fLog << "MeanRms() :" << endl; + for (Int_t i=0; iGetNumPixels(); + *fLog << "MCerPhotPix : pix.GetNumPhotons()" << endl; + for (Int_t i=0; iGetNumPixels(); + *fLog << "MCerPhotPix : pix.GetErrorPhot()" << endl; + for (Int_t i=0; iGetRunNumber() < 1025 ) + //if ( fRawRun->GetRunNumber() < 16279 ) //{ // rc = 1; @@ -128,8 +205,21 @@ Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - if (Theta > 45.0 || !SwTrigger() ) - { - //*fLog << "MSelBasic::Process; Theta = " << Theta << endl; + if ( Theta < 0.0 ) + { + *fLog << "MSelBasic::Process; Run, Event, Theta = " + << fRawRun->GetRunNumber()<< ", " + << fMcEvt->GetEvtNumber() << ", " << Theta << endl; rc = 1; + } + + else if ( Theta > 45.0 ) + { + rc = 2; + } + + else if ( !SwTrigger() ) + { + //*fLog << "MSelBasic::Process; SwTrigger = " << SwTrigger << endl; + rc = 3; } @@ -223,5 +313,9 @@ *fLog << GetDescriptor() << " execution statistics:" << endl; *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Basic selections are not fullfilled" << endl; + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle < 0" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Zenith angle too high" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Software trigger not fullfilled" << endl; *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Basic selections!" << endl; @@ -230,2 +324,3 @@ return kTRUE; } + Index: trunk/MagicSoft/Mars/manalysis/MSelBasic.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelBasic.h (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelBasic.h (revision 1809) @@ -30,5 +30,5 @@ const MRawRunHeader *fRawRun; - Int_t fErrors[2]; + Int_t fErrors[4]; public: Index: trunk/MagicSoft/Mars/manalysis/MSelFinal.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelFinal.cc (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelFinal.cc (revision 1809) @@ -41,4 +41,5 @@ #include "MHillas.h" +#include "MHillasExt.h" #include "MHillasSrc.h" #include "MCerPhotEvt.h" @@ -57,5 +58,5 @@ // Default constructor. // -MSelFinal::MSelFinal(MHillas *parhil, MHillasSrc *parhilsrc, +MSelFinal::MSelFinal(const char *HilName, const char *HilSrcName, const char *name, const char *title) { @@ -63,6 +64,9 @@ fTitle = title ? title : "Task to evaluate the Final Cuts"; - fHil = parhil; - fHilsrc = parhilsrc; + fHilName = HilName; + fHilSrcName = HilSrcName; + + fHadronnessCut = 0.2; + fAlphaCut = 100.0; //degrees } @@ -75,4 +79,19 @@ Bool_t MSelFinal::PreProcess(MParList *pList) { + fHil = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt"); + if (!fHil) + { + *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; + return kFALSE; + } + + fHadronness = (MHadronness*)pList->FindObject("MHadronness"); if (!fHadronness) @@ -90,22 +109,4 @@ } - fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt"); - if (!fEvt) - { - *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - - fCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fCam) - { - *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; - return kFALSE; - } - fMm2Deg = fCam->GetConvMm2Deg(); - - *fLog << "fMm2Deg = " << fMm2Deg << endl; - memset(fErrors, 0, sizeof(fErrors)); @@ -122,16 +123,26 @@ Bool_t MSelFinal::Process() { + //*fLog << "Entry MSelFinal; fHilSrc = " << fHilSrc << endl; + + + Int_t rc = 0; - Double_t alphacut = 20.0; + Double_t modalpha = fabs( fHilSrc->GetAlpha() ); - Double_t modalpha = fabs( fHilsrc->GetAlpha() ); Double_t h = fHadronness->GetHadronness(); - if ( h>0.5 || modalpha > alphacut ) + if ( h>fHadronnessCut ) { //*fLog << "MSelFinal::Process; h, alpha = " << h << ", " - // << fHilsrc->GetAlpha() << endl; + // << fHilSrc->GetAlpha() << endl; rc = 1; + } + + else if ( modalpha > fAlphaCut ) + { + //*fLog << "MSelFinal::Process; h, alpha = " << h << ", " + // << fHilSrc->GetAlpha() << endl; + rc = 2; } @@ -153,9 +164,22 @@ *fLog << GetDescriptor() << " execution statistics:" << endl; *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Final selections are not fullfilled" << endl; + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) + << (int)(fErrors[1]*100/GetNumExecutions()) + << "%) Evts skipped due to: g/h separation cut (" << fHadronnessCut + << ")" << endl; - *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Final selections!" << endl; + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) + << (int)(fErrors[2]*100/GetNumExecutions()) + << "%) Evts skipped due to: cut in ALPHA (" << fAlphaCut + << " degrees)" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Evts survived Final selections!" << endl; *fLog << endl; return kTRUE; } + + + Index: trunk/MagicSoft/Mars/manalysis/MSelFinal.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelFinal.h (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelFinal.h (revision 1809) @@ -28,12 +28,17 @@ MMcEvt *fMcEvt; MHillas *fHil; - MHillasSrc *fHilsrc; + MHillasSrc *fHilSrc; MHadronness *fHadronness; Double_t fMm2Deg; // conversion mm to degrees in camera - Int_t fErrors[2]; + Int_t fErrors[3]; + TString fHilName; + TString fHilSrcName; + + Float_t fHadronnessCut; + Float_t fAlphaCut; public: - MSelFinal(MHillas *fHil, MHillasSrc *fHilsrc, + MSelFinal(const char *HilName, const char *HilSrcName, const char *name=NULL, const char *title=NULL); @@ -42,4 +47,6 @@ Bool_t PostProcess(); + void SetHadronnessCut(Float_t hadcut) { fHadronnessCut = hadcut; } + void SetAlphaCut(Float_t alpha) { fAlphaCut = alpha; } ClassDef(MSelFinal, 0) // Task to evaluate final cuts Index: trunk/MagicSoft/Mars/manalysis/MSelStandard.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelStandard.cc (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelStandard.cc (revision 1809) @@ -1,2 +1,3 @@ + /* ======================================================================== *\ ! @@ -39,4 +40,5 @@ #include "MHillas.h" +#include "MHillasExt.h" #include "MHillasSrc.h" #include "MCerPhotEvt.h" @@ -54,5 +56,5 @@ // Default constructor. // -MSelStandard::MSelStandard(const MHillas *parhil, const MHillasSrc *parhilsrc, +MSelStandard::MSelStandard(const char *HilName, const char *HilSrcName, const char *name, const char *title) { @@ -60,6 +62,6 @@ fTitle = title ? title : "Task to evaluate the Standard Cuts"; - fHil = parhil; - fHilsrc = parhilsrc; + fHilName = HilName; + fHilSrcName = HilSrcName; } @@ -72,4 +74,19 @@ Bool_t MSelStandard::PreProcess(MParList *pList) { + fHil = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt"); + if (!fHil) + { + *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl; + return kFALSE; + } + + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); + if (!fHilSrc) + { + *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); if (!fMcEvt) @@ -95,5 +112,5 @@ fMm2Deg = fCam->GetConvMm2Deg(); - *fLog << "fMm2Deg = " << fMm2Deg << endl; + //*fLog << "fMm2Deg = " << fMm2Deg << endl; memset(fErrors, 0, sizeof(fErrors)); @@ -114,7 +131,7 @@ Int_t rc = 0; - //Double_t fLength = fHil->GetLength() * fMm2Deg; - //Double_t fWidth = fHil->GetWidth() * fMm2Deg; - Double_t fDist = fHilsrc->GetDist()* fMm2Deg; + Double_t fLength = fHil->GetLength() * fMm2Deg; + Double_t fWidth = fHil->GetWidth() * fMm2Deg; + Double_t fDist = fHilSrc->GetDist()* fMm2Deg; //Double_t fDelta = fHil->GetDelta() * kRad2Deg; Double_t fSize = fHil->GetSize(); @@ -122,6 +139,5 @@ Int_t fNumCorePixels = fHil->GetNumCorePixels(); - if ( fSize <= 60.0 || fDist< 0.4 || fDist > 1.1 - || fNumUsedPixels >= 92 || fNumCorePixels < 4) + if ( fNumUsedPixels >= 92 || fNumCorePixels < 4 ) { //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = " @@ -129,4 +145,19 @@ // << fNumCorePixels << endl; rc = 1; + } + + else if ( fSize <= 60.0 || fDist< 0.4 || fDist > 1.1 ) + { + //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = " + // << fSize << ", " << fDist << ", " << fNumUsedPixels << ", " + // << fNumCorePixels << endl; + rc = 2; + } + + else if ( fLength <= 0.0 || fWidth <= 0.0 ) + { + //*fLog << "MSelStandard::Process; fLength, fWidth = " + // << fLength << ", " << fWidth << endl; + rc = 3; } @@ -149,5 +180,9 @@ *fLog << GetDescriptor() << " execution statistics:" << endl; *fLog << dec << setfill(' '); - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Standard selections are not fullfilled" << endl; + *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) << (int)(fErrors[1]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on no.of used or core pxels not fullfilled" << endl; + + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) << (int)(fErrors[2]*100/GetNumExecutions()) << "%) Evts skipped due to: Requirements on SIZE or DIST not fullfilled" << endl; + + *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) << (int)(fErrors[3]*100/GetNumExecutions()) << "%) Evts skipped due to: Length or Width is <= 0" << endl; *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) << "%) Evts survived Standard selections!" << endl; @@ -156,2 +191,4 @@ return kTRUE; } + + Index: trunk/MagicSoft/Mars/manalysis/MSelStandard.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSelStandard.h (revision 1808) +++ trunk/MagicSoft/Mars/manalysis/MSelStandard.h (revision 1809) @@ -23,15 +23,17 @@ { private: - const MGeomCam *fCam; // Camera Geometry - const MCerPhotEvt *fEvt; // Cerenkov Photon Event - const MMcEvt *fMcEvt; - const MHillas *fHil; - const MHillasSrc *fHilsrc; + MGeomCam *fCam; // Camera Geometry + MCerPhotEvt *fEvt; // Cerenkov Photon Event + MMcEvt *fMcEvt; + MHillas *fHil; + MHillasSrc *fHilSrc; - Double_t fMm2Deg; // conversion mm to degrees in camera - Int_t fErrors[2]; + Double_t fMm2Deg; // conversion mm to degrees in camera + Int_t fErrors[4]; + TString fHilName; + TString fHilSrcName; public: - MSelStandard(const MHillas *fHil, const MHillasSrc *fHilsrc, + MSelStandard(const char *HilName, const char *HilSrcName, const char *name=NULL, const char *title=NULL); Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 1808) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 1809) @@ -122,4 +122,5 @@ { const char *name = gSystem->ExpandPathName(txt); + TString fname(name); delete [] name; @@ -336,6 +337,15 @@ // float frms_pedsig_phot[iMAXNUMPIX]; // standard deviation of the calibrated signals from the pedestal run */ fPedest->InitSize(iMAXNUMPIX); + *fLog << "PedestalRMS : "; for (Int_t i=0; iSetReadyToSave(); @@ -373,10 +383,24 @@ * be treated further. */ + *fLog << "outpars.bmontecarlo = " << outpars.bmontecarlo << endl; + *fLog << "outpars.imcparticle = " << outpars.imcparticle << endl; + *fLog << "outpars.dsourcera_hours = " << outpars.dsourcera_hours << endl; + *fLog << "outpars.dsourcedec_deg = " << outpars.dsourcedec_deg << endl; + + //*fLog << "File is a "; + //*fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data"); + //*fLog << " file." << endl; + + + // Next statement commented out because bmontecarlo was set wrongly + //fIsMcFile = outpars.bmontecarlo==TRUE; + fIsMcFile = (outpars.dsourcera_hours == 0.0 && + outpars.dsourcedec_deg == 0.0 && + outpars.imcparticle != 0 ); + *fLog << "File is a "; - *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data"); + *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data"); *fLog << " file." << endl; - - fIsMcFile = outpars.bmontecarlo==TRUE; - + *fLog << " " << endl; } @@ -433,6 +457,8 @@ ProcessRunHeader(outpars); + //rwagner: ReInit whenever new run commences // rc==-1 means: ReInit didn't work out + if (!fTaskList->ReInit(fParList)) return -1; @@ -453,4 +479,5 @@ TString m = cEND_EVENTS_TEMPLATE; Int_t p = m.First('%'); + if (!s.BeginsWith(m(0,p))) @@ -515,5 +542,5 @@ // - // Check for the existance of a next file to read + // Check for the existence of a next file to read // TNamed *file = (TNamed*)fFileNames->First(); @@ -538,9 +565,13 @@ if (!CheckHeader(fname)) - return kFALSE; + { + *fLog << "OpenNextFile : CheckHeader(fname) is FALSE" << endl; + return kFALSE; + } fIn = new ifstream(fname); *fLog << inf << "-----------------------------------------------------------------------" << endl; + switch (ReadRunHeader()) @@ -553,6 +584,11 @@ return kFALSE; default: + *fLog << "OpenNextFile : after ReadRunHeader, FnumEventsInRun = " + << fNumEventsInRun << endl; return kTRUE; } + + + } @@ -752,4 +788,32 @@ void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event) { + if (fRawRunHeader->GetRunNumber() == 1) + { + *fLog << "eventrecord" << endl; + *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl; + *fLog << "isecfrac_200ns = " << event.isecfrac_200ns << endl; + *fLog << "snot_ok_flags = " << event.snot_ok_flags << endl; + *fLog << "ialt_arcs = " << event.ialt_arcs << endl; + *fLog << "iaz_arcs = " << event.iaz_arcs << endl; + *fLog << "ipreproc_alt_arcs = " << event.ipreproc_alt_arcs << endl; + *fLog << "ipreproc_az_arcs = " << event.ipreproc_az_arcs << endl; + *fLog << "ifieldrot_arcs = " << event.ifieldrot_arcs << endl; + + *fLog << "srate_millihz = " << event.srate_millihz << endl; + *fLog << "fhourangle = " << event.fhourangle << endl; + *fLog << "fmcenergy_tev = " << event.fmcenergy_tev << endl; + *fLog << "fmcsize_phel = " << event.fmcsize_phel << endl; + *fLog << "imcimpact_m = " << event.imcimpact_m << endl; + *fLog << "imcparticle = " << event.imcparticle << endl; + *fLog << "imctriggerflag = " << event.imctriggerflag << endl; + } + + + for (Int_t i=0; iSetReadyToSave(); @@ -801,5 +870,5 @@ // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds) - fMcEvt->Fill(0, /*fEvtNum*/ + fMcEvt->Fill(event.isecs_since_midday, //0, /*fEvtNum*/ fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/ fIsMcFile ? event.fmcenergy_tev*1000 : 0, @@ -867,5 +936,11 @@ // an event // - if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0]) + + // "goto TONI" was introduced in order to check for a footer record + // after reading a run header; this is necessary for runs with + // zero events after the filter + TONI: + + if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0]) return kTRUE; @@ -874,4 +949,5 @@ // must be an event // + switch (ReadRunFooter()) { @@ -906,4 +982,5 @@ *fLog << "-----------------------------------------------------------------------" << endl; + if (ReadRunHeader() < 0) { @@ -911,4 +988,7 @@ return kFALSE; } + + goto TONI; + return kTRUE; } Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 1808) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 1809) @@ -48,4 +48,7 @@ UInt_t fNumFilterEvts; // number of events mentioned in the runs footers + Float_t savePedRMS[127]; + + Bool_t OpenNextFile(); Index: trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1808) +++ trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1809) @@ -118,4 +118,6 @@ Bool_t MFEventSelector::PreProcess(MParList *plist) { + memset(fErrors, 0, sizeof(fErrors)); + fNumSelectedEvts = 0; if (fSelRatio>0) @@ -153,4 +155,6 @@ Bool_t MFEventSelector::Process() { + Int_t rc; + const Float_t evt = gRandom->Uniform(); @@ -161,6 +165,15 @@ if (fResult) + { fNumSelectedEvts++; + rc = 0; + fErrors[rc]++; + return kTRUE; + } + + rc = 1; + fErrors[rc]++; + return kTRUE; } @@ -172,4 +185,21 @@ Bool_t MFEventSelector::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()) + << "%) Events not selected" << endl; + + *fLog << " " << fErrors[0] << " (" + << (int)(fErrors[0]*100/GetNumExecutions()) + << "%) Events selected!" << endl; + *fLog << endl; + } + + //--------------------------------- if (TestBit(kNumTotalFromFile)) fNumTotalEvts = -1; @@ -199,3 +229,2 @@ } - Index: trunk/MagicSoft/Mars/mfilter/MFEventSelector.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1808) +++ trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1809) @@ -37,4 +37,6 @@ */ + Int_t fErrors[2]; + public: // MFEventSelector(); Index: trunk/MagicSoft/Mars/mhist/MHHadronness.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 1808) +++ trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 1809) @@ -235,5 +235,11 @@ } - return val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5); + Float_t retValue; + if (val2x-val1x != 0.0) + retValue = val1y - (val2y-val1y)/(val2x-val1x) * (val1x-0.5); + else + retValue = 0.0; + + return retValue; } @@ -252,6 +258,41 @@ fQfac->Set(n); - const Stat_t sumg = fGhness->Integral(1, n+1); - const Stat_t sump = fPhness->Integral(1, n+1); + Stat_t sumg; + Stat_t sump; + + sumg = fGhness->Integral(1, n); + sump = fPhness->Integral(1, n); + + if (sumg == 0.0 || sump == 0.0) + { + *fLog << "MHHadronness::Finalize; sumg or sump is zero; sumg, sump = " + << sumg << ", " << sump << ". Cannot calculate hadronness" + << endl; + } + + + // Normalize photon distribution + Stat_t con; + if (sumg > 0.0) + for (Int_t i=1; i<=n; i++) + { + con = (fGhness->GetBinContent(i)) / sumg; + fGhness->SetBinContent(i, con); + } + + // Normalize hadron distribution + if (sump > 0.0) + for (Int_t i=1; i<=n; i++) + { + con = (fPhness->GetBinContent(i)) / sump; + fPhness->SetBinContent(i, con); + } + + // Calculate acceptances + sumg = fGhness->Integral(1, n); + sump = fPhness->Integral(1, n); + + *fLog << "MHHadronness::Finalize; sumg, sump = " << sumg << ", " + << sump << endl; Float_t max=0; @@ -259,6 +300,15 @@ for (Int_t i=1; i<=n; i++) { - const Stat_t ip = fPhness->Integral(1, i)/sump; - const Stat_t ig = fGhness->Integral(1, i)/sumg; + Stat_t ip; + if (sump != 0.0) + ip = fPhness->Integral(1, i)/sump; + else + ip = 0; + + Stat_t ig; + if (sumg != 0.0) + ig = fGhness->Integral(1, i)/sumg; + else + ig = 0; fIntPhness->SetBinContent(i, ip); @@ -409,10 +459,10 @@ c.cd(1); - gStyle->SetOptStat(10); + //gStyle->SetOptStat(10); Getghness()->DrawCopy(); Getphness()->SetLineColor(kRed); Getphness()->DrawCopy("same"); c.cd(2); - gStyle->SetOptStat(0); + //gStyle->SetOptStat(0); Getighness()->DrawCopy(); Getiphness()->SetLineColor(kRed); @@ -483,10 +533,10 @@ gPad->cd(1); - gStyle->SetOptStat(10); + //gStyle->SetOptStat(10); Getghness()->Draw(); Getphness()->SetLineColor(kRed); Getphness()->Draw("same"); gPad->cd(2); - gStyle->SetOptStat(0); + //gStyle->SetOptStat(0); Getighness()->Draw(); Getiphness()->SetLineColor(kRed); Index: trunk/MagicSoft/Mars/mhist/MHMatrix.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1808) +++ trunk/MagicSoft/Mars/mhist/MHMatrix.cc (revision 1809) @@ -504,4 +504,6 @@ } +// -------------------------------------------------------------------------- +// void MHMatrix::Reassign() { @@ -691,6 +693,6 @@ // // Define the reference matrix -// refcolumn number of the column containing the variable, for which a -// target distribution may be given; +// refcolumn number of the column (starting at 1)containing the variable, +// for which a target distribution may be given; // if refcolumn is negative the target distribution will be set // equal to the real distribution; the events in the reference @@ -748,14 +750,15 @@ //--------------------------------------------------------- // - // if refcol < 0 : select reference events randomly - // i.e. set the normaliztion factotrs equal to 1.0 + // if refcolumn < 0 : select reference events randomly + // i.e. set the normaliztion factotrs equal to 1.0 + // refcol is the column number starting at 0; it is >= 0 if (refcolumn<0) { - frefcol = -refcolumn; + frefcol = -refcolumn - 1; } else { - frefcol = refcolumn; + frefcol = refcolumn - 1; } @@ -790,15 +793,4 @@ // << ", " << fM(j-1,frefcol) << endl; } - - // if refcolumn<0 set target distribution equal to the real distribution - // in order to obtain normalization factors = 1.0 - //if (refcolumn<0) - //{ - // for (Int_t j=1; j<=fnbins; j++) - // { - // float cont = fHth-> GetBinContent(j); - // fHthsh->SetBinContent(j, cont); - // } - //} //--------------------------------------------------------- @@ -816,5 +808,5 @@ // if refcolumn < 0 set the correction factors equal to 1 - if ( refcolumn>=0.0 ) + if ( refcolumn>=0 ) b = fHthsh->GetBinContent(j); else @@ -948,5 +940,5 @@ { *fLog <GetBinContent(j); - if (a>0) *fLog << j << " "<< a << " "; + *fLog << j << " "<< a << " "; } *fLog <cd(1); - ((TH1F*)fHthsh)->Draw(); // target + ((TH1F*)fHthsh)->DrawCopy(); // target th1->cd(2); - ((TH1F*)fHth)->Draw(); // real histogram before + ((TH1F*)fHth)->DrawCopy(); // real histogram before th1->cd(3); - ((TH1F*)fHthd) -> Draw(); // correction factors + ((TH1F*)fHthd)->DrawCopy(); // correction factors th1->cd(4); - ((TH1F*)fHthaft) -> Draw(); // histogram after - - //--------------------------------------------------------- - // --- write onto output file - // - //TFile *outfile = new TFile(fileNameout, "RECREATE", ""); - //Mnew.Write(fileNameout); - //outfile->Close(); + ((TH1F*)fHthaft)->DrawCopy(); // histogram after + + //--------------------------------------------------------- return kTRUE; Index: trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1808) +++ trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 1809) @@ -224,5 +224,5 @@ TObject *MHSigmaTheta::DrawClone(Option_t *opt) { - TCanvas &c = *MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", + TCanvas &c = *MakeDefCanvas("SigmaThetaPlot", "Sigmabar vs. Theta", 900, 900); c.Divide(3, 3); @@ -236,9 +236,10 @@ c.cd(1); h = ((TH2*)&fSigmaTheta)->ProjectionX("ProjX-Theta", -1, 9999, "E"); + h->SetDirectory(NULL); h->SetTitle("Distribution of \\Theta"); h->SetXTitle("\\Theta [\\circ]"); h->SetYTitle("No.of events"); - h->Draw(opt); + h->DrawCopy(opt); h->SetBit(kCanDelete);; gPad->SetLogy(); @@ -246,9 +247,10 @@ c.cd(4); h = ((TH2*)&fSigmaTheta)->ProjectionY("ProjY-sigma", -1, 9999, "E"); + h->SetDirectory(NULL); h->SetTitle("Distribution of Sigmabar"); h->SetXTitle("Sigmabar"); h->SetYTitle("No.of events"); - h->Draw(opt); + h->DrawCopy(opt); h->SetBit(kCanDelete);; @@ -263,18 +265,20 @@ c.cd(2); l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zx"); + l->SetDirectory(NULL); l->SetTitle("Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)"); l->SetXTitle("\\Theta [\\circ]"); l->SetYTitle("Sigma^2-Sigmabar^2"); - l->Draw("box"); + l->DrawCopy("box"); l->SetBit(kCanDelete);; c.cd(5); l = (TH2D*) ((TH3*)&fDiffPixTheta)->Project3D("zy"); + l->SetDirectory(NULL); l->SetTitle("Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)"); l->SetXTitle("pixel"); l->SetYTitle("Sigma^2-Sigmabar^2"); - l->Draw("box"); + l->DrawCopy("box"); l->SetBit(kCanDelete);; @@ -290,18 +294,20 @@ c.cd(3); k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zx"); + k->SetDirectory(NULL); k->SetTitle("Sigma vs. \\Theta (all pixels)"); k->SetXTitle("\\Theta [\\circ]"); k->SetYTitle("Sigma"); - k->Draw("box"); + k->DrawCopy("box"); k->SetBit(kCanDelete);; c.cd(6); k = (TH2D*) ((TH3*)&fSigmaPixTheta)->Project3D("zy"); + k->SetDirectory(NULL); k->SetTitle("Sigma vs. pixel number (all \\Theta)"); k->SetXTitle("pixel"); k->SetYTitle("Sigma"); - k->Draw("box"); + k->DrawCopy("box"); k->SetBit(kCanDelete);; Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 1808) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 1809) @@ -88,4 +88,5 @@ void Print(Option_t *opt=NULL) const; + UInt_t GetEvtNumber() const { return fEvtNumber; } //Get Event Number Short_t GetPartId() const { return fPartId; } //Get Type of particle Float_t GetEnergy() const { return fEnergy; } //Get Energy