Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.cc (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.cc (revision 6232) @@ -645,7 +645,11 @@ // The default is not cut, i.e. all values (0-1) are taken - fThetaMin = 0; // in miliradians // FIXME: change name - fThetaMax = 90; // new: in deg; old: (1570) in miliradians // FIXME: change name - + //fThetaMin = 0; // in miliradians // FIXME: change name + //fThetaMax = 1570; // in miliradians // FIXME: change name + + fThetaMin = 0; // degrees (new in magic) + fThetaMax = 90; // degrees (new in magic) + + fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1 @@ -770,14 +774,14 @@ // Angular cuts not yet implemented ... - Double_t ConvMMToDeg = 0.00337034; - - fDistCutLow = 0.4/ConvMMToDeg; - fDistCutUp = 1.5/ConvMMToDeg; - - fLengthCutLow = 0.1/ConvMMToDeg; - fLengthCutUp = 1/ConvMMToDeg; - - fWidthCutLow = 0.07/ConvMMToDeg; - fWidthCutUp = 1/ConvMMToDeg; + fConvMMToDeg = 0.00337034; + + fDistCutLow = 0.4/fConvMMToDeg; + fDistCutUp = 1.5/fConvMMToDeg; + + fLengthCutLow = 0.1/fConvMMToDeg; + fLengthCutUp = 1/fConvMMToDeg; + + fWidthCutLow = 0.07/fConvMMToDeg; + fWidthCutUp = 1/fConvMMToDeg; // ENDTMP @@ -1141,6 +1145,12 @@ // in previous functions... - fThetaMin = int(ThetaMin*1000.0); - fThetaMax = int(ThetaMax*1000.0); + // fThetaMin = int(ThetaMin*1000.0); + // fThetaMax = int(ThetaMax*1000.0); + + + // in new magic theta is given in deg + // 0.5 is added so that the rounding to integer is correct + fThetaMin = int(ThetaMin + 0.5); + fThetaMax = int(ThetaMax + 0.5); fThetaRangeString = ("ThetaRange"); @@ -1148,5 +1158,7 @@ fThetaRangeString += ("_"); fThetaRangeString += (fThetaMax); - fThetaRangeString += ("mRad"); + // fThetaRangeString += ("mRad"); + fThetaRangeString += ("_Degrees"); // new in magic + return kTRUE; @@ -1168,4 +1180,21 @@ + + return kTRUE; +} + + +Bool_t MFindSupercutsONOFF::SetFilters(Double_t LeakageMax, Double_t DistMax, Double_t DistMin) +{ + + fDistCutLow = DistMin/fConvMMToDeg; + fDistCutUp = DistMax/fConvMMToDeg; + fLeakageMax = LeakageMax; + + *fLog << "MFindSupercutsONOFF::SetSizeRange" << endl + << "Data matrices will be filled with events whose MHillasSrc.fDist " << endl + << "is in the range " + << fDistCutLow <<"-"<< fDistCutUp << " degrees" << endl + << "and fLeakage " << "< " << fLeakageMax << endl; return kTRUE; @@ -1682,5 +1711,5 @@ - // TMP + // Cuts in Size, @@ -1696,8 +1725,21 @@ SizeCutMax.SetInverted(); - - - - // ENDTMP + // Cuts in Dist + TString DistCutMinString ("MHillasSrc.fDist>"); + DistCutMinString += fDistCutLow; + MContinue DistCutMin(DistCutMinString); + DistCutMin.SetInverted(); + + TString DistCutMaxString ("MHillasSrc.fDist<"); + DistCutMaxString += fDistCutUp; + MContinue DistCutMax(DistCutMaxString); + DistCutMax.SetInverted(); + + + // Cuts in leakage + TString LeakCutMaxString ("MNewImagePar.fLeakage1<"); + LeakCutMaxString += fLeakageMax; + MContinue LeakCutMax(LeakCutMaxString); + LeakCutMax.SetInverted(); @@ -1706,4 +1748,7 @@ Double_t prob = whichfractiontrain /(whichfractiontrain+whichfractiontest); + + + MFRandomSplit split(prob); @@ -1735,12 +1780,16 @@ tlist.AddToList(&read); -// tlist.AddToList(&ThetaCutMin); -// tlist.AddToList(&ThetaCutMax); - - //TMP + tlist.AddToList(&ThetaCutMin); + tlist.AddToList(&ThetaCutMax); + + tlist.AddToList(&SizeCutMin); tlist.AddToList(&SizeCutMax); - //ENDTMP + tlist.AddToList(&DistCutMin); + tlist.AddToList(&DistCutMax); + tlist.AddToList(&LeakCutMax); + + tlist.AddToList(&split); tlist.AddToList(&filltrain); @@ -1848,5 +1897,5 @@ - // TMP + // Cuts in Size, @@ -1863,7 +1912,27 @@ - - - // ENDTMP + // Cuts in Dist + TString DistCutMinString ("MHillasSrc.fDist>"); + DistCutMinString += fDistCutLow; + MContinue DistCutMin(DistCutMinString); + DistCutMin.SetInverted(); + + TString DistCutMaxString ("MHillasSrc.fDist<"); + DistCutMaxString += fDistCutUp; + MContinue DistCutMax(DistCutMaxString); + DistCutMax.SetInverted(); + + + // Cuts in leakage + TString LeakCutMaxString ("MNewImagePar.fLeakage1<"); + LeakCutMaxString += fLeakageMax; + MContinue LeakCutMax(LeakCutMaxString); + LeakCutMax.SetInverted(); + + + + + + @@ -1871,5 +1940,7 @@ Double_t prob = whichfractiontrain/(whichfractiontrain + whichfractiontest); - + + + MFRandomSplit split(prob); @@ -1903,13 +1974,17 @@ tlist.AddToList(&read); -// tlist.AddToList(&ThetaCutMin); -// tlist.AddToList(&ThetaCutMax); + tlist.AddToList(&ThetaCutMin); + tlist.AddToList(&ThetaCutMax); - //TMP + tlist.AddToList(&SizeCutMin); tlist.AddToList(&SizeCutMax); - //ENDTMP + tlist.AddToList(&DistCutMin); + tlist.AddToList(&DistCutMax); + tlist.AddToList(&LeakCutMax); + + tlist.AddToList(&split); @@ -1975,4 +2050,82 @@ + +/// **********/// + +// -------------------------------------------------------------------------- +// +// Read only training matrices ON and OFF +// +// + +Bool_t MFindSupercutsONOFF::ReadMatrixTrain(const TString &filetrainON, const TString &filetrainOFF) +{ + //-------------------------- + // read in training matrix ON sample + + TFile filetr(filetrainON); + fMatrixTrain->Read("MatrixTrain"); + fMatrixTrain->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrixTrain; Training matrix for ON sample was read in from file '" + << filetrainON << "'" << endl; + filetr.Close(); + + + // read in training matrix OFF sample + + TFile filetrOFF(filetrainOFF); + fMatrixTrainOFF->Read("MatrixTrainOFF"); + fMatrixTrainOFF->Print("SizeCols"); + + *fLog << "MFindSupercutsONOFF::ReadMatrixTrain; Training matrix for OFF sample was read in from file '" + << filetrainOFF << "'" << endl; + filetrOFF.Close(); + + return kTRUE; + +} + +// -------------------------------------------------------------------------- +// +// Read only test matrices ON and OFF +// +// + +Bool_t MFindSupercutsONOFF::ReadMatrixTest(const TString &filetestON, const TString &filetestOFF) +{ + //-------------------------- + // read in testing matrix ON sample + +//-------------------------- + // read in test matrix for ON sample + + TFile filete(filetestON); + fMatrixTest->Read("MatrixTest"); + fMatrixTest->Print("SizeCols"); + + *fLog << "MFindSupercuts::ReadMatrixTest; Test matrix for ON sample was read in from file '" + << filetestON << "'" << endl; + filete.Close(); + + + //-------------------------- + // read in test matrix for OFF sample + + TFile fileteOFF(filetestOFF); + fMatrixTestOFF->Read("MatrixTestOFF"); + fMatrixTestOFF->Print("SizeCols"); + + *fLog << "MFindSupercutsONOFF::ReadMatrixTest; Test matrix for OFF sample was read in from file '" + << filetestOFF << "'" << endl; + filete.Close(); + + + + return kTRUE; + +} + + /// **********/// @@ -1999,5 +2152,5 @@ //-------------------------- // read in test matrix - + TFile filete(filetest); fMatrixTest->Read("MatrixTest"); @@ -2007,5 +2160,5 @@ << filetest << "'" << endl; filete.Close(); - + return kTRUE; } @@ -2032,5 +2185,5 @@ //-------------------------- // read in test matrix - + TFile filete(filetest); fMatrixTestOFF->Read("MatrixTestOFF"); @@ -2040,5 +2193,5 @@ << filetest << "'" << endl; filete.Close(); - + return kTRUE; } Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.h (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFF.h (revision 6232) @@ -59,4 +59,7 @@ Double_t fWidthCutUp; + Double_t fLeakageMax; + + Double_t fConvMMToDeg; // ENDTMP @@ -85,6 +88,6 @@ Double_t fAlphaBkgMax; - Int_t fThetaMin; // Cuts in ThetaOrig.fVal (in mili radians!!!) - Int_t fThetaMax; // Cuts in ThetaOrig.fVal (in mili radians !!!) + Int_t fThetaMin; // Cut in Theta value (in degrees) + Int_t fThetaMax; // Cuts in Theta value (in degrees) TString fThetaRangeString; @@ -314,4 +317,8 @@ Bool_t CheckAlphaSigBkg(); + + + + void SetUseFittedQuantities (Bool_t b) {fUseFittedQuantities = b;} @@ -371,4 +378,10 @@ Bool_t ReadMatrixOFF( const TString &filetrainOFF, const TString &filetestOFF); + + Bool_t ReadMatrixTrain(const TString &filetrainON, const TString &filetrainOFF); + Bool_t ReadMatrixTest(const TString &filetestON, const TString &filetestOFF); + + + Bool_t ComputeNormFactorTrain (); Bool_t ComputeNormFactorTest (); @@ -395,4 +408,7 @@ Bool_t SetSizeRange(Double_t SizeMin, Double_t SizeMax); + + Bool_t SetFilters(Double_t LeakageMax, Double_t DistMax, Double_t DistMin); + Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.cc (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.cc (revision 6232) @@ -110,11 +110,21 @@ //--------------------------- - + // ********************************* + // ** NOT TRUE ANY LONGER... ** // // Cuts (low and up) in variable ThetaOrig.fVal // The default is not cut, i.e. all values (0-1) are taken // fThetaOrig.fVal is measured in Radians; thus 1 = 57 degrees. + // fThetaMin = 0.0; + // fThetaMax = 1.0; + // ********************************* + + + // For MAGIC new (July 2004) Theta is measured in degrees + fThetaMin = 0.0; - fThetaMax = 1.0; + fThetaMax = 90.0; + + fAlphaSig = 20; // By default, signal is expected in alpha<20 for CT1 @@ -545,18 +555,5 @@ - // Give a "soft warning" if Vector has increasing values - // of CosThetas; which means decreasing values of Thetas - /* - if (fCosThetaRangeVector[0] <= fCosThetaRangeVector[1]) - { - *fLog << "MFindSupercutsONOFFThetaLoop::SetCosThetaRangeVector; " - << endl - << "Values in vector fCosThetaRangeVector are in increasing " - << "order, i.e., values in theta will be in decreasing order" - << endl - << "DO NOT forget it !!" << endl; - } - - */ + // fCosThetaBinCenterVector is defined @@ -622,8 +619,18 @@ // Now costheta values have to be converted to Radians // Function TMath::ACos(x) gives angle in Radians - + fThetaMin = TMath::ACos(fThetaMin); fThetaMax = TMath::ACos(fThetaMax); + + + // For new MAGIC (July 2004), angles are given in degrees. + // So I convert angles into degrees by multiplying by a the + // conversion rad/degree + + fThetaMin = fThetaMin * TMath::RadToDeg(); + fThetaMax = fThetaMax * TMath::RadToDeg(); + + // If fThetaMin > fThetaMax means that fCosThetaRangeVector @@ -644,5 +651,5 @@ << "MFindSupercutsONOFFThetaLoop::SetThetaRange; " << endl << "Theta range for this iteration is set to " - << fThetaMin << "-" << fThetaMax << " Radians." << endl + << fThetaMin << "-" << fThetaMax << " degrees." << endl << "-------------------------------------------------"<< endl; @@ -932,19 +939,5 @@ - // Histograms references are removed from the current directory - - /* IT DOES NOT WORK !!! THESE STATEMENTS PRODUCE SEGMENTATION FAULT. - fNormFactorTrainHist -> SetDirectory(NULL); - - fNormFactorTestHist -> SetDirectory(NULL); - - fSigmaLiMaTrainHist -> SetDirectory(NULL); - - fSigmaLiMaTestHist -> SetDirectory(NULL); - - fSigmaLiMaTestHist -> SetDirectory(NULL); - - fNexTestHist -> SetDirectory(NULL); - */ + @@ -1059,4 +1052,43 @@ Double_t fofftest) { + + // Check that fractions set by user are within the expected limits + + if (fontrain < 0.00 || fontrain > 1.00 || + fontest < 0.00 || fontest > 1.00 || + fofftrain < 0.00 || fofftrain > 1.00 || + fofftest < 0.00 || fontest > 1.00) + { + *fLog << err + << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: " + << endl + << "Fractions of ON/OFF data for TRAIN/TEST sample which were specified by user " << endl + << "are outside the natural limits 0.0-1.0. " << endl + << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted." + << endl; + + exit(1); + + + } + + + if((fontrain+fontest) > 1.00 || + (fofftrain+fofftest) >1.00) + + { + *fLog << err + << "MFindSupercutsONOFFThetaLoop::SetFractionTrainTestOnOffEvents: " + << endl + << "The sum of the fractions of TRAIN/TEST for ON or OFF data which were specified by user " << endl + << "are larger than 1.0. " << endl + << "Set reasonable fraction numbers for TRAIN/TEST please. The program will be aborted." + << endl; + + exit(1); + + + } + fWhichFractionTrain = fontrain; @@ -1154,5 +1186,5 @@ // Path + "OptSCParametersONOFF" + "ThetaRange" - // + int {(fThetaMin_fThetaMax)*1000} + ".root" + // + int {(fThetaMin_fThetaMax)} + ".root" @@ -1160,5 +1192,5 @@ // Path + "ON(OFF)Data" + "Train(Test)" + "Matrix" + Fraction - // + "ThetaRange" + int {(fThetaMin_fThetaMax)*1000} + ".root" + // + "ThetaRange" + int {(fThetaMin_fThetaMax)} + ".root" @@ -1176,6 +1208,13 @@ - Int_t ThetaMinTimes1000 = 0; - Int_t ThetaMaxTimes1000 = 0; + // Int_t ThetaMinTimes1000 = 0; + // Int_t ThetaMaxTimes1000 = 0; + + // For MAGIC new (July 2004) Theta is given in deg + + Int_t ThetaMinInt = 0; + Int_t ThetaMaxInt = 0; + + Int_t FractionONTrainTimes100 = 0; Int_t FractionONTestTimes100 = 0; @@ -1200,8 +1239,19 @@ - ThetaMinTimes1000 = int(fThetaMin*1000); - ThetaMaxTimes1000 = int(fThetaMax*1000); - sprintf(ThetaRangeString, "ThetaRange%d_%dmRad", - ThetaMinTimes1000, ThetaMaxTimes1000); + //ThetaMinTimes1000 = int(fThetaMin*1000); + //ThetaMaxTimes1000 = int(fThetaMax*1000); + //sprintf(ThetaRangeString, "ThetaRange%d_%dmRad", + // ThetaMinTimes1000, ThetaMaxTimes1000); + + // For New MAGIC (July 2004) + // Theta is given in deg, no need to multiply by 1000 + + // 0.5 is added so that the rounding to integer is correct + ThetaMinInt = int(fThetaMin + 0.5); + ThetaMaxInt = int(fThetaMax + 0.5); + sprintf(ThetaRangeString, "ThetaRange%d_%d_Degrees", + ThetaMinInt, ThetaMaxInt); + + fThetaRangeStringVector[i] = (ThetaRangeString); @@ -1211,24 +1261,12 @@ // endtmp - Double_t tmpdouble = 0.0; - tmpdouble = fWhichFractionTrain*100 - int(fWhichFractionTrain*100); - - FractionONTrainTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTrain*100): - int(fWhichFractionTrain*100) + 1; - - - tmpdouble = fWhichFractionTest*100 - int(fWhichFractionTest*100); - - FractionONTestTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTest*100): - int(fWhichFractionTest*100) + 1; - - tmpdouble = fWhichFractionTrainOFF*100 - int(fWhichFractionTrainOFF*100); - - FractionOFFTrainTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTrainOFF*100): - int(fWhichFractionTrainOFF*100) + 1; - - tmpdouble = fWhichFractionTestOFF*100 - int(fWhichFractionTestOFF*100); - FractionOFFTestTimes100 = tmpdouble < 0.5 ? int(fWhichFractionTestOFF*100): - int(fWhichFractionTestOFF*100) + 1; + // 0.5 is added so that the rounding to integer is correct + + FractionONTrainTimes100 = int(fWhichFractionTrain*100 + 0.5); + FractionONTestTimes100 = int(fWhichFractionTest*100 + 0.5); + FractionOFFTrainTimes100 = int(fWhichFractionTrainOFF*100 + 0.5); + FractionOFFTestTimes100 = int(fWhichFractionTestOFF*100 + 0.5); + + @@ -1525,4 +1563,5 @@ FindSupercuts.SetSizeRange(fSizeCutLow, fSizeCutUp); + FindSupercuts.SetFilters(fLeakMax, fDistMax, fDistMin); @@ -1627,8 +1666,19 @@ << "Reading Data matrices from root files... " << endl; + /* FindSupercuts.ReadMatrix(fTrainMatrixONFilenameVector[i], fTestMatrixONFilenameVector[i]); FindSupercuts.ReadMatrixOFF(fTrainMatrixOFFFilenameVector[i], fTestMatrixOFFFilenameVector[i]); + */ + + FindSupercuts.ReadMatrixTrain(fTrainMatrixONFilenameVector[i], + fTrainMatrixOFFFilenameVector[i]); + + if (fTestParameters) + { + FindSupercuts.ReadMatrixTest(fTestMatrixONFilenameVector[i], + fTestMatrixOFFFilenameVector[i]); + } } @@ -2325,6 +2375,6 @@ // TEMP - // cout << NexVSAlphaNameTmp << endl; - // cout << SignificanceVSAlphaNameTmp << endl; + //cout << NexVSAlphaNameTmp << endl; + //cout << SignificanceVSAlphaNameTmp << endl; // ENDTEMP @@ -2388,4 +2438,7 @@ AlphaOFFAfterCuts = (TH1F*) rootfile.Get(TmpAlphaOFFHistoName); + + + @@ -2885,4 +2938,18 @@ { + + *fLog << endl << endl + << "****************************************************************************" << endl + << "Combining number of excess events and significance for the " << endl + << "selected bins in theta angle using the function " << endl + << " MFindSupercutsONOFFThetaLoop::ComputeOverallSignificance() " << endl + << "****************************************************************************" << endl + << endl; + + + + + + // First of all let's check that we have the proper the "drinks" for the party @@ -2917,17 +2984,23 @@ // Remember that histogram index START AT 1 !!! + + + // Histograms that will contain the alpha distributions + // for the combined theta bins + + TH1F* CombinedAlphaTrainON = NULL; - TString CombinedAlphaTrainONName = ("CombinedAlphaTrainON"); + TString CombinedAlphaTrainONName = ("TrainSampleAlphaONCombinedThetaBins"); TH1F* CombinedAlphaTrainOFF = NULL; - TString CombinedAlphaTrainOFFName = ("CombinedAlphaTrainOFF"); + TString CombinedAlphaTrainOFFName = ("TrainSampleAlphaOFFCombinedThetaBins"); TH1F* CombinedAlphaTestON = NULL; - TString CombinedAlphaTestONName = ("CombinedAlphaTestON"); + TString CombinedAlphaTestONName = ("TestSampleAlphaONCombinedThetaBins"); TH1F* CombinedAlphaTestOFF = NULL; - TString CombinedAlphaTestOFFName = ("CombinedAlphaTestOFF"); + TString CombinedAlphaTestOFFName = ("TestSampleAlphaONCombinedThetaBins"); - - + + @@ -3089,4 +3162,12 @@ { + + + *fLog << endl + << "****************************************************************************" << endl + << "Combining data for the TRAIN sample " << endl + << "****************************************************************************" << endl + << endl; + // Retrieving NormFactorTrainHisto from root file. @@ -3489,9 +3570,13 @@ CombinedAlphaTrainOFF->SetEntries(tmpint); - + // Check that the number of events in the combined normalized histogram is + // that of the initial sample + + /* cout << "SILLY INFO FOR TRAIN SAMPLE: Total number of events Before nomralization and " << " re-re normalized" << endl << TotalNumberOFFTrain << " - " << EventCounter << endl; + */ } @@ -3506,4 +3591,11 @@ if (CombineTestData) { + + *fLog << endl + << "****************************************************************************" << endl + << "Combining data for the TEST sample " << endl + << "****************************************************************************" << endl + << endl; + // Retrieving NormFactorTestHisto from root file. @@ -3930,7 +4022,14 @@ + // Check that the number of events in the combined normalized histogram is + // that of the initial sample + + /* cout << "SILLY INFO FOR TEST SAMPLE: Total number of events Before nomralization and " << " re-re normalized" << endl << TotalNumberOFFTest << " - " << EventCounter << endl; + + */ + /* @@ -3982,4 +4081,12 @@ if (SuccessfulThetaBinsVector.GetSize() >= 1) { + + + *fLog << endl + << "**************************************************************************************" << endl + << "Computation of excess events and singnificance for the combined alpha histograms " << endl + << "**************************************************************************************" << endl + << endl; + // There is, at least, ONE theta bin for which optimization of // supercuts was SUCCESSFUL, and thus, it is worth to @@ -4026,5 +4133,7 @@ // Name of psfile where Alpha plot will be stored - TString psfilename = ("TRAINSampleCombined"); + + TString psfilename = (fPathForFiles); + psfilename += ("TRAINSampleCombinedThetaBins"); findsig.FindSigmaONOFF(CombinedAlphaTrainON, CombinedAlphaTrainOFF, @@ -4053,5 +4162,6 @@ // Name of psfile where Alpha plot will be stored - TString psfilename = ("TESTSampleCombined"); + TString psfilename = (fPathForFiles); + psfilename += ("TESTSampleCombinedThetaBins"); findsig.FindSigmaONOFF(CombinedAlphaTestON, CombinedAlphaTestOFF, @@ -4075,4 +4185,77 @@ + *fLog << endl + << "**************************************************************************************" << endl + << "Alpha distributions (Train/Test and ON/OFF) for the combined " << endl + << " theta bins are stored in the root file " << endl + << fAlphaDistributionsRootFilename << endl + << "**************************************************************************************" << endl + << endl; + + + + TFile rootfiletowrite(fAlphaDistributionsRootFilename, "UPDATE"); + + + + // Histograms that will contain the Norm factors for the + // combined TRAIN/TEST theta bins are created + + + TString CombinedNormFactorTrainName = ("TrainSampleNormFactorCombinedThetaBins"); + TH1F CombinedNormFactorTrain (CombinedNormFactorTrainName, + CombinedNormFactorTrainName, + 1, 0.0, 1.0); + + TString CombinedNormFactorTestName = ("TestSampleNormFactorCombinedThetaBins"); + TH1F CombinedNormFactorTest (CombinedNormFactorTestName, + CombinedNormFactorTestName, + 1, 0.0, 1.0); + + + if (CombineTrainData) + { + + // Histogram that will contain the Norm factors for the + // combined TRAIN theta bins are filled + + + CombinedNormFactorTrain.Fill(0.5, MeanNormFactorTrain); + + + // Train histos are written to root file + + CombinedAlphaTrainON -> Write(); + CombinedAlphaTrainOFF -> Write(); + CombinedNormFactorTrain.Write(); + + + + } + + + if (CombineTestData) + { + + // Histogram that will contain the Norm factors for the + // combined TEST theta bins are filled + + + CombinedNormFactorTest.Fill(0.5, MeanNormFactorTest); + + // Test histos are written to root file + + CombinedAlphaTestON -> Write(); + CombinedAlphaTestOFF -> Write(); + CombinedNormFactorTest.Write(); + + + } + + + // Root file is closed + + rootfiletowrite.Close(); + } Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.h (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MFindSupercutsONOFFThetaLoop.h (revision 6232) @@ -73,4 +73,7 @@ Double_t fSizeCutUp; + Double_t fLeakMax; + Double_t fDistMax; + Double_t fDistMin; // Variables for binning of alpha plots @@ -373,4 +376,8 @@ void SetSizeRange(Double_t SizeMin, Double_t SizeMax) {fSizeCutLow = SizeMin; fSizeCutUp = SizeMax; } + +//DM: + void SetFilters(Double_t LeakageMax, Double_t DistMax, Double_t DistMin) + {fLeakMax = LeakageMax; fDistMax = DistMax; fDistMin = DistMin; } // Function that loops over the theta ranges defined by Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MHFindSignificanceONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MHFindSignificanceONOFF.cc (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MHFindSignificanceONOFF.cc (revision 6232) @@ -132,6 +132,5 @@ //----------------------------- - // ignore unwanted points - if (error > 1.e19) + // ignore unwanted points if (error > 1.e19) continue; @@ -373,5 +372,7 @@ // allow rebinning of the alpha plot - fRebin = kTRUE; +// fRebin = kTRUE; + + fRebin = kFALSE; // allow reducing the degree of the polynomial @@ -396,4 +397,5 @@ fPrintResultsOntoAlphaPlot = kTRUE; + //fPrintResultsOntoAlphaPlot = kFALSE; @@ -872,6 +874,4 @@ if (fDraw) { - - // TEMPORALLY I will plot fHistOFF and fHistOFFNormalized (with the fits) @@ -1046,13 +1046,10 @@ // count bins with zero entry - if (content <= 0.0) - { - fNzeroOFF++; - // The error of the bin is set to a huge number, - // so that it does not have any weight in the fit - fHistOFF->SetBinError(i, dummy); - } + if (content < 0.0001) + fNzeroOFF++; + // set minimum error - if (content < 9.0) + // modified by K. Mase + if (content < 0.0) { fMlowOFF += 1; @@ -1081,7 +1078,8 @@ mean /= ((Double_t) fMbinsOFF); rms /= ((Double_t) fMbinsOFF); - } - - rms = sqrt( rms - mean*mean ); + rms = sqrt( rms - mean*mean ); + } + + // if there are no events in the background region @@ -1155,20 +1153,5 @@ fConstantBackg = kFALSE; - - // Condition for disabling the fitting procedure and - // assuming a constant background (before Nov 2004) - - // if ( fNzeroOFF > 0 || (Double_t)fMlowOFF>0.05*(Double_t)fMbinsOFF ) - - - // Condition for disabling the fitting procedure and - // assuming a constant background (After Nov 01 2004) - // I softened the condition to allow the fit also in situations - // where the reduction of the background is such that very - // few events survived; which is - // Specially frequent with Random Forest at high Sizes) - - if ( (Double_t)fNzeroOFF > 0.1*(Double_t)fMbinsOFF || - (Double_t)fMlowOFF > 0.2*(Double_t)fMbinsOFF ) + if ( fNzeroOFF > 0 || (Double_t)fMlowOFF>0.05*(Double_t)fMbinsOFF ) { *fLog << "MHFindSignificanceONOFF::FitPolynomialOFF; polynomial fit not possible, fNzeroOFF, fMlowOFF, fMbinsOFF = " @@ -1201,5 +1184,6 @@ Double_t val, err; val = mean; - err = rms; // sqrt( mean / (Double_t)fMbinsOFF ); + err = rms; + fPolyOFF->SetParameter(0, val); @@ -1571,18 +1555,13 @@ // count bins with zero entry - if (content <= 0.0) - { - fNzero++; - // The error of the bin is set to a huge number, - // so that it does not have any weight in the fit - fHistOFF->SetBinError(i, dummy); - } - + /*if (content <= 0.0) + fNzero++;*/ + // set minimum error - if (content < 9.0) + /*if (content < 9.0) { fMlow += 1; fHist->SetBinError(i, 3.0); - } + }*/ //*fLog << "Take : i, content, error = " << i << ", " @@ -1681,21 +1660,5 @@ fConstantBackg = kFALSE; - - // Condition for disabling the fitting procedure and - // assuming a constant background (before Nov 2004) - - // if ( fNzero > 0 || (Double_t)fMlow>0.05*(Double_t)fMbins ) - - - // Condition for disabling the fitting procedure and - // assuming a constant background (After Nov 01 2004) - // I softened the condition to allow the fit also in situations - // where the reduction of the background is such that very - // few events survived; which is - // Specially frequent with Random Forest at high Sizes) - - if ( (Double_t)fNzero > 0.1*(Double_t)fMbins || - (Double_t)fMlow > 0.2*(Double_t)fMbins ) - + if ( fNzero > 0 || (Double_t)fMlow>0.05*(Double_t)fMbins ) { *fLog << "MHFindSignificanceONOFF::FitPolynomial; polynomial fit not possible, fNzero, fMlow, fMbins = " @@ -1728,5 +1691,5 @@ Double_t val, err; val = mean; - err = rms; // sqrt( mean / (Double_t)fMbins ); + err = sqrt( mean / (Double_t)fMbins ); fPoly->SetParameter(0, val); @@ -2776,5 +2739,5 @@ } sum *= fac*fac; - + if (sum < 0.0) { @@ -3231,7 +3194,5 @@ CanvasHistOFFNormalized -> Update(); - - - */ +*/ return kTRUE; @@ -3292,4 +3253,8 @@ fpoly->DrawCopy("same"); } + + fHistOFFNormalized->SetLineColor(kRed); + fHistOFFNormalized->DrawCopy("same"); + } @@ -3316,8 +3281,8 @@ { // 10, 1 is white and solid - fbackg->SetLineColor(10); + fbackg->SetLineColor(kRed); fbackg->SetLineStyle(1); fbackg->SetLineWidth(4); - // fbackg->DrawCopy("same"); I do not want to draw it... already too many things. + //fbackg->DrawCopy("same"); // I do not want to draw it... already too many things. } } @@ -3495,5 +3460,5 @@ // fPsFilenameString - filename += ("AlphaPlotAfterSupercuts.ps"); + filename += ("AlphaPlot.ps"); cout << filename << endl; fCanvas -> SaveAs(filename); Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.cc (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.cc (revision 6232) @@ -133,4 +133,7 @@ fWidthUpperLimit = 0.4; + // Temp upper limit to get rid of triangular events + fLog10ConcUpperLimit = -0.35; + fLeakage1UpperLimit = 0.25; @@ -138,6 +141,6 @@ fDistLowerLimit = 0.1; - fLengthLowerLimit = 0.09; - fWidthLowerLimit = 0.05; + fLengthLowerLimit = 0.05; // values at October 21th 2004 + fWidthLowerLimit = 0.03; // values at October 21th 2004 @@ -237,4 +240,81 @@ } + + +// Function implementing a filter for muon induced images +// (Values from Keichi cuts October 18th 2004) + +// Function returns kTRUE if the event does NOT survive the +// muon filter. So kTRUE means event MUST be rejected + + +Bool_t MSupercutsCalcONOFF::MuonFilter (Double_t Size, + Double_t Length) +{ + // Parametrization for cut in length. + // Event is rejected if + // Log10(Length) < k1*(log10(Size)) + k2 (Size<10000photons) + // k1 = 0.286; k2 = -1.91 (October 18th 2004) + + // Length < 0.2 deg (Size > 10000) + + + // Parametrization for cut in Length/Size + // Event is rejected if + // log10(LenOverSize) < k3*(log10(Size)) + k4 (Size<4000photons) + // k3 = -0.780; k4 = -1.71 + + + Double_t LengthUpperLimit = 0.0; + Double_t LengthOverSizeUpperLimit = 0.0; + + + Double_t SizeUpperLimitForCutInLength = 10000; // in photons + Double_t SizeUpperLimitForCutInLengthOverSize = 4000; // in photons + + Double_t k1 = 0.286; + Double_t k2 = -1.91; + Double_t k3 = -0.78; + Double_t k4 = -1.71; + + + + + + if (Size <= SizeUpperLimitForCutInLength) + { + LengthUpperLimit = k1 * TMath::Log10(Size) + k2; + LengthUpperLimit = TMath::Power(10, LengthUpperLimit); + } + else + { + LengthUpperLimit = 0.2; + } + + // Apply cut in Length + + if (Length < LengthUpperLimit) + return kTRUE; + + + + if (Size < SizeUpperLimitForCutInLengthOverSize) + { + LengthOverSizeUpperLimit = k3 * TMath::Log10(Size) + k4; + LengthOverSizeUpperLimit = TMath::Power(10, LengthOverSizeUpperLimit); + + if (Length/Size < LengthOverSizeUpperLimit) + return kTRUE; + } + + + return kFALSE; + +} + + + + + @@ -340,5 +420,8 @@ // ls2: ls^2 // ct: Cos(ZA.) - Cos(fThetaOffset) - // dd2: (DIST- fDistOffset)^2 + // dd2: DIST^2 - fDistOffset^2 + + // lconc: log10CONC + Double_t limit; @@ -488,4 +571,14 @@ + + // The parametrization of upper cut in CONC is just provisional. + // The aim is to get rid of triangular events. + // The parametrization is + // Log10(Conc)Uppercut = m*(log(Size) - log(SizeOffset)) + b + + + Double_t log10conc = TMath::Log10(conc); + + // in the parameterization of the cuts // the dist parameter used is the one computed from the source position, @@ -508,5 +601,13 @@ const Double_t asym = asym0 * fMm2Deg; - + + // MuuonLikeEvent is a variable which is set + // to kFALSE/kTRUE by the filter cuts implemented + // in MuonFilter function. If kTRUE, the event + // will be rejected. + +// DM: + //Bool_t MuonLikeEvent = MuonFilter (size, length); + // computation of the cut limits @@ -524,6 +625,6 @@ Double_t asymlow = CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2); - Double_t concup = CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2); - Double_t conclow = CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2); + Double_t log10concup = CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2); + Double_t log10conclow = CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2); Double_t leakageup = CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2); @@ -533,10 +634,16 @@ Double_t lengthoverwidth = length/width; + + // Cut in CONC is parametrized + Double_t hadronness = 0.0; - // If the cut values computed before for Dist, Length and Width exceed the upper limits set - // by the user through function MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits, - // (or the default values defined in constructor); such cut values are replaced by the + // If the cut values computed before for Dist, + // Length and Width exceed the upper limits set + // by the user through function + // MSupercutsCalcONOFF::SetHillasDistLengthWidthUpperLowerLimits, + // (or the default values defined in constructor); such cut values + // are replaced by the // the limit values @@ -576,4 +683,9 @@ } + if (log10concup > fLog10ConcUpperLimit) + { + log10concup = fLog10ConcUpperLimit; + } + if (leakageup > fLeakage1UpperLimit) @@ -582,13 +694,9 @@ } - - + // Upper cut in leakage 1 also implemented - - - if (// newdist > distup || @@ -600,5 +708,5 @@ leakage > leakageup || leakage < leakagelow || - lengthoverwidth < fLengthOverWidthUpperLimit + lengthoverwidth < fLengthOverWidthUpperLimit || // //asym < asymup && @@ -608,7 +716,9 @@ //dist > distlow && - //conc < concup && - //conc > conclow && + log10conc > log10concup + // log10conc < log10conclow // + //log10conc > log10concup || + // DM MuonLikeEvent // if KTRUE, event is rejected ) @@ -641,5 +751,5 @@ ShowerParams[3] = newdist; ShowerParams[4] = asym; - ShowerParams[5] = conc; + ShowerParams[5] = log10conc; ShowerParams[6] = leakage; ShowerParams[7] = theta; @@ -666,6 +776,6 @@ SuperCutParams[6] = asymup; SuperCutParams[7] = asymlow; - SuperCutParams[8] = concup; - SuperCutParams[9] = conclow; + SuperCutParams[8] = log10concup; + SuperCutParams[9] = log10conclow; SuperCutParams[10] = leakageup; SuperCutParams[11] = leakagelow; Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.h (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MSupercutsCalcONOFF.h (revision 6232) @@ -94,4 +94,6 @@ Double_t fWidthUpperLimit; + Double_t fLog10ConcUpperLimit; + Double_t fDistLowerLimit; Double_t fLengthLowerLimit; @@ -149,4 +151,12 @@ Double_t widthlow); + // Function implementing a filter for muon induced images +// (Keichi cuts October 18th 2004) + + Bool_t MuonFilter (Double_t Size, + Double_t Length); + + + void InitMapping(MHMatrix *mat); // use quantity ThetaOrig.fVal at theta Index: trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MTSupercutsApplied.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MTSupercutsApplied.cc (revision 5329) +++ trunk/MagicSoft/Mars/mtemp/mmpi/SupercutsONOFFClasses/MTSupercutsApplied.cc (revision 6232) @@ -130,6 +130,6 @@ StringVector[6] = ("AsymUp/D"); StringVector[7] = ("AsymLow/D"); - StringVector[8] = ("ConcUp/D"); - StringVector[9] = ("ConcLow/D"); + StringVector[8] = ("Log10ConcUp/D"); + StringVector[9] = ("Log10ConcLow/D"); StringVector[10] = ("LeakageUp/D"); StringVector[11] = ("LeakageLow/D"); @@ -155,5 +155,5 @@ StringVector2[3] = ("NewDist/D"); StringVector2[4] = ("Asym/D"); - StringVector2[5] = ("Conc/D"); + StringVector2[5] = ("Log10Conc/D"); StringVector2[6] = ("Leakage/D"); StringVector2[7] = ("Theta/D");