Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4583) +++ /trunk/MagicSoft/Mars/Changelog (revision 4584) @@ -19,4 +19,36 @@ -*-*- END OF LINE -*-*- + + + + 2004/08/12 : Wolfgang Wittek + + * manalysis/MSigmabarCalc.[h,cc] + - replace MMcEvt by MPointingPos + + * manalysis/MSigmabar.[h,cc] + - in member function Calc() return fSigmabarInner, + not fSigmabar + - update comments + - sigmabar is the sqrt of the average (pedRMS^2/area) + + * manalysis/MPad.[h,cc] + - replace MMcEvt by MPointingPos + - remove bugs + + * mfilter/MFSelBasic.[h,cc] + - replace MMcEvt by MPointingPos + + * mfilter/Makefile + - add -I../mpointing + + + * mhist/MHSigmaTheta.[h,cc] + - replace MMcEvt by MPointingPos + - replace 'MCerPhotPix cerpix' by 'MCerPhotPix &cerpix' + - add plot "Sigmabar(Outer) versus Theta" + + * macros/ONOFFAnalysis.C + - Job A : got the padding working, work in progress Index: /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C =================================================================== --- /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 4583) +++ /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 4584) @@ -171,6 +171,6 @@ //----------------------------------------------- //TString tag = "040126"; - TString tag = "040127"; - //TString tag = "040215"; + //TString tag = "040127"; + TString tag = "040215"; //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; @@ -183,5 +183,6 @@ //const char *offfile = "*.OFF"; // 15 Feb 04 - const char *offfile = "*.OFF"; + //const char *offfile = "*.OFF"; + const char *offfile = "174*.OFF"; @@ -191,5 +192,5 @@ //const char *onfile = "MCerPhot_output"; //const char *onfile = "*.ON"; - const char *onfile = "1216*.ON"; + //const char *onfile = "1216*.ON"; // 27 Jan 04 //const char *onfile = "12*.ON"; @@ -198,9 +199,13 @@ // 15 Feb 04 //const char *onfile = "*.ON"; - - - const char *mcfile = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root"; + const char *onfile = "174*.ON"; + + + //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_3/gh/0/0/G_M0_00_0_550*.root"; //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_4/gh/0/0/G_M0_00_0_550*.root"; //const char *mcfile = "/data/MAGIC/mc_eth/magLQE_5/gh/0/0/G_M0_00_0_550*.root"; + //const char *mcfile = "calibrated_gammas"; + const char *mcfile = "calibrated_data_david*"; + //----------------------------------------------- @@ -224,12 +229,15 @@ // 15 Feb 04, David - if (tag == "040215") - TString inPath = "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/"; - + //if (tag == "040215") + // TString inPath = "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/"; + if (tag == "040215") + TString inPathON = "/.magic/magicserv01/scratch/Daniel/CalibratedData/MispointingTest/2004_02_15ForgottenData/finesam20/"; + TString inPathOFF = "/.magic/magicserv01/scratch/David/CalibratedData/Mkn421/2004_02_15/"; + TString inPathMC = "/.magic/magicserv01/scratch/David/MCData/MCApril2004/Data/gammas_highnoise/"; // path for output from Mars //TString outPath = "~wittek/datacrab_feb04/"; //TString outPath = "~wittek/datacrab_01march04/"; - TString outPath = "~wittek/datacrab_19mar04/"; + TString outPath = "~wittek/datacrab_27april04/"; //----------------------------------------------- @@ -248,9 +256,12 @@ Bool_t JobA = kTRUE; - Bool_t GPad = kFALSE; // generate padding histograms? - Bool_t WPad = kFALSE; // write out padding histograms ? - Bool_t RPad = kFALSE; // read in padding histograms ? - Bool_t Wout = kTRUE; // write out root file ON1.root - // (or OFF1.root or MC1.root)? + Bool_t GPadON = kFALSE; // \ generate padding histograms + Bool_t GPadOFF = kFALSE; // | and write them onto disk + Bool_t GPadMC = kFALSE; // / + Bool_t Merge = kFALSE; // merge padding histograms + // and write them onto disk + Bool_t Wout = kTRUE; // read in merged padding histograms and + // write out root file of padded data + // (ON1.root, OFF1.root or MC1.root) @@ -310,6 +321,6 @@ // 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 + // - read MC root file + // - calculate eff. collection area // - optimize energy estimator // - read ON root file @@ -346,9 +357,11 @@ gLog << "Macro ONOFFAnalysis : Start of Job A" << endl; gLog << "" << endl; - gLog << "Macro ONOFFAnalysis : JobA, WPad, RPad, Wout = " - << (JobA ? "kTRUE" : "kFALSE") << ", " - << (WPad ? "kTRUE" : "kFALSE") << ", " - << (RPad ? "kTRUE" : "kFALSE") << ", " - << (Wout ? "kTRUE" : "kFALSE") << endl; + gLog << "Macro ONOFFAnalysis : JobA, GPadON, GPadOFF, GPadMC, Merge, Wout = " + << (JobA ? "kTRUE" : "kFALSE") << ", " + << (GPadON ? "kTRUE" : "kFALSE") << ", " + << (GPadOFF ? "kTRUE" : "kFALSE") << ", " + << (GPadMC ? "kTRUE" : "kFALSE") << ", " + << (Merge ? "kTRUE" : "kFALSE") << ", " + << (Wout ? "kTRUE" : "kFALSE") << endl; @@ -358,13 +371,13 @@ - TString fileON = inPath; + TString fileON = inPathON; fileON += onfile; fileON += ".root"; - TString fileOFF = inPath; + TString fileOFF = inPathOFF; fileOFF += offfile; fileOFF += ".root"; - TString fileMC = mcfile; + TString fileMC = inPathMC; fileMC += mcfile; fileMC += ".root"; @@ -373,5 +386,5 @@ << fileOFF << ", " << fileMC << endl; - // name of file to conatin the histograms for the padding + // name of file to conatin the merged histograms for the padding TString outNameSigTh = outPath; outNameSigTh += "SigmaTheta"; @@ -379,8 +392,22 @@ //-------------------------------------------------- + // name of files to contain the paddding histograms of ON, OFF and MC data + TString NamePadON(outPath); + NamePadON += "PadON"; + NamePadON += ".root"; + + TString NamePadOFF(outPath); + NamePadOFF += "PadOFF"; + NamePadOFF += ".root"; + + TString NamePadMC(outPath); + NamePadMC += "PadMC"; + NamePadMC += ".root"; + + //-------------------------------------------------- // type of data to be padded - TString typeInput = "ON"; - //TString typeInput = "OFF"; - //TString typeInput = "MC"; + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -416,6 +443,6 @@ // // read ON, OFF and MC data files - // generate (or read in) the padding histograms for ON and OFF data - // and merge these histograms + // generate (or read in) the padding histograms for ON, OFF and MC data + // MPad pad; @@ -423,11 +450,16 @@ pad.SetDataType(typeInput); + if (GPadON || GPadOFF || GPadMC) + { // generate the padding histograms - if (GPad) - { gLog << "=====================================================" << endl; gLog << "Start generating the padding histograms" << endl; + } + //----------------------------------------- // ON events + + if (GPadON) + { gLog << "-----------" << endl; @@ -438,14 +470,17 @@ MParList pliston; + MObservatory observon; + MReadMarsFile readON("Events", fileON); readON.DisableAutoScheme(); - //MCT1ReadPreProc readON(fileON); - - //MFSelBasic selthetaon; - //selthetaon.SetCuts(-100.0, 29.5, 35.5); - //MContinue contthetaon(&selthetaon); - - MBlindPixelCalc blindon; - blindon.SetUseBlindPixels(); + + MGeomApply apply; + + MSourcePosfromStarPos sourcefromstaron; + sourcefromstaron.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); + + MBlindPixelsCalc2 blindon; + //blindon.SetUseBlindPixels(); + blindon.SetCheckPedestalRms(); MFSelBasic selbasicon; @@ -459,5 +494,5 @@ MHSigmaTheta sigthON("SigmaThetaON"); - MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MMcEvt"); + MFillH fillsigthetaON ("SigmaThetaON[MHSigmaTheta]", "MPointingPos"); fillsigthetaON.SetName("FillSigTheta"); @@ -467,4 +502,5 @@ pliston.AddToList(&tliston); InitBinnings(&pliston); + pliston.AddToList(&observon); pliston.AddToList(&blindON); pliston.AddToList(&sigthON); @@ -475,5 +511,6 @@ tliston.AddToList(&readON); - //tliston.AddToList(&contthetaon); + tliston.AddToList(&apply); + tliston.AddToList(&sourcefromstaron); tliston.AddToList(&blindon); @@ -489,6 +526,6 @@ evtloopon.SetProgressBar(&baron); - Int_t maxeventson = -1; - //Int_t maxeventson = 10000; + //Int_t maxeventson = -1; + Int_t maxeventson = 10000; if ( !evtloopon.Eventloop(maxeventson) ) return; @@ -507,6 +544,28 @@ TH2D *blindnthon = blindON.GetBlindN(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : write padding plots for ON data from file " + << NamePadON << endl; + + TFile filewriteon(NamePadON, "RECREATE", ""); + sigthon->Write(); + sigpixthon->Write(); + diffpixthon->Write(); + blindidthon->Write(); + blindnthon->Write(); + filewriteon.Close(); + + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for ON data written onto file " + << NamePadON << endl; + } + // end of GPadON + + //----------------------------------------- // OFF events + + if (GPadOFF) + { gLog << "------------" << endl; @@ -517,14 +576,15 @@ MParList plistoff; + MObservatory observoff; + MReadMarsFile readOFF("Events", fileOFF); readOFF.DisableAutoScheme(); - // MCT1ReadPreProc readOFF(fileOFF); - - MFSelBasic selthetaoff; - selthetaoff.SetCuts(-100.0, 29.5, 35.5); - MContinue contthetaoff(&selthetaoff); - - MBlindPixelCalc blindoff; - blindoff.SetUseBlindPixels(); + + MSourcePosfromStarPos sourcefromstaroff; + sourcefromstaroff.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); + + MBlindPixelsCalc2 blindoff; + //blindoff.SetUseBlindPixels(); + blindoff.SetCheckPedestalRms(); MFSelBasic selbasicoff; @@ -538,5 +598,5 @@ MHSigmaTheta sigthOFF("SigmaThetaOFF"); - MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MMcEvt"); + MFillH fillsigthetaOFF ("SigmaThetaOFF[MHSigmaTheta]", "MPointingPos"); fillsigthetaOFF.SetName("FillSigThetaOFF"); @@ -546,4 +606,5 @@ plistoff.AddToList(&tlistoff); InitBinnings(&plistoff); + plistoff.AddToList(&observoff); plistoff.AddToList(&blindOFF); plistoff.AddToList(&sigthOFF); @@ -554,5 +615,5 @@ tlistoff.AddToList(&readOFF); - //tlistoff.AddToList(&contthetaoff); + tlistoff.AddToList(&sourcefromstaroff); tlistoff.AddToList(&blindoff); @@ -568,6 +629,6 @@ evtloopoff.SetProgressBar(&baroff); - Int_t maxeventsoff = -1; - //Int_t maxeventsoff = 20000; + //Int_t maxeventsoff = -1; + Int_t maxeventsoff = 10000; if ( !evtloopoff.Eventloop(maxeventsoff) ) return; @@ -586,7 +647,29 @@ TH2D *blindnthoff = blindOFF.GetBlindN(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : write padding plots for OFF data from file " + << NamePadOFF << endl; + + TFile filewriteoff(NamePadOFF, "RECREATE", ""); + sigthoff->Write(); + sigpixthoff->Write(); + diffpixthoff->Write(); + blindidthoff->Write(); + blindnthoff->Write(); + filewriteoff.Close(); + + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for OFF data written onto file " + << NamePadOFF << endl; + } + // end of GPadOFF + + //----------------------------------------- // MC events + + if (GPadMC) + { gLog << "------------" << endl; @@ -597,14 +680,14 @@ MParList plistmc; + MObservatory observmc; + MReadMarsFile readMC("Events", fileMC); readMC.DisableAutoScheme(); - // MCT1ReadPreProc readMC(fileMC); - - MFSelBasic selthetamc; - selthetamc.SetCuts(-100.0, 29.5, 35.5); - MContinue contthetamc(&selthetamc); - - MBlindPixelCalc blindmc; - blindmc.SetUseBlindPixels(); + + MSourcePosfromStarPos sourcefromstarmc; + + MBlindPixelsCalc2 blindmc; + //blindmc.SetUseBlindPixels(); + blindmc.SetCheckPedestalRms(); MFSelBasic selbasicmc; @@ -618,5 +701,5 @@ MHSigmaTheta sigthMC("SigmaThetaMC"); - MFillH fillsigthetaMC ("SigmaThetaMC[MHSigmaTheta]", "MMcEvt"); + MFillH fillsigthetaMC ("SigmaThetaMC[MHSigmaTheta]", "MPointingPos"); fillsigthetaMC.SetName("FillSigThetaMC"); @@ -626,4 +709,5 @@ plistmc.AddToList(&tlistmc); InitBinnings(&plistmc); + plistmc.AddToList(&observmc); plistmc.AddToList(&blindMC); plistmc.AddToList(&sigthMC); @@ -634,5 +718,5 @@ tlistmc.AddToList(&readMC); - //tlistmc.AddToList(&contthetamc); + tlistmc.AddToList(&sourcefromstarmc); tlistmc.AddToList(&blindmc); @@ -648,6 +732,6 @@ evtloopmc.SetProgressBar(&barmc); - Int_t maxeventsmc = -1; - //Int_t maxeventsmc = 20000; + //Int_t maxeventsmc = -1; + Int_t maxeventsmc = 10000; if ( !evtloopmc.Eventloop(maxeventsmc) ) return; @@ -666,30 +750,150 @@ TH2D *blindnthmc = blindMC.GetBlindN(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : write padding plots for MC data from file " + << NamePadMC << endl; + + TFile filewritemc(NamePadMC, "RECREATE", ""); + sigthmc->Write(); + sigpixthmc->Write(); + diffpixthmc->Write(); + blindidthmc->Write(); + blindnthmc->Write(); + filewritemc.Close(); + + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for MC data written onto file " + << NamePadMC << endl; + } + // end of GPadMC //----------------------------------------- + if (GPadON || GPadOFF || GPadMC) + { + gLog << "" << endl; gLog << "End of generating the padding histograms" << endl; gLog << "=====================================================" << endl; - - pad.MergeONOFFMC(sigthmc, diffpixthmc, sigpixthmc, - blindidthmc, blindnthmc, - sigthon, diffpixthon, sigpixthon, - blindidthon, blindnthon, - sigthoff, diffpixthoff, sigpixthoff, - blindidthoff, blindnthoff); - - - if (WPad) - { - // write the padding histograms onto a file --------- - pad.WritePaddingDist(outNameSigTh); - } } - // read the padding histograms --------------------------- - if (RPad) + if (Merge) { - pad.ReadPaddingDist(outNameSigTh); + + //----------------------------------- + TH2D sigon; + TH3D sigpixon; + TH3D diffpixon; + TH2D blindidon; + TH2D blindnon; + + TFile filereadon(NamePadON); + filereadon.ls(); + sigon.Read("2D-ThetaSigmabar(Inner)"); + sigpixon.Read("3D-ThetaPixSigma"); + diffpixon.Read("3D-ThetaPixDiff"); + blindidon.Read("2D-IdBlindPixels"); + blindnon.Read("2D-NBlindPixels"); + + TH2D *sigthon = new TH2D( (const TH2D&)sigon ); + TH3D *sigpixthon = new TH3D( (const TH3D&)sigpixon ); + TH3D *diffpixthon = new TH3D( (const TH3D&)diffpixon ); + TH2D *blindidthon = new TH2D( (const TH2D&)blindidon ); + TH2D *blindnthon = new TH2D( (const TH2D&)blindnon ); + + /* + TH2D *sigthon = (TH2D*)sigon.Clone(); + TH3D *sigpixthon = (TH3D*)sigpixon.Clone(); + TH3D *diffpixthon = (TH3D*)diffpixon.Clone(); + TH2D *blindidthon = (TH2D*)blindidon.Clone(); + TH2D *blindnthon = (TH2D*)blindnon.Clone(); + */ + + //filereadon.Close(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for ON data read from file " + << NamePadON << endl; + + + //----------------------------------- + TH2D sigoff; + TH3D sigpixoff; + TH3D diffpixoff; + TH2D blindidoff; + TH2D blindnoff; + + TFile filereadoff(NamePadOFF); + filereadoff.ls(); + sigoff.Read("2D-ThetaSigmabar(Inner)"); + sigpixoff.Read("3D-ThetaPixSigma"); + diffpixoff.Read("3D-ThetaPixDiff"); + blindidoff.Read("2D-IdBlindPixels"); + blindnoff.Read("2D-NBlindPixels"); + + TH2D *sigthoff = new TH2D( (const TH2D&)sigoff ); + TH3D *sigpixthoff = new TH3D( (const TH3D&)sigpixoff ); + TH3D *diffpixthoff = new TH3D( (const TH3D&)diffpixoff ); + TH2D *blindidthoff = new TH2D( (const TH2D&)blindidoff ); + TH2D *blindnthoff = new TH2D( (const TH2D&)blindnoff ); + + /* + TH2D *sigthoff = (TH2D*)sigoff.Clone(); + TH3D *sigpixthoff = (TH3D*)sigpixoff.Clone(); + TH3D *diffpixthoff = (TH3D*)diffpixoff.Clone(); + TH2D *blindidthoff = (TH2D*)blindidoff.Clone(); + TH2D *blindnthoff = (TH2D*)blindnoff.Clone(); + */ + + //filereadoff.Close(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for OFF data read from file " + << NamePadOFF << endl; + + + //----------------------------------- + TH2D sigmc; + TH3D sigpixmc; + TH3D diffpixmc; + TH2D blindidmc; + TH2D blindnmc; + + TFile filereadmc(NamePadMC); + filereadmc.ls(); + sigmc.Read("2D-ThetaSigmabar(Inner)"); + sigpixmc.Read("3D-ThetaPixSigma"); + diffpixmc.Read("3D-ThetaPixDiff"); + blindidmc.Read("2D-IdBlindPixels"); + blindnmc.Read("2D-NBlindPixels"); + + TH2D *sigthmc = new TH2D( (const TH2D&)sigmc ); + TH3D *sigpixthmc = new TH3D( (const TH3D&)sigpixmc ); + TH3D *diffpixthmc = new TH3D( (const TH3D&)diffpixmc ); + TH2D *blindidthmc = new TH2D( (const TH2D&)blindidmc ); + TH2D *blindnthmc = new TH2D( (const TH2D&)blindnmc ); + + /* + TH2D *sigthmc = (TH2D*)sigmc.Clone(); + TH3D *sigpixthmc = (TH3D*)sigpixmc.Clone(); + TH3D *diffpixthmc = (TH3D*)diffpixmc.Clone(); + TH2D *blindidthmc = (TH2D*)blindidmc.Clone(); + TH2D *blindnthmc = (TH2D*)blindnmc.Clone(); + */ + + //filereadmc.Close(); + gLog << "" << endl; + gLog << "Macro ONOFFAnalysis : padding plots for MC data read from file " + << NamePadMC << endl; + + pad.MergeONOFFMC(*sigthmc, *diffpixthmc, *sigpixthmc, + *blindidthmc, *blindnthmc, + *sigthon, *diffpixthon, *sigpixthon, + *blindidthon, *blindnthon, + *sigthoff, *diffpixthoff,*sigpixthoff, + *blindidthoff,*blindnthoff); + + // write the target padding histograms onto a file --------- + pad.WritePaddingDist(outNameSigTh); } + // end of Merge + @@ -698,4 +902,8 @@ if (Wout) { + // read the target padding histograms --------------------------- + pad.ReadPaddingDist(outNameSigTh); + + gLog << "=====================================================" << endl; gLog << "Start the padding" << endl; @@ -723,58 +931,47 @@ MGeomApply apply; - //MPedestalWorkaround waround; - // a way to find out whether one is dealing with MC : - MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC - f1.SetName("Select MC"); - MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data - f2.SetName("Select Data"); - - //if (typeInput == "ON") - //{ - // MCT1PointingCorrCalc pointcorr(sourceName, "MCT1PointingCorrCalc", - // "MCT1PointingCorrCalc"); - //} + //MFDataMember f1("MRawRunHeader.fRunType", '>', 255.5); // MC + //f1.SetName("Select MC"); + //MFDataMember f2("MRawRunHeader.fRunType", '<', 255.5); // data + //f2.SetName("Select Data"); + MSourcePosfromStarPos sourcefromstar; if (typeInput == "ON") { - sourcefromstar.SetSourceAndStarPosition( - "Crab", 22, 0, 52, 5, 34, 32.0, - "Zeta-Tau", 21, 8, 33, 5, 37, 38.7); - sourcefromstar.AddStar("Tau114", 21, 56, 13, 5, 27, 38.1); - //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOn.4.txt",0); - - if(inPath == "/.magic/magicserv01/scratch/calibrated/") - { - sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions2stars.txt", 0); - //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsNoStar.txt", 0); - } - else if(inPath == "/.magic/magicserv01/scratch/calibrated26/") - sourcefromstar.AddFile("~wittek/datacrab_26feb04/stars_2004_01_26", 0); - - else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/") - sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); + //sourcefromstar.SetSourceAndStarPosition( + // "Crab", 22, 0, 52, 5, 34, 32.0, + // "Zeta-Tau", 21, 8, 33, 5, 37, 38.7); + //sourcefromstar.AddStar("Tau114", 21, 56, 13, 5, 27, 38.1); + ////sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOn.4.txt",0); + + //if(inPath == "/.magic/magicserv01/scratch/calibrated/") + //{ + //sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions2stars.txt", 0); + ////sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsNoStar.txt", 0); + //} + //else if(inPath == "/.magic/magicserv01/scratch/calibrated26/") + // sourcefromstar.AddFile("~wittek/datacrab_26feb04/stars_2004_01_26", 0); + + //else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15/") + sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); } else if (typeInput == "OFF") { - if(inPath == "/.magic/magicserv01/scratch/calibrated/") - sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOff.txt", 0); - else if(inPath == "/.magic/magicserv01/scratch/calibrated26/") - sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions_2004_01_26", 0); - else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15") - sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); + //if(inPath == "/.magic/magicserv01/scratch/calibrated/") + // sourcefromstar.AddFile("~wittek/datacrab_26feb04/positionsOff.txt", 0); + //else if(inPath == "/.magic/magicserv01/scratch/calibrated26/") + // sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions_2004_01_26", 0); + //else if(inPath == "/mnt/magicserv01/scratch/David/CalibratedData/Crab/2004_02_15") + sourcefromstar.AddFile("~wittek/datacrab_26feb04/positions040215.txt", 0); } - - - //MBlindPixelCalc blindbeforepad; + else if (typeInput == "MC") + { + } + + MBlindPixelsCalc2 blindbeforepad; //blindbeforepad.SetUseBlindPixels(); - //blindbeforepad.SetName("BlindBeforePadding"); - - //MBlindPixelCalc blind; - //blind.SetUseBlindPixels(); - //blind.SetUseInterpolation(); - //blind.SetName("BlindAfterPadding"); - - MSigmabarCalc sigbar; + blindbeforepad.SetCheckPedestalRms(); + blindbeforepad.SetName("BlindBeforePadding"); //MBadPixelCalcRms blind; @@ -794,5 +991,5 @@ MSigmabarCalc sigbarcalc; - MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MMcEvt"); + MFillH fillsigtheta ("SigmaTheta[MHSigmaTheta]", "MPointingPos"); fillsigtheta.SetName("HSigmaTheta"); @@ -845,6 +1042,5 @@ //write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE); //write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); - //write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MPointingPos", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); @@ -868,22 +1064,17 @@ tliston.AddToList(&read); - //tliston.AddToList(&f1); //tliston.AddToList(&f2); tliston.AddToList(&apply); - //tliston.AddToList(&waround); tliston.AddToList(&sourcefromstar); - //tliston.AddToList(&blindbeforepad); - // tliston.AddToList(&pad); - // if (typeInput == "ON") - // tliston.AddToList(&pointcorr); - - tliston.AddToList(&sigbar); + tliston.AddToList(&blindbeforepad); + + tliston.AddToList(&contbasic); + tliston.AddToList(&pad); + tliston.AddToList(&blind); - tliston.AddToList(&contbasic); - tliston.AddToList(&fillblind); - //tliston.AddToList(&sigbarcalc); + tliston.AddToList(&sigbarcalc); tliston.AddToList(&fillsigtheta); tliston.AddToList(&clean); @@ -913,6 +1104,6 @@ // evtloop.Write(); - Int_t maxevents = -1; - //Int_t maxevents = 1000; + //Int_t maxevents = -1; + Int_t maxevents = 2000; if ( !evtloop.Eventloop(maxevents) ) return; @@ -1109,5 +1300,5 @@ matrixg.EnableGraphicalOutput(); - matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)"); + matrixg.AddColumn("cos((MPointingPos.fZd)/kRad2Deg)"); matrixg.AddColumn("MSigmabar.fSigmabar"); matrixg.AddColumn("log10(MHillas.fSize)"); @@ -1198,5 +1389,5 @@ bing.SetEdges(10, 0., 1.0); - MH3 gref("cos(MMcEvt.fTelescopeTheta)"); + MH3 gref("cos((MPointingPos.fZd)/kRad2Deg)"); gref.SetName(mgname); MH::SetBinning(&gref.GetHist(), &bing); @@ -1231,5 +1422,5 @@ writetraing.AddContainer("MRawRunHeader", "RunHeaders"); writetraing.AddContainer("MTime", "Events"); - writetraing.AddContainer("MMcEvt", "Events"); + writetraing.AddContainer("MPointingPos", "Events"); writetraing.AddContainer("ThetaOrig", "Events"); writetraing.AddContainer("MSrcPosCam", "Events"); @@ -1248,5 +1439,5 @@ writetestg.AddContainer("MRawRunHeader", "RunHeaders"); writetestg.AddContainer("MTime", "Events"); - writetestg.AddContainer("MMcEvt", "Events"); + writetestg.AddContainer("MPointingPos", "Events"); writetestg.AddContainer("ThetaOrig", "Events"); writetestg.AddContainer("MSrcPosCam", "Events"); @@ -1331,5 +1522,5 @@ binh.SetEdges(10, 0., 1.0); - //MH3 href("cos(MMcEvt.fTelescopeTheta)"); + //MH3 href("cos((MPointingPos.fZd)/kRad2Deg)"); //href.SetName(mhname); //MH::SetBinning(&href.GetHist(), &binh); @@ -1366,5 +1557,5 @@ writetrainh.AddContainer("MRawRunHeader", "RunHeaders"); writetrainh.AddContainer("MTime", "Events"); - writetrainh.AddContainer("MMcEvt", "Events"); + writetrainh.AddContainer("MPointingPos", "Events"); writetrainh.AddContainer("ThetaOrig", "Events"); writetrainh.AddContainer("MSrcPosCam", "Events"); @@ -1382,5 +1573,5 @@ writetesth.AddContainer("MRawRunHeader", "RunHeaders"); writetesth.AddContainer("MTime", "Events"); - writetesth.AddContainer("MMcEvt", "Events"); + writetesth.AddContainer("MPointingPos", "Events"); writetesth.AddContainer("ThetaOrig", "Events"); writetesth.AddContainer("MSrcPosCam", "Events"); @@ -1667,5 +1858,5 @@ write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MPointingPos", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); @@ -1992,5 +2183,5 @@ bin.SetEdges(10, 0., 1.0); - MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); + MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)"); mh3.SetName(mname); MH::SetBinning(&mh3.GetHist(), &bin); @@ -2355,7 +2546,6 @@ write.AddContainer("MRawRunHeader", "RunHeaders"); - //write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); - //write.AddContainer("ThetaOrig", "Events"); + write.AddContainer("MTime", "Events"); + write.AddContainer("MPointingPos", "Events"); write.AddContainer("MSrcPosCam", "Events"); write.AddContainer("MSigmabar", "Events"); @@ -3137,5 +3327,5 @@ - MFillH filler("MHMcCollectionArea", "MMcEvt"); + MFillH filler("MHMcCollectionArea", "MPointingPos"); filler.SetName("CollectionArea"); @@ -3326,5 +3516,5 @@ write.AddContainer("MRawRunHeader", "RunHeaders"); write.AddContainer("MTime", "Events"); - write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MPointingPos", "Events"); write.AddContainer("ThetaOrig", "Events"); write.AddContainer("MSrcPosCam", "Events"); @@ -3400,8 +3590,8 @@ // new from Robert - MFillH hfill6("MHTimeDiffTheta", "MMcEvt"); + MFillH hfill6("MHTimeDiffTheta", "MPointingPos"); hfill6.SetName("HTimeDiffTheta"); - MFillH hfill6a("MHTimeDiffTime", "MMcEvt"); + MFillH hfill6a("MHTimeDiffTime", "MPointingPos"); hfill6a.SetName("HTimeDiffTime"); @@ -3611,24 +3801,24 @@ - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MPad.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 4583) +++ /trunk/MagicSoft/Mars/manalysis/MPad.cc (revision 4584) @@ -64,8 +64,9 @@ #include #include +#include #include "MBinning.h" #include "MSigmabar.h" -#include "MMcEvt.hxx" +#include "MPointingPos.h" #include "MLog.h" #include "MLogManip.h" @@ -75,4 +76,6 @@ #include "MCerPhotPix.h" #include "MCerPhotEvt.h" + +#include "MH.h" #include "MPedPhotCam.h" @@ -190,4 +193,6 @@ // from sigma_k to sigma_j + //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl; + Double_t eps = 1.e-10; @@ -205,7 +210,16 @@ Double_t v, s, w, t, x, u, a, b, arg; + //*fLog << "Content of histcp : " << endl; + for (Int_t j=nbinssig; j >= 1; j--) { + //*fLog << "j, hista, histb , histap : " << j << ", " + // << hista->GetBinContent(j) << ", " + // << histb->GetBinContent(j) << ", " + // << histap->GetBinContent(j) << endl; + v = histcp->GetBinContent(j); + //*fLog << "cp : j, v = " << j << ", " << v << endl; + if ( fabs(v) < eps ) continue; if (v >= 0.0) @@ -255,6 +269,6 @@ } - *fLog << all << "k, j, arg = " << k << ", " << j - << ", " << arg << endl; + //*fLog << all << "k, j, arg = " << k << ", " << j + // << ", " << arg << endl; //...................................... @@ -302,4 +316,7 @@ // check results + + //*fLog << "Content of histap, histbp, histcp : " << endl; + for (Int_t j=1; j<=nbinssig; j++) { @@ -308,4 +325,7 @@ Double_t c = histcp->GetBinContent(j); + //*fLog << "j, a, b, c = " << j << ": " << a << ", " << b << ", " + // << c << endl; + if( fabs(a-b)>3.0*eps || fabs(c)>3.0*eps ) *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = " @@ -314,22 +334,28 @@ //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("Merging", "", 600, 600)); - c.Divide(2, 2); + TCanvas *pad = MH::MakeDefCanvas("Merging", "", 600, 600); gROOT->SetSelectedPad(NULL); - c.cd(1); + pad->Divide(2, 2); + + pad->cd(1); + gPad->SetBorderMode(0); hista->SetDirectory(NULL); - hista->DrawCopy(); + hista->Draw(); hista->SetBit(kCanDelete); - c.cd(2); + pad->cd(2); + gPad->SetBorderMode(0); histb->SetDirectory(NULL); - histb->DrawCopy(); + histb->Draw(); histb->SetBit(kCanDelete); - c.cd(3); + pad->cd(3); + gPad->SetBorderMode(0); histap->SetDirectory(NULL); - histap->DrawCopy(); + histap->Draw(); histap->SetBit(kCanDelete); + + pad->DrawClone(); //-------------------------------------------------------------------- @@ -339,4 +365,6 @@ delete histcp; + //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl; + return kTRUE; } @@ -357,69 +385,125 @@ Bool_t MPad::MergeONOFFMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon, - TH2D *sigthoff, TH3D *diffpixthoff, TH3D *sigmapixthoff, - TH2D *blindidthoff, TH2D *blindnthoff) + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon, + TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff, + TH2D& blindidthoff, TH2D& blindnthoff) { + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + // for testing + /* + TAxis *xa = sigthon.GetXaxis(); + Int_t nbitheta = xa->GetNbins(); + + TAxis *ya = sigthon.GetYaxis(); + Int_t nbisig = ya->GetNbins(); + for (Int_t i=1; i<=nbitheta; i++) + { + for (Int_t j=1; j<=nbisig; j++) + { + if (i>0) + { + sigthmc.SetBinContent( i, j, (Float_t) (625000.0+(nbisig*nbisig*nbisig-j*j*j)) ); + sigthon.SetBinContent( i, j, (Float_t)(1.0) ); + sigthoff.SetBinContent( i, j, (Float_t)( (0.5*nbisig-j)*(0.5*nbisig-j) ) ); + } + else + { + sigthmc.SetBinContent( i, j, 0.0); + sigthon.SetBinContent( i, j, 0.0); + sigthoff.SetBinContent(i, j, 0.0); + } + } + } + */ + //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + *fLog << all << "----------------------------------------------------------------------------------" << endl; *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding" << endl; - fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); + fHSigmaTheta = new TH2D( (const TH2D&) sigthon ); fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); + *fLog << "after booking fHSigmaTheta" << endl; + + fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon ); fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); + *fLog << "after booking fHDiffPixTheta" << endl; + + fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon ); fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); + *fLog << "after booking fHSigmaPixTheta" << endl; + + fHBlindPixNTheta = new TH2D( (const TH2D&) blindnthon ); fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); - fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); + *fLog << "after booking fHBlindPixNTheta" << endl; + + fHBlindPixIdTheta = new TH2D( (const TH2D&) blindidthon ); fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + *fLog << "after booking fHBlindPixIdTheta" << endl; + *fLog << all << "Histograms for the merged padding plots were booked" + << endl; + //-------------------------- - fHSigmaThetaMC = sigthmc; + fHSigmaThetaMC = &sigthmc; fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = sigthon; + fHSigmaThetaON = &sigthon; fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - fHSigmaThetaOFF = sigthoff; + fHSigmaThetaOFF = &sigthoff; fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)"); + //*fLog << all << "addresses of SigmaTheta padding plots were copied" + // << endl; + //-------------------------- - fHBlindPixNThetaMC = blindnthmc; + fHBlindPixNThetaMC = &blindnthmc; fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); - fHBlindPixNThetaON = blindnthon; + fHBlindPixNThetaON = &blindnthon; fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); - fHBlindPixNThetaOFF = blindnthoff; + fHBlindPixNThetaOFF = &blindnthoff; fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)"); + //*fLog << all << "addresses of BlindPixTheta padding plots were copied" + // << endl; + //-------------------------- - fHBlindPixIdThetaMC = blindidthmc; + fHBlindPixIdThetaMC = &blindidthmc; fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); - fHBlindPixIdThetaON = blindidthon; + fHBlindPixIdThetaON = &blindidthon; fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); - fHBlindPixIdThetaOFF = blindidthoff; + fHBlindPixIdThetaOFF = &blindidthoff; fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)"); + //*fLog << all << "addresses of BlindPixIdTheta padding plots were copied" + // << endl; + //-------------------------- - fHDiffPixThetaMC = diffpixthmc; + fHDiffPixThetaMC = &diffpixthmc; fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = diffpixthon; + fHDiffPixThetaON = &diffpixthon; fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - fHDiffPixThetaOFF = diffpixthoff; + fHDiffPixThetaOFF = &diffpixthoff; fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)"); + //*fLog << all << "addresses of DiffPixTheta padding plots were copied" + // << endl; + //-------------------------- - fHSigmaPixThetaMC = sigmapixthmc; + fHSigmaPixThetaMC = &sigmapixthmc; fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = sigmapixthon; + fHSigmaPixThetaON = &sigmapixthon; fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - fHSigmaPixThetaOFF = sigmapixthoff; + fHSigmaPixThetaOFF = &sigmapixthoff; fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)"); + + //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" + // << endl; + //-------------------------- @@ -427,5 +511,553 @@ // get binning for fHgON, fHgOFF and fHgMC from sigthon // binning in Theta - TAxis *ax = sigthon->GetXaxis(); + TAxis *ax = sigthon.GetXaxis(); + Int_t nbinstheta = ax->GetNbins(); + TArrayD edgesx; + edgesx.Set(nbinstheta+1); + for (Int_t i=0; i<=nbinstheta; i++) + { + edgesx[i] = ax->GetBinLowEdge(i+1); + //*fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] + // << endl; + } + MBinning binth; + binth.SetEdges(edgesx); + + // binning in sigmabar + TAxis *ay = sigthon.GetYaxis(); + Int_t nbinssig = ay->GetNbins(); + TArrayD edgesy; + edgesy.Set(nbinssig+1); + for (Int_t i=0; i<=nbinssig; i++) + { + edgesy[i] = ay->GetBinLowEdge(i+1); + //*fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] + // << endl; + } + MBinning binsg; + binsg.SetEdges(edgesy); + + + fHgMC = new TH3D; + MH::SetBinning(fHgMC, &binth, &binsg, &binsg); + fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); + + fHgON = new TH3D; + MH::SetBinning(fHgON, &binth, &binsg, &binsg); + fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); + + fHgOFF = new TH3D; + MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); + fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); + + //*fLog << all << "fHg histograms were booked" << endl; + + //............ loop over Theta bins ........................... + + + + //*fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; + for (Int_t l=1; l<=nbinstheta; l++) + { + //------------------------------------------- + // merge ON and OFF distributions + // input : hista is the normalized 1D distr. of sigmabar for ON data + // histb is the normalized 1D distr. of sigmabar for OFF data + // output : histap will be the merged distribution (ON-OFF) + // fHga(k,j) will tell which fraction of ON events should be + // padded from bin k to bin j + // fHgb(k,j) will tell which fraction of OFF events should be + // padded from bin k to bin j + + TH2D * fHga = new TH2D; + MH::SetBinning(fHga, &binsg, &binsg); + TH2D * fHgb = new TH2D; + MH::SetBinning(fHgb, &binsg, &binsg); + + TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); + hista->SetNameTitle("ON-orig", "ON (orig)"); + Stat_t suma = hista->Integral(); + if (suma != 0.0) hista->Scale(1./suma); + + TH1D *histap = new TH1D( (const TH1D&)*hista ); + histap->SetNameTitle("ON-OFF-merged)", "ON-OFF (merged)"); + + TH1D *histb = sigthoff.ProjectionY("sigoff_y", l, l, ""); + histb->SetNameTitle("OFF-orig", "OFF (orig)"); + Stat_t sumb = histb->Integral(); + if (sumb != 0.0) histb->Scale(1./sumb); + + if (suma == 0.0 || sumb == 0.0) + { + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + *fLog << err + << "cannot call Merge2Distributions(ON, OFF) for theta bin l = " + << l << "; NevON, NevOFF = " << suma << ", " << sumb << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(ON, OFF) for theta bin l = " + // << l << endl; + + Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); + + // fill fHgON and fHgOFF + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t a = fHga->GetBinContent(k,j); + fHgON->SetBinContent(l, k, j, a); + + Double_t b = fHgb->GetBinContent(k,j); + fHgOFF->SetBinContent(l, k, j, b); + } + + //------------------------------------------- + // now merge ON-OFF and MC distributions + // input : histe is the result of the merging of ON and OFF distributions + // histf is the normalized 1D distr. of sigmabar for MC data + // output : histep will be the merged distribution (target distribution) + // fHge(k,j) will tell which fraction of ON, OFF events should be + // padded from bin k to bin j + // fHgf(k,j) will tell which fraction of MC events should be + // padded from bin k to bin j + + TH2D * fHge = new TH2D; + MH::SetBinning(fHge, &binsg, &binsg); + TH2D * fHgf = new TH2D; + MH::SetBinning(fHgf, &binsg, &binsg); + + TH1D *histe = new TH1D( (const TH1D&)*histap); + histe->SetNameTitle("ON-OFF-merged", "ON-OFF (merged)"); + + TH1D *histep = new TH1D( (const TH1D&)*histe); + histep->SetNameTitle("ON-OFF-MC-merged)", "ON-OFF-MC (merged)"); + + TH1D *histf = sigthmc.ProjectionY("sigmc_y", l, l, ""); + histf->SetNameTitle("MC-orig", "MC (orig)"); + Double_t sumf = histf->Integral(); + if (sumf != 0.0) histf->Scale(1./sumf); + + if (sumf == 0.0) + { + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + delete histe; + delete histf; + delete fHge; + delete fHgf; + delete histep; + + *fLog << err + << "cannot call Merge2Distributions(ON-OFF,MC) for theta bin l = " + << l << "; NevMC = " << sumf << endl; + + // go to next theta bin (l) + continue; + } + + //*fLog << "call Merge2Distributions(ON-OFF, MC) for theta bin l = " + // << l << endl; + + Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig); + + // update fHgON and fHgOFF + UpdateHg(fHga, histap, fHge, fHgON, nbinssig, l); + UpdateHg(fHgb, histap, fHge, fHgOFF, nbinssig, l); + + // fill fHgMC + for (Int_t k=1; k<=nbinssig; k++) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t f = fHgf->GetBinContent(k,j); + fHgMC->SetBinContent(l, k, j, f); + } + + // fill target distribution fHSigmaTheta + // (only for plotting, it will not be used in the padding) + for (Int_t j=1; j<=nbinssig; j++) + { + Double_t ep = histep->GetBinContent(j); + fHSigmaTheta->SetBinContent(l, j, ep); + } + //------------------------------------------- + + delete hista; + delete histb; + delete fHga; + delete fHgb; + delete histap; + + delete histe; + delete histf; + delete fHge; + delete fHgf; + delete histep; + } + //............ end of loop over Theta bins .................... + + + // Define the target distributions + // fHDiffPixTheta, fHSigmaPixTheta + // (they are calculated as + // averages of the ON and OFF distributions) + // fHBlindPixNTheta, fHBlindPixIdTheta + // (they are calculated as + // the sum of the ON and OFF distributions) + // (fHDiffPixTheta will be used in the padding of MC events) + + //............ start of new loop over Theta bins .................... + for (Int_t j=1; j<=nbinstheta; j++) + { + // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC + + Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); + + Double_t normON; + Double_t normOFF; + + // set ranges for 2D-projections of 3D-histograms + TAxis *ax = diffpixthon.GetXaxis(); + ax->SetRange(j, j); + ax = diffpixthoff.GetXaxis(); + ax->SetRange(j, j); + + TH2D *hist; + TH2D *histOFF; + + TH1D *hist1; + TH1D *hist1OFF; + + // weights for ON and OFF distrubtions when + // calculating the weighted average + Double_t facON = 0.5 * (1. - fracON + fracOFF); + Double_t facOFF = 0.5 * (1. + fracON - fracOFF); + + *fLog << all << "Theta bin j : fracON, fracOFF, fracMC, facON, facOFF = " + << j << ", " << fracON << ", " << fracOFF << ", " + << fracMC << ", " << facON << ", " << facOFF << endl; + + //------------------------------------------------------------------ + // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' + ay = diffpixthon.GetYaxis(); + Int_t nbinspixel = ay->GetNbins(); + + TAxis *az = diffpixthon.GetZaxis(); + Int_t nbinsdiff = az->GetNbins(); + + hist = (TH2D*)diffpixthon.Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)diffpixthoff.Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinsdiff; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHDiffPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + //------------------------------------------------------------------ + // define target distribution 'sigma vs. pixel number' + ay = sigmapixthon.GetYaxis(); + nbinspixel = ay->GetNbins(); + + az = sigmapixthon.GetZaxis(); + Int_t nbinssigma = az->GetNbins(); + + hist = (TH2D*)sigmapixthon.Project3D("zy"); + hist->SetName("dummy"); + histOFF = (TH2D*)sigmapixthoff.Project3D("zy"); + + normON = hist->Integral(); + normOFF = histOFF->Integral(); + if (normON != 0.0) hist->Scale(1.0/normON); + if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); + + // weighted average of ON and OFF distributions + hist->Scale(facON); + hist->Add(histOFF, facOFF); + + for (Int_t k=1; k<=nbinspixel; k++) + for (Int_t l=1; l<=nbinssigma; l++) + { + Double_t cont = hist->GetBinContent(k,l); + fHSigmaPixTheta->SetBinContent(j, k, l, cont); + } + + delete hist; + delete histOFF; + + + //------------------------------------------------------------------ + // define target distribution 'number of blind pixels per event' + ay = blindnthon.GetYaxis(); + Int_t nbinsn = ay->GetNbins(); + + hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, ""); + + normON = hist1->Integral(); + normOFF = hist1OFF->Integral(); + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + Stat_t sum1 = hist1->Integral(); + if (sum1 != 0.0) hist1->Scale( 1.0/sum1 ); + + for (Int_t k=1; k<=nbinsn; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixNTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + + //------------------------------------------------------------------ + // define target distribution 'id of blind pixel' + ay = blindidthon.GetYaxis(); + Int_t nbinsid = ay->GetNbins(); + + hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); + hist1->SetName("dummy"); + hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, ""); + + // divide by the number of events (from fHBlindPixNTheta) + if (normON != 0.0) hist1->Scale(1.0/normON); + if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); + + // sum of ON and OFF distributions + hist1->Add(hist1OFF, 1.0); + + for (Int_t k=1; k<=nbinsid; k++) + { + Double_t cont = hist1->GetBinContent(k); + fHBlindPixIdTheta->SetBinContent(j, k, cont); + } + + delete hist1; + delete hist1OFF; + } + //............ end of new loop over Theta bins .................... + + *fLog << all + << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" + << endl; + *fLog << all << "----------------------------------------------------------" + << endl; + //-------------------------------------------- + + fHSigmaThetaMC->SetDirectory(NULL); + fHSigmaThetaON->SetDirectory(NULL); + fHSigmaThetaOFF->SetDirectory(NULL); + + fHDiffPixThetaMC->SetDirectory(NULL); + fHDiffPixThetaON->SetDirectory(NULL); + fHDiffPixThetaOFF->SetDirectory(NULL); + + fHSigmaPixThetaMC->SetDirectory(NULL); + fHSigmaPixThetaON->SetDirectory(NULL); + fHSigmaPixThetaOFF->SetDirectory(NULL); + + fHBlindPixIdThetaMC->SetDirectory(NULL); + fHBlindPixIdThetaON->SetDirectory(NULL); + fHBlindPixIdThetaOFF->SetDirectory(NULL); + + fHBlindPixNThetaMC->SetDirectory(NULL); + fHBlindPixNThetaON->SetDirectory(NULL); + fHBlindPixNThetaOFF->SetDirectory(NULL); + + fHgMC->SetDirectory(NULL); + fHgON->SetDirectory(NULL); + fHgOFF->SetDirectory(NULL); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Update the matrix fHgA +// which contains the final padding prescriptions +// + +Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, + Int_t nbinssig, Int_t l) +{ + // histap target distribution after ON-OFF merging + // fHga padding prescription for ON (or OFF) data to get to histap + // fHge padding prescription to get from histap to the target + // distribution after the ON-OFF-MC merging + // fHgA updated padding prescription for ON (or OFF) data to get + // from the original ON (or OFF) histogram to the target + // distribution after the ON-OFF-MC merging + // l Theta bin + + Double_t eps = 1.e-10; + + for (Int_t k=1; k<=nbinssig; k++) + { + Double_t na = fHga->Integral(1, nbinssig, k, k, " "); + Double_t ne = fHge->Integral(k, k, 1, nbinssig, " "); + Double_t nap = histap->GetBinContent(k); + + if (ne <= eps) + { + // go to next k + continue; + } + + Double_t r = nap - na; + + if (r >= ne-eps) + { + for (Int_t j=k+1; j<=nbinssig; j++) + { + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,k,j); + fHgA->SetBinContent(l, k, j, A + e); + } + // go to next k + continue; + } + + for (Int_t j=k+1; j<=nbinssig; j++) + { + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,k,j); + fHgA->SetBinContent(l, k, j, A + r*e/ne); + } + + if (na <= eps) + { + // go to next k + continue; + } + + for (Int_t i=1; iGetBinContent(i,k); + Double_t A = fHgA->GetBinContent(l,i,k); + fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na); + } + + for (Int_t i=1; i<=nbinssig; i++) + { + if (i == k) continue; + for (Int_t j=i+1; j<=nbinssig; j++) + { + if (j == k) continue; + Double_t a = fHga->GetBinContent(i,k); + Double_t e = fHge->GetBinContent(k,j); + Double_t A = fHgA->GetBinContent(l,i,j); + fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne) ); + } + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Merge the distributions +// fHSigmaTheta 2D-histogram (Theta, sigmabar) +// +// ===> of ON and MC data <============== +// +// and define the matrices fHgMC and fHgON +// +// These matrices tell which fraction of events should be padded +// from bin k of sigmabar to bin j +// + +Bool_t MPad::MergeONMC( + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon) +{ + *fLog << all << "----------------------------------------------------------------------------------" << endl; + *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" + << endl; + + fHSigmaTheta = new TH2D( (TH2D&) sigthon ); + fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); + + fHDiffPixTheta = new TH3D( (TH3D&) diffpixthon ); + fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); + + fHSigmaPixTheta = new TH3D( (TH3D&) sigmapixthon ); + fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); + + fHBlindPixNTheta = new TH2D( (TH2D&) blindnthon ); + fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); + + fHBlindPixIdTheta = new TH2D( (TH2D&) blindidthon ); + fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); + + //-------------------------- + fHSigmaThetaMC = &sigthmc; + fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); + fHSigmaThetaON = &sigthon; + fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); + + //-------------------------- + fHBlindPixNThetaMC = &blindnthmc; + fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); + fHBlindPixNThetaON = &blindnthon; + fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); + + //-------------------------- + fHBlindPixIdThetaMC = &blindidthmc; + fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); + fHBlindPixIdThetaON = &blindidthon; + fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); + + //-------------------------- + fHDiffPixThetaMC = &diffpixthmc; + fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); + fHDiffPixThetaON = &diffpixthon; + fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); + + //-------------------------- + fHSigmaPixThetaMC = &sigmapixthmc; + fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); + fHSigmaPixThetaON = &sigmapixthon; + fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); + + //-------------------------- + + + // get binning for fHgON and fHgMC from sigthon + // binning in Theta + TAxis *ax = sigthon.GetXaxis(); Int_t nbinstheta = ax->GetNbins(); TArrayD edgesx; @@ -441,5 +1073,5 @@ // binning in sigmabar - TAxis *ay = sigthon->GetYaxis(); + TAxis *ay = sigthon.GetYaxis(); Int_t nbinssig = ay->GetNbins(); TArrayD edgesy; @@ -463,411 +1095,4 @@ fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - fHgOFF = new TH3D; - MH::SetBinning(fHgOFF, &binth, &binsg, &binsg); - fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF"); - - //............ loop over Theta bins ........................... - - - - *fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl; - for (Int_t l=1; l<=nbinstheta; l++) - { - TH2D * fHga = new TH2D; - MH::SetBinning(fHga, &binsg, &binsg); - TH2D * fHgb = new TH2D; - MH::SetBinning(fHgb, &binsg, &binsg); - - //------------------------------------------- - // merge ON and OFF distributions - // input : hista is the normalized 1D distr. of sigmabar for ON data - // histb is the normalized 1D distr. of sigmabar for OFF data - // output : histap will be the merged distribution (ON-OFF) - // fHga(k,j) will tell which fraction of ON events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of OFF events should be - // padded from bin k to bin j - - TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); - Stat_t suma = hista->Integral(); - hista->Scale(1./suma); - - TH1D *histap = new TH1D( (const TH1D&)*hista ); - - TH1D *histb = sigthoff->ProjectionY("sigoff_y", l, l, ""); - Stat_t sumb = histb->Integral(); - histb->Scale(1./sumb); - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - delete hista; - delete histb; - - // fill fHgON and fHgOFF - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent(k,j); - fHga->SetBinContent (k, j, 0.0); - fHgON->SetBinContent(l, k, j, a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgb->SetBinContent (k, j, 0.0); - fHgOFF->SetBinContent(l, k, j, b); - } - - //------------------------------------------- - // now merge ON-OFF and MC distributions - // input : hista is the result of the merging of ON and OFF distributions - // histb is the normalized 1D distr. of sigmabar for MC data - // output : histap will be the merged distribution (target distribution) - // fHga(k,j) will tell which fraction of ON, OFF events should be - // padded from bin k to bin j - // fHgb(k,j) will tell which fraction of MC events should be - // padded from bin k to bin j - - hista = new TH1D( (const TH1D&)*histap); - - histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); - sumb = histb->Integral(); - histb->Scale(1./sumb); - - Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); - delete hista; - delete histb; - - // update fHgON and fHgOFF and fill fHgMC - for (Int_t k=1; k<=nbinssig; k++) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = fHga->GetBinContent (k,j); - - Double_t aON = fHgON->GetBinContent(l,k,j); - fHgON->SetBinContent(l, k, j, aON+a); - - Double_t aOFF = fHgOFF->GetBinContent(l,k,j); - fHgOFF->SetBinContent(l, k, j, aOFF+a); - - Double_t b = fHgb->GetBinContent(k,j); - fHgMC->SetBinContent(l, k, j, b); - } - - // fill target distribution fHSigmaTheta - // (only for plotting, it will not be used in the padding) - for (Int_t j=1; j<=nbinssig; j++) - { - Double_t a = histap->GetBinContent(j); - fHSigmaTheta->SetBinContent(l, j, a); - } - - delete fHga; - delete fHgb; - - delete histap; - } - //............ end of loop over Theta bins .................... - - - // Define the target distributions - // fHDiffPixTheta, fHSigmaPixTheta - // (they are calculated as - // averages of the ON and OFF distributions) - // fHBlindPixNTheta, fHBlindPixIdTheta - // (they are calculated as - // the sum of the ON and OFF distributions) - // (fHDiffPixTheta will be used in the padding of MC events) - - //............ start of new loop over Theta bins .................... - for (Int_t j=1; j<=nbinstheta; j++) - { - // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC - - Double_t fracON = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - Double_t fracMC = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, ""); - - Double_t normON; - Double_t normOFF; - - // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon->GetXaxis(); - ax->SetRange(j, j); - ax = diffpixthoff->GetXaxis(); - ax->SetRange(j, j); - - TH2D *hist; - TH2D *histOFF; - - TH1D *hist1; - TH1D *hist1OFF; - - // weights for ON and OFF distrubtions when - // calculating the weighted average - Double_t facON = 0.5 * (1. - fracON + fracOFF); - Double_t facOFF = 0.5 * (1. + fracON - fracOFF); - - *fLog << all << "Theta bin j : fracON, fracOFF, facON, facOFF = " - << j << ", " << fracON << ", " << fracOFF << ", " - << fracMC << ", " << facON << ", " << facOFF << endl; - - //------------------------------------------------------------------ - // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon->GetYaxis(); - Int_t nbinspixel = ay->GetNbins(); - - TAxis *az = diffpixthon->GetZaxis(); - Int_t nbinsdiff = az->GetNbins(); - - hist = (TH2D*)diffpixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)diffpixthoff->Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinsdiff; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHDiffPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - //------------------------------------------------------------------ - // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon->GetYaxis(); - nbinspixel = ay->GetNbins(); - - az = sigmapixthon->GetZaxis(); - Int_t nbinssigma = az->GetNbins(); - - hist = (TH2D*)sigmapixthon->Project3D("zy"); - hist->SetName("dummy"); - histOFF = (TH2D*)sigmapixthoff->Project3D("zy"); - - normON = hist->Integral(); - normOFF = histOFF->Integral(); - if (normON != 0.0) hist->Scale(1.0/normON); - if (normOFF != 0.0) histOFF->Scale(1.0/normOFF); - - // weighted average of ON and OFF distributions - hist->Scale(facON); - hist->Add(histOFF, facOFF); - - for (Int_t k=1; k<=nbinspixel; k++) - for (Int_t l=1; l<=nbinssigma; l++) - { - Double_t cont = hist->GetBinContent(k,l); - fHSigmaPixTheta->SetBinContent(j, k, l, cont); - } - - delete hist; - delete histOFF; - - - //------------------------------------------------------------------ - // define target distribution 'number of blind pixels per event' - ay = blindnthon->GetYaxis(); - Int_t nbinsn = ay->GetNbins(); - - hist1 = fHBlindPixNThetaON->ProjectionY("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, ""); - - normON = hist1->Integral(); - normOFF = hist1OFF->Integral(); - if (normON != 0.0) hist1->Scale(1.0/normON); - if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); - - // sum of ON and OFF distributions - hist1->Add(hist1OFF, 1.0); - hist1->Scale( 1.0/hist1->Integral() ); - - for (Int_t k=1; k<=nbinsn; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixNTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - - //------------------------------------------------------------------ - // define target distribution 'id of blind pixel' - ay = blindidthon->GetYaxis(); - Int_t nbinsid = ay->GetNbins(); - - hist1 = fHBlindPixIdThetaON->ProjectionY("", j, j, ""); - hist1->SetName("dummy"); - hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, ""); - - // divide by the number of events (from fHBlindPixNTheta) - if (normON != 0.0) hist1->Scale(1.0/normON); - if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF); - - // sum of ON and OFF distributions - hist1->Add(hist1OFF, 1.0); - - for (Int_t k=1; k<=nbinsid; k++) - { - Double_t cont = hist1->GetBinContent(k); - fHBlindPixIdTheta->SetBinContent(j, k, cont); - } - - delete hist1; - delete hist1OFF; - } - //............ end of new loop over Theta bins .................... - - *fLog << all - << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" - << endl; - *fLog << all << "----------------------------------------------------------" - << endl; - //-------------------------------------------- - - fHSigmaThetaMC->SetDirectory(NULL); - fHSigmaThetaON->SetDirectory(NULL); - fHSigmaThetaOFF->SetDirectory(NULL); - - fHDiffPixThetaMC->SetDirectory(NULL); - fHDiffPixThetaON->SetDirectory(NULL); - fHDiffPixThetaOFF->SetDirectory(NULL); - - fHSigmaPixThetaMC->SetDirectory(NULL); - fHSigmaPixThetaON->SetDirectory(NULL); - fHSigmaPixThetaOFF->SetDirectory(NULL); - - fHBlindPixIdThetaMC->SetDirectory(NULL); - fHBlindPixIdThetaON->SetDirectory(NULL); - fHBlindPixIdThetaOFF->SetDirectory(NULL); - - fHBlindPixNThetaMC->SetDirectory(NULL); - fHBlindPixNThetaON->SetDirectory(NULL); - fHBlindPixNThetaOFF->SetDirectory(NULL); - - fHgMC->SetDirectory(NULL); - fHgON->SetDirectory(NULL); - fHgOFF->SetDirectory(NULL); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Merge the distributions -// fHSigmaTheta 2D-histogram (Theta, sigmabar) -// -// ===> of ON and MC data <============== -// -// and define the matrices fHgMC and fHgON -// -// These matrices tell which fraction of events should be padded -// from bin k of sigmabar to bin j -// - -Bool_t MPad::MergeONMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon) -{ - *fLog << all << "----------------------------------------------------------------------------------" << endl; - *fLog << all << "MPad::MergeONMC(); Merge the ON and MC histograms to obtain the target distributions for the padding" - << endl; - - fHSigmaTheta = new TH2D( (TH2D&)*sigthon ); - fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)"); - - fHDiffPixTheta = new TH3D( (TH3D&)*diffpixthon ); - fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)"); - - fHSigmaPixTheta = new TH3D( (TH3D&)*sigmapixthon ); - fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)"); - - fHBlindPixNTheta = new TH2D( (TH2D&)*blindnthon ); - fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)"); - - fHBlindPixIdTheta = new TH2D( (TH2D&)*blindidthon ); - fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)"); - - //-------------------------- - fHSigmaThetaMC = sigthmc; - fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)"); - fHSigmaThetaON = sigthon; - fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)"); - - //-------------------------- - fHBlindPixNThetaMC = blindnthmc; - fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)"); - fHBlindPixNThetaON = blindnthon; - fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)"); - - //-------------------------- - fHBlindPixIdThetaMC = blindidthmc; - fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)"); - fHBlindPixIdThetaON = blindidthon; - fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)"); - - //-------------------------- - fHDiffPixThetaMC = diffpixthmc; - fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)"); - fHDiffPixThetaON = diffpixthon; - fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)"); - - //-------------------------- - fHSigmaPixThetaMC = sigmapixthmc; - fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)"); - fHSigmaPixThetaON = sigmapixthon; - fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)"); - - //-------------------------- - - - // get binning for fHgON and fHgMC from sigthon - // binning in Theta - TAxis *ax = sigthon->GetXaxis(); - Int_t nbinstheta = ax->GetNbins(); - TArrayD edgesx; - edgesx.Set(nbinstheta+1); - for (Int_t i=0; i<=nbinstheta; i++) - { - edgesx[i] = ax->GetBinLowEdge(i+1); - *fLog << all << "i, theta low edge = " << i << ", " << edgesx[i] - << endl; - } - MBinning binth; - binth.SetEdges(edgesx); - - // binning in sigmabar - TAxis *ay = sigthon->GetYaxis(); - Int_t nbinssig = ay->GetNbins(); - TArrayD edgesy; - edgesy.Set(nbinssig+1); - for (Int_t i=0; i<=nbinssig; i++) - { - edgesy[i] = ay->GetBinLowEdge(i+1); - *fLog << all << "i, sigmabar low edge = " << i << ", " << edgesy[i] - << endl; - } - MBinning binsg; - binsg.SetEdges(edgesy); - - - fHgMC = new TH3D; - MH::SetBinning(fHgMC, &binth, &binsg, &binsg); - fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC"); - - fHgON = new TH3D; - MH::SetBinning(fHgON, &binth, &binsg, &binsg); - fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON"); - //............ loop over Theta bins ........................... @@ -893,13 +1118,13 @@ // padded from bin k to bin j - TH1D *hista = sigthon->ProjectionY("sigon_y", l, l, ""); + TH1D *hista = sigthon.ProjectionY("sigon_y", l, l, ""); Stat_t suma = hista->Integral(); - hista->Scale(1./suma); + if (suma != 0.0) hista->Scale(1./suma); TH1D *histap = new TH1D( (const TH1D&)*hista ); - TH1D *histb = sigthmc->ProjectionY("sigmc_y", l, l, ""); + TH1D *histb = sigthmc.ProjectionY("sigmc_y", l, l, ""); Stat_t sumb = histb->Integral(); - histb->Scale(1./sumb); + if (sumb != 0.0) histb->Scale(1./sumb); Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig); @@ -958,7 +1183,7 @@ // set ranges for 2D-projections of 3D-histograms - TAxis *ax = diffpixthon->GetXaxis(); + TAxis *ax = diffpixthon.GetXaxis(); ax->SetRange(j, j); - ax = diffpixthmc->GetXaxis(); + ax = diffpixthmc.GetXaxis(); ax->SetRange(j, j); @@ -980,13 +1205,13 @@ //------------------------------------------------------------------ // define target distribution 'sigma^2-sigmavar^2 vs. pixel number' - ay = diffpixthon->GetYaxis(); + ay = diffpixthon.GetYaxis(); Int_t nbinspixel = ay->GetNbins(); - TAxis *az = diffpixthon->GetZaxis(); + TAxis *az = diffpixthon.GetZaxis(); Int_t nbinsdiff = az->GetNbins(); - hist = (TH2D*)diffpixthon->Project3D("zy"); + hist = (TH2D*)diffpixthon.Project3D("zy"); hist->SetName("dummy"); - histMC = (TH2D*)diffpixthmc->Project3D("zy"); + histMC = (TH2D*)diffpixthmc.Project3D("zy"); normON = hist->Integral(); @@ -1011,18 +1236,18 @@ //------------------------------------------------------------------ // define target distribution 'sigma vs. pixel number' - ay = sigmapixthon->GetYaxis(); + ay = sigmapixthon.GetYaxis(); nbinspixel = ay->GetNbins(); - az = sigmapixthon->GetZaxis(); + az = sigmapixthon.GetZaxis(); Int_t nbinssigma = az->GetNbins(); - hist = (TH2D*)sigmapixthon->Project3D("zy"); + hist = (TH2D*)sigmapixthon.Project3D("zy"); hist->SetName("dummy"); - histMC = (TH2D*)sigmapixthmc->Project3D("zy"); + histMC = (TH2D*)sigmapixthmc.Project3D("zy"); normON = hist->Integral(); normMC = histMC->Integral(); if (normON != 0.0) hist->Scale(1.0/normON); - if (normMC != 0.0) histMC->Scale(1.0/normMC); + if (normMC != 0.0) histMC->Scale(1.0/normMC); // weighted average of ON and MC distributions @@ -1043,5 +1268,5 @@ //------------------------------------------------------------------ // define target distribution 'number of blind pixels per event' - ay = blindnthon->GetYaxis(); + ay = blindnthon.GetYaxis(); Int_t nbinsn = ay->GetNbins(); @@ -1053,9 +1278,10 @@ normMC = hist1MC->Integral(); if (normON != 0.0) hist1->Scale(1.0/normON); - if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); // sum of ON and MC distributions hist1->Add(hist1MC, 1.0); - hist1->Scale( 1.0/hist1->Integral() ); + Stat_t sum1 = hist1->Integral(); + if (sum1 != 0.0) hist1->Scale( 1.0/sum1 ); for (Int_t k=1; k<=nbinsn; k++) @@ -1070,5 +1296,5 @@ //------------------------------------------------------------------ // define target distribution 'id of blind pixel' - ay = blindidthon->GetYaxis(); + ay = blindidthon.GetYaxis(); Int_t nbinsid = ay->GetNbins(); @@ -1079,5 +1305,5 @@ // divide by the number of events (from fHBlindPixNTheta) if (normON != 0.0) hist1->Scale(1.0/normON); - if (normMC != 0.0) hist1MC->Scale(1.0/normMC); + if (normMC != 0.0) hist1MC->Scale(1.0/normMC); // sum of ON and MC distributions @@ -1138,4 +1364,6 @@ //------------------------------------ + + /* fHSigmaThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC"); @@ -1147,7 +1375,11 @@ return kFALSE; } + */ + fHSigmaThetaMC = new TH2D; + fHSigmaThetaMC->Read("2D-ThetaSigmabarMC"); *fLog << all << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl; + /* fHSigmaThetaON = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON"); @@ -1159,7 +1391,11 @@ return kFALSE; } + */ + fHSigmaThetaON = new TH2D; + fHSigmaThetaON->Read("2D-ThetaSigmabarON"); *fLog << all << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl; + /* fHSigmaThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF"); @@ -1171,8 +1407,12 @@ return kFALSE; } + */ + fHSigmaThetaOFF = new TH2D; + fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF"); *fLog << all << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl; //------------------------------------ + /* fHDiffPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixDiff"); @@ -1184,7 +1424,11 @@ return kFALSE; } + */ + fHDiffPixTheta = new TH3D; + fHDiffPixTheta->Read("3D-ThetaPixDiff"); *fLog << all << "MPad : Object '3D-ThetaPixDiff' was read in" << endl; + /* fHDiffPixThetaMC = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC"); @@ -1196,7 +1440,11 @@ return kFALSE; } + */ + fHDiffPixThetaMC = new TH3D; + fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC"); *fLog << all << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl; + /* fHDiffPixThetaON = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON"); @@ -1208,7 +1456,11 @@ return kFALSE; } + */ + fHDiffPixThetaON = new TH3D; + fHDiffPixThetaON->Read("3D-ThetaPixDiffON"); *fLog << all << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl; + /* fHDiffPixThetaOFF = (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF"); @@ -1220,4 +1472,7 @@ return kFALSE; } + */ + fHDiffPixThetaOFF = new TH3D; + fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF"); *fLog << all << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl; @@ -1225,4 +1480,5 @@ //------------------------------------ + /* fHSigmaPixTheta = (TH3D*) gROOT->FindObject("3D-ThetaPixSigma"); @@ -1234,7 +1490,11 @@ return kFALSE; } + */ + fHSigmaPixTheta = new TH3D; + fHSigmaPixTheta->Read("3D-ThetaPixSigma"); *fLog << all << "MPad : Object '3D-ThetaPixSigma' was read in" << endl; + /* fHSigmaPixThetaMC = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC"); @@ -1246,7 +1506,11 @@ return kFALSE; } + */ + fHSigmaPixThetaMC = new TH3D; + fHSigmaPixThetaMC->Read("3D-ThetaPixSigmaMC"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl; + /* fHSigmaPixThetaON = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON"); @@ -1258,7 +1522,11 @@ return kFALSE; } + */ + fHSigmaPixThetaON = new TH3D; + fHSigmaPixThetaON->Read("3D-ThetaPixSigmaON"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl; + /* fHSigmaPixThetaOFF = (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF"); @@ -1270,8 +1538,12 @@ return kFALSE; } + */ + fHSigmaPixThetaOFF = new TH3D; + fHSigmaPixThetaOFF->Read("3D-ThetaPixSigmaOFF"); *fLog << all << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl; //------------------------------------ + /* fHgMC = (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC"); @@ -1283,7 +1555,11 @@ return kFALSE; } + */ + fHgMC = new TH3D; + fHgMC->Read("3D-PaddingMatrixMC"); *fLog << all << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl; + /* fHgON = (TH3D*) gROOT->FindObject("3D-PaddingMatrixON"); @@ -1295,7 +1571,11 @@ return kFALSE; } + */ + fHgON = new TH3D; + fHgON->Read("3D-PaddingMatrixON"); *fLog << all << "MPad : Object '3D-PaddingMatrixON' was read in" << endl; + /* fHgOFF = (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF"); @@ -1307,8 +1587,28 @@ return kFALSE; } + */ + fHgOFF = new TH3D; + fHgOFF->Read("3D-PaddingMatrixOFF"); *fLog << all << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl; //------------------------------------ + /* + fHBlindPixIdTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindId"); + if (!fHBlindPixIdTheta) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindId' not found on root file" + << endl; + return kFALSE; + } + */ + fHBlindPixIdTheta = new TH2D; + fHBlindPixIdTheta->Read("2D-ThetaBlindId"); + *fLog << all + << "MPad : Object '2D-ThetaBlindId' was read in" << endl; + + /* fHBlindPixIdThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC"); @@ -1320,7 +1620,11 @@ return kFALSE; } + */ + fHBlindPixIdThetaMC = new TH2D; + fHBlindPixIdThetaMC->Read("2D-ThetaBlindIdMC"); *fLog << all << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + /* fHBlindPixIdThetaON = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON"); @@ -1332,7 +1636,11 @@ return kFALSE; } + */ + fHBlindPixIdThetaON = new TH2D; + fHBlindPixIdThetaON->Read("2D-ThetaBlindIdON"); *fLog << all << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl; + /* fHBlindPixIdThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF"); @@ -1344,8 +1652,28 @@ return kFALSE; } + */ + fHBlindPixIdThetaOFF = new TH2D; + fHBlindPixIdThetaOFF->Read("2D-ThetaBlindIdOFF"); *fLog << all - << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl; + << "MPad : Object '2D-ThetaBlindIdOFF' was read in" << endl; //------------------------------------ + /* + fHBlindPixNTheta = + (TH2D*) gROOT->FindObject("2D-ThetaBlindN"); + if (!fHBlindPixNTheta) + { + *fLog << all + << "MPad : Object '2D-ThetaBlindN' not found on root file" + << endl; + return kFALSE; + } + */ + fHBlindPixNTheta = new TH2D; + fHBlindPixNTheta->Read("2D-ThetaBlindN"); + *fLog << all + << "MPad : Object '2D-ThetaBlindN' was read in" << endl; + + /* fHBlindPixNThetaMC = (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC"); @@ -1357,7 +1685,11 @@ return kFALSE; } + */ + fHBlindPixNThetaMC = new TH2D; + fHBlindPixNThetaMC->Read("2D-ThetaBlindNMC"); *fLog << all << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl; + /* fHBlindPixNThetaON = (TH2D*) gROOT->FindObject("2D-ThetaBlindNON"); @@ -1369,7 +1701,11 @@ return kFALSE; } + */ + fHBlindPixNThetaON = new TH2D; + fHBlindPixNThetaON->Read("2D-ThetaBlindNON"); *fLog << all << "MPad : Object '2D-ThetaBlindNON' was read in" << endl; + /* fHBlindPixNThetaOFF = (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF"); @@ -1381,4 +1717,7 @@ return kFALSE; } + */ + fHBlindPixNThetaOFF = new TH2D; + fHBlindPixNThetaOFF->Read("2D-ThetaBlindNOFF"); *fLog << all << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl; @@ -1402,8 +1741,10 @@ TFile outfile(namefileout, "RECREATE"); + fHBlindPixNTheta->Write(); fHBlindPixNThetaMC->Write(); fHBlindPixNThetaON->Write(); fHBlindPixNThetaOFF->Write(); + fHBlindPixIdTheta->Write(); fHBlindPixIdThetaMC->Write(); fHBlindPixIdThetaON->Write(); @@ -1470,11 +1811,11 @@ } - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) + fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointPos) { - *fLog << err << dbginf << "MPad : MMcEvt not found... aborting." + *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." << endl; return kFALSE; - } + } fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam"); @@ -1510,6 +1851,6 @@ } - fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); - if (!fBlinds) + fBlindPix = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); + if (!fBlindPix) { *fLog << err << "MPad : MBlindPixels not found... aborting." @@ -1532,5 +1873,5 @@ // plot pedestal sigmas fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", - 100, 0.0, 5.0, 100, 0.0, 5.0); + 100, 0.0, 75.0, 100, 0.0, 75.0); fHSigmaPedestal->SetXTitle("Pedestal sigma before padding"); fHSigmaPedestal->SetYTitle("Pedestal sigma after padding"); @@ -1538,10 +1879,10 @@ // 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); + 100, -100.0, 300.0, 100, -100, 300); 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 = new TH1D("NSB","Distribution of added NSB", 100, -100.0, 200.0); fHNSB->SetXTitle("no.of added NSB photons"); fHNSB->SetYTitle("no.of pixels"); @@ -1550,5 +1891,5 @@ //-------------------------------------------------------------------- - fIter = 20; + fIter = 10; memset(fErrors, 0, sizeof(fErrors)); @@ -1570,5 +1911,5 @@ Int_t MPad::Process() { - *fLog << all << "Entry MPad::Process();" << endl; + //*fLog << all << "Entry MPad::Process();" << endl; // this is the conversion factor from the number of photons @@ -1576,10 +1917,16 @@ // later to be taken from MCalibration // fPEperPhoton = PW * LG * QE * 1D - Double_t fPEperPhoton = 0.95 * 0.90 * 0.13 * 0.9; + Double_t fPEperPhoton = 0.92 * 0.94 * 0.23 * 0.9; //------------------------------------------------ - // add pixels to MCerPhotEvt which are not yet in; - // set their number of photons equal to zero - + // Add pixels to MCerPhotEvt which are not yet in; + // set their number of photons equal to zero. + // Purpose : by the padding the number of photons is changed + // so that cleaning and cuts may be affected. + // However, no big effect is expectec unless the padding is strong + // (strong increase of the noise level) + // At present comment out this code + + /* UInt_t imaxnumpix = fCam->GetNumPixels(); @@ -1603,5 +1950,5 @@ fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms()); } - + */ @@ -1615,12 +1962,22 @@ // Calculate average sigma of the event // - Double_t sigbarold_phot = fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPad::Process(); before padding : " << endl; - fSigmabar->Print(""); - Double_t sigbarold = sigbarold_phot * fPEperPhoton; - Double_t sigbarold2 = sigbarold*sigbarold; - - const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); - *fLog << all << "theta = " << theta << endl; + + fSigmabar->Calc(*fCam, *fPed, *fEvt); + //*fLog << all << "MPad::Process(); before padding : " << endl; + //fSigmabar->Print(""); + + // inner sigmabar/sqrt(area) + Double_t sigbarInnerold_phot = fSigmabar->GetSigmabarInner(); + Double_t sigbarInnerold = sigbarInnerold_phot * fPEperPhoton; + Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold; + + // outer sigmabar/sqrt(area) + Double_t sigbarOuterold_phot = fSigmabar->GetSigmabarOuter(); + Double_t sigbarOuterold = sigbarOuterold_phot * fPEperPhoton; + Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold; + + const Double_t theta = fPointPos->GetZd(); + + //*fLog << all << "theta = " << theta << endl; Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta); @@ -1629,6 +1986,5 @@ *fLog << warn << "MPad::Process(); binNumber out of range : theta, binTheta = " - << theta << ", " << binTheta << "; Skip event " << endl; - // event cannot be padded; skip event + << theta << ", " << binTheta << "; aborting " << endl; rc = 2; @@ -1637,7 +1993,10 @@ } + //*fLog << all << "binTheta = " << binTheta << endl; + + //------------------------------------------- // get number of events in this theta bin - // and number of events in this sigbarold_phot bin + // and number of events in this sigbarInnerold_phot bin Double_t nTheta; @@ -1657,16 +2016,24 @@ TH1D *hn; hn = st->ProjectionY("", binTheta, binTheta, ""); + + //*fLog << "Title of histogram : " << st->GetTitle() << endl; + //for (Int_t i=1; i<=hn->GetNbinsX(); i++) + //{ + // *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i) + // << endl; + //} + nTheta = hn->Integral(); - Int_t binSigma = hn->FindBin(sigbarold_phot); + Int_t binSigma = hn->FindBin(sigbarInnerold_phot); nSigma = hn->GetBinContent(binSigma); - *fLog << all - << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " - << theta << ", " << sigbarold_phot << ", " << binTheta << ", " - << binSigma << ", " << nTheta << ", " << nSigma << endl; + //*fLog << all + // << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = " + // << theta << ", " << sigbarInnerold_phot << ", " << binTheta + // << ", " + // << binSigma << ", " << nTheta << ", " << nSigma << endl; delete hn; - - + //------------------------------------------- // for the current theta : @@ -1688,6 +2055,7 @@ if ( nblind->Integral() == 0.0 ) { - *fLog << warn << "MPad::Process(); projection for Theta bin " - << binTheta << " has no entries; Skip event " << endl; + *fLog << warn << "MPad::Process(); projection of '" + << fHBlindPixNThetaOFF->GetName() << "' for Theta bin " + << binTheta << " has no entries; aborting " << endl; // event cannot be padded; skip event delete nblind; @@ -1732,8 +2100,8 @@ nlist++; - fBlinds->SetPixelBlind(idBlind); - *fLog << all << "idBlind = " << idBlind << endl; + fBlindPix->SetPixelBlind(idBlind); + //*fLog << all << "idBlind = " << idBlind << endl; } - fBlinds->SetReadyToSave(); + fBlindPix->SetReadyToSave(); delete hblind; @@ -1752,7 +2120,7 @@ if ( nblind->Integral() == 0.0 ) { - *fLog << warn << "MPad::Process(); projection for Theta bin " + *fLog << warn << "MPad::Process(); projection of '" + << fHBlindPixNThetaON->GetName() << "' for Theta bin " << binTheta << " has no entries; Skip event " << endl; - // event cannot be padded; skip event delete nblind; @@ -1796,8 +2164,8 @@ nlist++; - fBlinds->SetPixelBlind(idBlind); - *fLog << all << "idBlind = " << idBlind << endl; + fBlindPix->SetPixelBlind(idBlind); + //*fLog << all << "idBlind = " << idBlind << endl; } - fBlinds->SetReadyToSave(); + fBlindPix->SetReadyToSave(); delete hblind; @@ -1810,7 +2178,7 @@ // - Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarold_phot); - *fLog << all << "binSig, sigbarold_phot = " << binSig << ", " - << sigbarold_phot << endl; + Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot); + //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ", " + // << sigbarInnerold_phot << endl; Double_t prob; @@ -1843,6 +2211,6 @@ prob = 0.0; - *fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma - << ", " << prob << endl; + //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ", " << nSigma + // << ", " << prob << endl; if ( prob <= 0.0 || gRandom->Uniform() > prob ) @@ -1850,5 +2218,5 @@ delete hpad; // event should not be padded - *fLog << all << "event is not padded" << endl; + //*fLog << all << "event is not padded" << endl; rc = 8; @@ -1857,5 +2225,5 @@ } // event should be padded - *fLog << all << "event will be padded" << endl; + //*fLog << all << "event will be padded" << endl; @@ -1864,6 +2232,6 @@ // for MC/ON/OFF data according to the matrix fHgMC/ON/OFF // - Double_t sigmabar_phot = 0; - Double_t sigmabar = 0; + Double_t sigmabarInner_phot = 0; + Double_t sigmabarInner = 0; @@ -1876,23 +2244,25 @@ //} - sigmabar_phot = hpad->GetRandom(); - sigmabar = sigmabar_phot * fPEperPhoton; - - *fLog << all << "sigmabar_phot = " << sigmabar_phot << endl; + sigmabarInner_phot = hpad->GetRandom(); + sigmabarInner = sigmabarInner_phot * fPEperPhoton; + + //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl; delete hpad; - const Double_t sigmabar2 = sigmabar*sigmabar; + // new inner sigmabar2/area + const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner; //------------------------------------------- - *fLog << all << "MPad::Process(); sigbarold, sigmabar = " - << sigbarold << ", "<< sigmabar << endl; - - // Skip event if target sigmabar is <= sigbarold - if (sigmabar <= sigbarold) + //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " + // << sigbarInnerold << ", "<< sigmabarInner << endl; + + // Stop if target sigmabar is <= sigbarold + if (sigmabarInner <= sigbarInnerold) { - *fLog << all << "MPad::Process(); target sigmabar is less than sigbarold : " - << sigmabar << ", " << sigbarold << ", Skip event" << endl; + *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : " + << sigmabarInner << ", " << sigbarInnerold << ", aborting" + << endl; rc = 4; @@ -1923,16 +2293,21 @@ } - // 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 + // In this Theta bin, get RMS of (Sigma^2-sigmabar^2) for inner pixels + // The average electronic noise (to be added) has to be in the order of + // this RMS, otherwise the electronic noise of an individual pixel + // (elNoise2Pix) may become negative + + // Attention : maximum ID of inner pixels hard coded !!! + Int_t idmaxpixinner = 396; + Int_t binmaxpixinner = + fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)idmaxpixinner ); TH1D *hnoise; if (fType == "MC") - hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "ON") - hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else if (fType == "OFF") - hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, 9999, ""); + hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, ""); else { @@ -1946,8 +2321,8 @@ delete hnoise; - elNoise2 = TMath::Min(RMS, sigmabar2 - sigbarold2); - *fLog << all << "elNoise2 = " << elNoise2 << endl; - - lambdabar = (sigmabar2 - sigbarold2 - elNoise2) / F2excess; + elNoise2 = TMath::Min(2.0*RMS, sigmabarInner2 - sigbarInnerold2); + //*fLog << all << "elNoise2 = " << elNoise2 << endl; + + lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess; // This value of lambdabar is the same for all pixels; @@ -1957,5 +2332,5 @@ // do the padding for each pixel // - // pad only pixels - which are used (before image cleaning) + // pad only pixels - which are used and not blind (before image cleaning) // @@ -1971,6 +2346,29 @@ + // throw the Sigma for the pixels from the distribution fHDiffPixTheta + // MC : from fHDiffPixTheta + // ON : from fHDiffPixThetaOFF + // OFF : from fHDiffPixThetaON + + TH3D *sp = NULL; + + if (fType == "MC") sp = fHDiffPixTheta; + else if (fType == "ON") sp = fHDiffPixThetaOFF; + else if (fType == "OFF") sp = fHDiffPixThetaON; + else + { + *fLog << err << "MPad::Process(); illegal data type... aborting" + << endl; + return kERROR; + } + + Double_t sigbarold2; + Double_t sigmabar2; + Double_t sigmabarOuter2; + + for (UInt_t i=0; iGetPixRatio(j); + if (ratioArea < 1.5) + { + sigbarold2 = sigbarInnerold2; + sigmabar2 = sigmabarInner2; + } + else + { + sigbarold2 = sigbarOuterold2; + //sigbarOuter2 = sigmabarInner2 * sigbarOuterold2 / sigbarInnerold2; + sigmabarOuter2 = sigmabarInner2 + (sigbarOuterold2 - sigbarInnerold2); + sigmabar2 = sigmabarOuter2; + } + MPedPhotPix &ppix = (*fPed)[j]; + + if ( fBlindPix != NULL && fBlindPix->IsBlind(j) ) + { + // this should never occur, because blind pixels should have + // been set unused by MBlindPixelsCalc2::UnMap() + //*fLog << all << "MPad::Process; blind pixel found which is used, j = " + // << j << "... go to next pixel." << endl; + continue; + } + + // count number of pixels treated + fWarnings[0]++; + + Double_t oldsigma_phot = ppix.GetRms(); Double_t oldsigma = oldsigma_phot * fPEperPhoton; - Double_t oldsigma2 = oldsigma*oldsigma; + Double_t oldsigma2 = oldsigma*oldsigma / ratioArea; //--------------------------------- @@ -2005,69 +2430,64 @@ TH1D *hdiff; - // throw the Sigma for this pixel from the distribution fHDiffPixTheta - // MC : from fHDiffPixTheta - // ON : from fHDiffPixThetaOFF - // OFF : from fHDiffPixThetaON - - TH3D *sp = NULL; - - if (fType == "MC") sp = fHDiffPixTheta; - else if (fType == "ON") sp = fHDiffPixThetaOFF; - else if (fType == "OFF") sp = fHDiffPixThetaON; - else - { - *fLog << err << "MPad::Process(); illegal data type... aborting" - << endl; - return kERROR; - } - - hdiff = sp->ProjectionZ("", binTheta, binTheta, - binPixel, binPixel, ""); - - if ( hdiff->GetEntries() == 0 ) - { - *fLog << warn << "MPad::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; + + hdiff = sp->ProjectionZ("", binTheta, binTheta, + binPixel, binPixel, ""); + Double_t integral = hdiff->Integral(); + // if there are no entries in hdiff, diff cannot be thrown + // in this case diff will be set to zero + if ( integral == 0 ) + { + //*fLog << warn << "MPad::Process(); fType = " << fType + // << ", projection of '" + // << sp->GetName() << "' for Theta bin " + // << binTheta << " and pixel " << j + // << " has no entries; set diff equal to the old difference " + // << endl; + + diff = TMath::Max(oldsigma2 - sigbarold2, + lambdabar*F2excess + oldsigma2 - sigmabar2); + + rc = 2; + fWarnings[rc]++; + } + + // start of else ------------------- + else + { + count = 0; + ok = kFALSE; + for (Int_t m=0; mGetRandom(); + + //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : " + // << m << ", " << diff_phot << endl; + + 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; - - rw = 1; - fWarnings[rw]++; - } - else - { - rw = 0; - fWarnings[rw]++; - } + // the following condition ensures that elNoise2Pix > 0.0 + if ( (diff + sigmabar2 - oldsigma2 + - lambdabar*F2excess) > 0.0 ) + { + ok = kTRUE; + break; + } + } + + if (!ok) + { + //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratioArea, lambdabar = " + // << theta << ", " + // << j << ", " << count << ", " << sigmabar2 << ", " + // << diff << ", " << oldsigma2 << ", " << ratioArea << ", " + // << lambdabar << endl; + diff = lambdabar*F2excess + oldsigma2 - sigmabar2; + + rw = 1; + fWarnings[rw]++; + } + } + // end of else -------------------- delete hdiff; @@ -2078,5 +2498,5 @@ // get the additional sigma^2 for this pixel (due to the padding) - addSig2 = sigma2*ratioArea - oldsigma2; + addSig2 = (sigma2 - oldsigma2) * ratioArea; addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton); @@ -2139,5 +2559,5 @@ - Double_t newsigma = sqrt( oldsigma2 + addSig2 ); + Double_t newsigma = sqrt( sigma2 * ratioArea ); Double_t newsigma_phot = newsigma / fPEperPhoton; ppix.SetRms( newsigma_phot ); @@ -2150,8 +2570,8 @@ fSigmabar->Calc(*fCam, *fPed, *fEvt); - *fLog << all << "MPad::Process(); after padding : " << endl; - fSigmabar->Print(""); - - *fLog << all << "Exit MPad::Process();" << endl; + //*fLog << all << "MPad::Process(); after padding : " << endl; + //fSigmabar->Print(""); + + //*fLog << all << "Exit MPad::Process();" << endl; rc = 0; @@ -2173,18 +2593,15 @@ *fLog << dec << setfill(' '); - int temp; - if (fWarnings[0] != 0.0) - temp = (int)(fWarnings[1]*100/fWarnings[0]); - else - temp = 100; + if (fWarnings[0] == 0 ) fWarnings[0] = 1; *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) - << temp - << "%) Warning : iteration to find acceptable sigma failed" + << (int)(fWarnings[1]*100/fWarnings[0]) + << "%) Pixel: iteration to find acceptable sigma failed" << ", fIter = " << fIter << endl; - *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) - << (int)(fErrors[1]*100/GetNumExecutions()) - << "%) Evts skipped due to: Sigmabar_old > 0" << endl; + *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3) + << (int)(fWarnings[2]*100/fWarnings[0]) + << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << endl; + *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) @@ -2192,19 +2609,7 @@ << "%) 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) @@ -2224,6 +2629,6 @@ //--------------------------------------------------------------- - TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 1500)); - c.Divide(3, 5); + TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 600)); + c.Divide(3, 2); gROOT->SetSelectedPad(NULL); @@ -2245,5 +2650,5 @@ //-------------------------------------------------------------------- - + /* c.cd(4); fHSigmaTheta->SetDirectory(NULL); @@ -2251,9 +2656,9 @@ fHSigmaTheta->DrawCopy(); fHSigmaTheta->SetBit(kCanDelete); - + */ //-------------------------------------------------------------------- - + /* c.cd(7); fHBlindPixNTheta->SetDirectory(NULL); @@ -2261,8 +2666,9 @@ fHBlindPixNTheta->DrawCopy(); fHBlindPixNTheta->SetBit(kCanDelete); + */ //-------------------------------------------------------------------- - + /* c.cd(10); fHBlindPixIdTheta->SetDirectory(NULL); @@ -2270,9 +2676,10 @@ fHBlindPixIdTheta->DrawCopy(); fHBlindPixIdTheta->SetBit(kCanDelete); - + */ //-------------------------------------------------------------------- // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2 + /* c.cd(5); TH2D *l1 = NULL; @@ -2296,8 +2703,10 @@ l2->DrawCopy("box"); l2->SetBit(kCanDelete);; + */ //-------------------------------------------------------------------- // draw the 3D histogram (target): Theta, pixel, Sigma + /* c.cd(6); TH2D *k1 = NULL; @@ -2321,13 +2730,14 @@ k2->DrawCopy("box"); k2->SetBit(kCanDelete);; - + */ //-------------------------------------------------------------------- - c.cd(13); + + c.cd(4); TH2D *m1; m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy"); m1->SetDirectory(NULL); - m1->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); + m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all \\Theta)"); m1->SetXTitle("Sigmabar-old"); m1->SetYTitle("Sigmabar-new"); @@ -2336,9 +2746,9 @@ m1->SetBit(kCanDelete);; - c.cd(14); + c.cd(5); TH2D *m2; m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy"); m2->SetDirectory(NULL); - m2->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); + m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all \\Theta)"); m2->SetXTitle("Sigmabar-old"); m2->SetYTitle("Sigmabar-new"); @@ -2347,9 +2757,9 @@ m2->SetBit(kCanDelete);; - c.cd(15); + c.cd(6); TH2D *m3; m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy"); m3->SetDirectory(NULL); - m3->SetTitle("(Target) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); + m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all \\Theta)"); m3->SetXTitle("Sigmabar-old"); m3->SetYTitle("Sigmabar-new"); @@ -2369,2 +2779,6 @@ + + + + Index: /trunk/MagicSoft/Mars/manalysis/MPad.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MPad.h (revision 4583) +++ /trunk/MagicSoft/Mars/manalysis/MPad.h (revision 4584) @@ -17,5 +17,5 @@ class MCerPhotEvt; class MPedPhotCam; -class MMcEvt; +class MPointingPos; class MSigmabar; class MParList; @@ -31,7 +31,7 @@ MCerPhotEvt *fEvt; MSigmabar *fSigmabar; - MMcEvt *fMcEvt; + MPointingPos *fPointPos; MPedPhotCam *fPed; - MBlindPixels *fBlinds; + MBlindPixels *fBlindPix; TString fType; // type of data to be padded @@ -41,5 +41,5 @@ Int_t fErrors[9]; - Int_t fWarnings[2]; + Int_t fWarnings[3]; //---------------------------------- @@ -89,7 +89,9 @@ - Bool_t Merge2Distributions(TH1D *hista, TH1D * histb, TH1D * histap, - TH2D *fHga, TH2D * fHgb, Int_t nbinssig); + Bool_t Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap, + TH2D *fHga, TH2D *fHgb, Int_t nbinssig); + Bool_t UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA, + Int_t nbinssig, Int_t l); public: @@ -98,16 +100,16 @@ Bool_t MergeONOFFMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon, - TH2D *sigthoff=NULL, TH3D *diffpixthoff=NULL,TH3D *sigmapixthoff=NULL, - TH2D *blindidthoff=NULL, TH2D *blindnthoff=NULL); + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D& diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon, + TH2D& sigthoff, TH3D& diffpixthoff,TH3D& sigmapixthoff, + TH2D& blindidthoff, TH2D& blindnthoff); Bool_t MergeONMC( - TH2D *sigthmc, TH3D *diffpixthmc, TH3D *sigmapixthmc, - TH2D *blindidthmc, TH2D *blindnthmc, - TH2D *sigthon, TH3D *diffpixthon, TH3D *sigmapixthon, - TH2D *blindidthon, TH2D *blindnthon); + TH2D& sigthmc, TH3D& diffpixthmc, TH3D& sigmapixthmc, + TH2D& blindidthmc, TH2D& blindnthmc, + TH2D& sigthon, TH3D&diffpixthon, TH3D& sigmapixthon, + TH2D& blindidthon, TH2D& blindnthon); Bool_t ReadPaddingDist(const char *filein); Index: /trunk/MagicSoft/Mars/manalysis/MSigmabar.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 4583) +++ /trunk/MagicSoft/Mars/manalysis/MSigmabar.cc (revision 4584) @@ -30,7 +30,10 @@ // // // This is the storage container to hold information about // -// "average" pedestal sigmas // -// // -// In calculating averages all sigmas are normalized to the area of pixel 0// +// "average" pedestal RMS // +// // +// The "average" pedestal RMS is calculated as // +// = sqrt( sum_i( (pedRMS_i)^2/area_i ) / no.of pixels ) // +// // +// which is the sqrt of the average (pedRMS^2 per area) // // // ///////////////////////////////////////////////////////////////////////////// @@ -109,5 +112,5 @@ // - // sum up sigma/sqrt(area) for each sector, + // sum up sigma^2/area for each sector, // separately for inner and outer region; // @@ -147,4 +150,6 @@ const Double_t ratio = geom.GetPixRatio(idx); + //*fLog << "pixel id, ratio = " << idx << ", " << ratio << endl; + const MGeomPix &gpix = geom[idx]; @@ -171,10 +176,10 @@ { innerPixels[sector]++; - innerSum[sector]+= sigma * sqrt(ratio); + innerSum[sector]+= sigma*sigma * ratio; } else { outerPixels[sector]++; - outerSum[sector]+= sigma * sqrt(ratio); + outerSum[sector]+= sigma*sigma * ratio; } } @@ -192,12 +197,12 @@ } - if (fInnerPixels > 0) fSigmabarInner = fSumInner / fInnerPixels; - if (fOuterPixels > 0) fSigmabarOuter = fSumOuter / fOuterPixels; - - // - // this is the average sigma/sqrt(area) (over all pixels) + if (fInnerPixels > 0) fSigmabarInner = sqrt(fSumInner / fInnerPixels); + if (fOuterPixels > 0) fSigmabarOuter = sqrt(fSumOuter / fOuterPixels); + + // + // this is the sqrt of the average sigma^2/area // fSigmabar = (fInnerPixels+fOuterPixels)<=0 ? 0: - (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels); + sqrt( (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels) ); for (UInt_t i=0; i<6; i++) @@ -209,7 +214,7 @@ const Double_t sum = ip + op; - fSigmabarInnerSector[i] = ip <=0 ? 0 : iss/ip; - fSigmabarOuterSector[i] = op <=0 ? 0 : oss/op; - fSigmabarSector[i] = sum<=0 ? 0 : (iss+oss)/sum; + fSigmabarInnerSector[i] = ip <=0 ? 0 : sqrt(iss/ip); + fSigmabarOuterSector[i] = op <=0 ? 0 : sqrt(oss/op); + fSigmabarSector[i] = sum<=0 ? 0 : sqrt((iss+oss)/sum); } @@ -217,5 +222,5 @@ //Print(opt); - return fSigmabar; + return fSigmabarInner; } Index: /trunk/MagicSoft/Mars/manalysis/MSigmabar.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 4583) +++ /trunk/MagicSoft/Mars/manalysis/MSigmabar.h (revision 4584) @@ -14,5 +14,5 @@ { private: - Float_t fSigmabar; // Sigmabar ("average" RMS) of pedestal + Float_t fSigmabar; // Sigmabar ( sqrt(average pedestalRMS^2) ) Float_t fSigmabarInner; // --only for inner pixels Float_t fSigmabarOuter; // --only for outer pixels Index: /trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 4583) +++ /trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc (revision 4584) @@ -37,5 +37,5 @@ // MPedPhotCam // MRawRunHeader -// MMcEvt (FIXME: Must be replaced by a 'real-data' container) +// MPointingPos // MCerPhotEvt // @@ -58,5 +58,5 @@ #include "MSigmabarParam.h" -#include "MMcEvt.hxx" +#include "MPointingPos.h" ClassImp(MSigmabarCalc); @@ -105,8 +105,8 @@ // This is needed for determining min/max Theta - fMcEvt = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt")); - if (!fMcEvt) + fPointPos = (MPointingPos*)pList->FindObject(AddSerialNumber("MPointingPos")); + if (!fPointPos) { - *fLog << warn << "MMcEvt not found... aborting." << endl; + *fLog << warn << "MPointingPos not found... aborting." << endl; fThetaMin = 0; fThetaMax = 90; @@ -148,11 +148,11 @@ if (rcGetTelescopeTheta()*kRad2Deg; + const Double_t theta = fPointPos->GetZd(); - if (theta>fSigmabarMax) fSigmabarMax=theta; - if (thetafThetaMax) fThetaMax=theta; + if (thetaFindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << dbginf << "MMcEvt not found... aborting." << endl; + fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointPos) + { + *fLog << dbginf << "MPointingPos not found... aborting." << endl; return kFALSE; } @@ -146,5 +146,5 @@ Int_t MFSelBasic::Process() { - const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta(); + const Double_t theta = fPointPos->GetZd(); fResult = kFALSE; Index: /trunk/MagicSoft/Mars/mfilter/MFSelBasic.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFSelBasic.h (revision 4583) +++ /trunk/MagicSoft/Mars/mfilter/MFSelBasic.h (revision 4584) @@ -14,5 +14,5 @@ #endif -class MMcEvt; +class MPointingPos; class MGeomCam; class MCerPhotEvt; @@ -23,5 +23,5 @@ { private: - const MMcEvt *fMcEvt; + const MPointingPos *fPointPos; const MGeomCam *fCam; // Camera Geometry const MCerPhotEvt *fEvt; // Cerenkov Photon Event Index: /trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mfilter/Makefile (revision 4583) +++ /trunk/MagicSoft/Mars/mfilter/Makefile (revision 4584) @@ -13,5 +13,5 @@ INCLUDES = -I. -I../mbase -I../mfbase -I../mraw -I../mmc -I../mdata \ -I../manalysis -I../mfileio -I../mgeom -I../mimage \ - -I../mhbase -I../mmain -I../mgui -I../msignal \ + -I../mhbase -I../mmain -I../mgui -I../msignal -I../mpointing \ -I../mpedestal Index: /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 4583) +++ /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc (revision 4584) @@ -26,7 +26,8 @@ ////////////////////////////////////////////////////////////////////////////// // // -// MHSigmaTheta (extension of Robert's MHSigmabarTheta) // +// MHSigmaTheta // // // -// calculates - the 2D-histogram sigmabar vs. Theta, and // +// calculates - the 2D-histogram sigmabar vs. Theta (Inner), and // +// - the 2D-histogram sigmabar vs. Theta (Outer) // // - the 3D-histogram sigma, pixel no., Theta // // - the 3D-histogram (sigma^2-sigmabar^2), pixel no., Theta // @@ -39,5 +40,5 @@ #include "MTime.h" -#include "MMcEvt.hxx" +#include "MPointingPos.h" #include "MBinning.h" @@ -46,4 +47,5 @@ #include "MGeomCam.h" +#include "MBlindPixels.h" #include "MPedPhotCam.h" @@ -70,8 +72,14 @@ fSigmaTheta.SetDirectory(NULL); - fSigmaTheta.SetName("2D-ThetaSigmabar"); + fSigmaTheta.SetName("2D-ThetaSigmabar(Inner)"); fSigmaTheta.SetTitle("2D: \\bar{\\sigma}, \\Theta"); fSigmaTheta.SetXTitle("\\Theta [\\circ]"); - fSigmaTheta.SetYTitle("Sigmabar"); + fSigmaTheta.SetYTitle("Sigmabar(Inner) / SQRT(Area)"); + + fSigmaThetaOuter.SetDirectory(NULL); + fSigmaThetaOuter.SetName("2D-ThetaSigmabar(Outer)"); + fSigmaThetaOuter.SetTitle("2D: \\bar{\\sigma}, \\Theta"); + fSigmaThetaOuter.SetXTitle("\\Theta [\\circ]"); + fSigmaThetaOuter.SetYTitle("Sigmabar(Outer) / SQRT(Area)"); fSigmaPixTheta.SetDirectory(NULL); @@ -80,5 +88,5 @@ fSigmaPixTheta.SetXTitle("\\Theta [\\circ]"); fSigmaPixTheta.SetYTitle("Pixel Id"); - fSigmaPixTheta.SetZTitle("\\sigma"); + fSigmaPixTheta.SetZTitle("Sigma"); fDiffPixTheta.SetDirectory(NULL); @@ -87,5 +95,5 @@ fDiffPixTheta.SetXTitle("\\Theta [\\circ]"); fDiffPixTheta.SetYTitle("Pixel Id"); - fDiffPixTheta.SetZTitle("{\\sigma}^{2}-\\bar{\\sigma}^{2}"); + fDiffPixTheta.SetZTitle("(Sigma2 - Sigmabar2)/Area"); // Set default binning @@ -102,7 +110,8 @@ binspix.SetEdges(578, -0.5, 577.5); - SetBinning(&fSigmaTheta, &binst, &binsb); - SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb); - SetBinning(&fDiffPixTheta, &binst, &binspix, &binsd); + SetBinning(&fSigmaTheta, &binst, &binsb); + SetBinning(&fSigmaThetaOuter, &binst, &binsb); + SetBinning(&fSigmaPixTheta, &binst, &binspix, &binsb); + SetBinning(&fDiffPixTheta, &binst, &binspix, &binsd); } @@ -120,7 +129,7 @@ } - fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); - if (!fMcEvt) - *fLog << warn << "MMcEvt not found... aborting." << endl; + fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); + if (!fPointPos) + *fLog << warn << "MPointingPos not found... aborting." << endl; @@ -132,4 +141,11 @@ } fPed->InitSize(fCam->GetNumPixels()); + + + fBlindPix = (MBlindPixels*)plist->FindObject("MBlindPixels"); + if (!fBlindPix) + { + *fLog << err << "MBlindPixels not found... continue. " << endl; + } @@ -177,6 +193,7 @@ if (binssigma) { - SetBinning(&fSigmaTheta, binstheta, binssigma); - SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma); + SetBinning(&fSigmaTheta, binstheta, binssigma); + SetBinning(&fSigmaThetaOuter, binstheta, binssigma); + SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma); } @@ -191,13 +208,18 @@ // Fill the histograms // +// ignore pixels if they are unused or blind +// Bool_t MHSigmaTheta::Fill(const MParContainer *par, const Stat_t w) { - Double_t theta = fMcEvt ? fMcEvt->GetTelescopeTheta()*kRad2Deg : 0; + Double_t theta = fPointPos->GetZd(); fSigmabar->Calc(*fCam, *fPed, *fEvt); - Double_t mysig = fSigmabar->GetSigmabarInner(); - - //*fLog << "theta, mysig = " << theta << ", " << mysig << endl; + Double_t mysig = fSigmabar->GetSigmabarInner(); + Double_t mysigouter = fSigmabar->GetSigmabarOuter(); + + //*fLog << "theta, mysig, mysigouter = " << theta << ", " << mysig + // << ", " << mysigouter << endl; fSigmaTheta.Fill(theta, mysig); + fSigmaThetaOuter.Fill(theta, mysigouter); const UInt_t npix = fEvt->GetNumPixels(); @@ -205,10 +227,24 @@ for (UInt_t i=0; i 0.5) + { + // inner pixel + diff = sigma*sigma*ratio - mysig*mysig; + } + else + { + // outer pixel + diff = sigma*sigma*ratio - mysigouter*mysigouter; + } + fDiffPixTheta.Fill(theta, (Double_t)id, diff); } @@ -262,5 +309,5 @@ AppendPad(""); - pad->Divide(3, 2); + pad->Divide(3, 3); // draw the 2D histogram Sigmabar versus Theta @@ -283,5 +330,5 @@ h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}"); + h->SetYTitle("(Sigma2 - Sigmabar2) / Area"); h->Draw("box"); h->SetBit(kCanDelete); @@ -293,5 +340,5 @@ h->SetTitle("\\sigma_{ped} vs. \\Theta (all pixels)"); h->SetXTitle("\\Theta [\\circ]"); - h->SetYTitle("\\sigma_{ped}"); + h->SetYTitle("Sigma"); h->Draw("box"); h->SetBit(kCanDelete); @@ -313,5 +360,5 @@ h->SetTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2} vs. pixel Id (all \\Theta)"); h->SetXTitle("Pixel Id"); - h->SetYTitle("\\sigma_{ped}^{2}-\\bar{\\sigma}_{ped}^{2}"); + h->SetYTitle("Sigma2 - sigmabar2) / Area"); h->Draw("box"); h->SetBit(kCanDelete); @@ -323,5 +370,5 @@ h->SetTitle("\\sigma_{ped} vs. pixel Id (all \\Theta)"); h->SetXTitle("Pixel Id"); - h->SetYTitle("\\sigma_{ped}"); + h->SetYTitle("Sigma"); h->Draw("box"); h->SetBit(kCanDelete); @@ -329,4 +376,7 @@ pad->cd(4); fSigmaTheta.Draw(opt); + + pad->cd(7); + fSigmaThetaOuter.Draw(opt); //pad->cd(8); Index: /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 4583) +++ /trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h (revision 4584) @@ -16,8 +16,9 @@ class MGeomCam; class MCerPhotEvt; -class MMcEvt; +class MPointingPos; class MPedPhotCam; class MSigmabar; class MParList; +class MBlindPixels; @@ -29,8 +30,11 @@ MCerPhotEvt *fEvt; //! MSigmabar *fSigmabar; //! - MMcEvt *fMcEvt; //! - - TH2D fSigmaTheta; // 2D-distribution sigmabar versus Theta; - // sigmabar is the average pedestasl sigma in an event + MPointingPos *fPointPos; //! + MBlindPixels *fBlindPix; //! + + // sigmabar is the average pedestal sigma + TH2D fSigmaTheta; // 2D-distribution sigmabar versus Theta (Inner) + TH2D fSigmaThetaOuter; // 2D-distribution sigmabar versus Theta (Outer) + TH3D fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma TH3D fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2 @@ -46,4 +50,7 @@ const TH2D *GetSigmaTheta() { return &fSigmaTheta; } const TH2D *GetSigmaTheta() const { return &fSigmaTheta; } + + const TH2D *GetSigmaThetaOuter() { return &fSigmaThetaOuter; } + const TH2D *GetSigmaThetaOuter() const { return &fSigmaThetaOuter; } const TH3D *GetSigmaPixTheta() { return &fSigmaPixTheta; }