/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Wolfgang Wittek 04/2003 ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MCT1SupercutsCalc // // // ///////////////////////////////////////////////////////////////////////////// #include "MCT1SupercutsCalc.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MHillasExt.h" #include "MHillasSrc.h" #include "MMcEvt.hxx" #include "MCerPhotEvt.h" #include "MGeomCam.h" #include "MHadronness.h" ClassImp(MCT1SupercutsCalc); void MCT1SupercutsCalc::InitParams() { //--------------------------------- // cut parameters fLengthUp[0] = 0.315585; fLengthUp[1] = 0.001455; fLengthUp[2] = 0.203198; fLengthUp[3] = 0.005532; fLengthUp[4] =-0.001670; fLengthUp[5] =-0.020362; fLengthUp[6] = 0.007388; fLengthUp[7] =-0.013463; fWidthUp[0] = 0.145412; fWidthUp[1] =-0.001771; fWidthUp[2] = 0.054462; fWidthUp[3] = 0.022280; fWidthUp[4] =-0.009893; fWidthUp[5] = 0.056353; fWidthUp[6] = 0.020711; fWidthUp[7] =-0.016703; fDistUp[0] = 1.787943; fDistUp[1] = 0.; fDistUp[2] = 2.942310; fDistUp[3] = 0.199815; fDistUp[4] = 0.; fDistUp[5] = 0.249909; fDistUp[6] = 0.189697; fDistUp[7] = 0.; fLengthLo[0] = 0.151530; fLengthLo[1] = 0.028323; fLengthLo[2] = 0.510707; fLengthLo[3] = 0.053089; fLengthLo[4] = 0.013708; fLengthLo[5] = 2.357993; fLengthLo[6] = 0.000080; fLengthLo[7] =-0.007157; fWidthLo[0] = 0.089187; fWidthLo[1] =-0.006430; fWidthLo[2] = 0.074442; fWidthLo[3] = 0.003738; fWidthLo[4] =-0.004256; fWidthLo[5] =-0.014101; fWidthLo[6] = 0.006126; fWidthLo[7] =-0.002849; fDistLo[0] = 0.589406; fDistLo[1] = 0.; fDistLo[2] =-0.083964; fDistLo[3] =-0.007975; fDistLo[4] = 0.; fDistLo[5] = 0.045374; fDistLo[6] =-0.001750; fDistLo[7] = 0.; fAsymUp[0] = 0.061267; fAsymUp[1] = 0.014462; fAsymUp[2] = 0.014327; fAsymUp[3] = 0.014540; fAsymUp[4] = 0.013391; fAsymUp[5] = 0.012319; fAsymUp[6] = 0.010444; fAsymUp[7] = 0.008328; fAsymLo[0] =-0.012055; fAsymLo[1] = 0.009157; fAsymLo[2] = 0.005441; fAsymLo[3] = 0.000399; fAsymLo[4] = 0.001433; fAsymLo[5] =-0.002050; fAsymLo[6] =-0.000104; fAsymLo[7] =-0.001188; fAlphaUp[0] = 13.123440; fAlphaUp[1] = 0.; fAlphaUp[2] = 0.; fAlphaUp[3] = 0.; fAlphaUp[4] = 0.; fAlphaUp[5] = 0.; fAlphaUp[6] = 0.; fAlphaUp[7] = 0.; //--------------------------------- } // -------------------------------------------------------------------------- // MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, const char *hilsrcname, const char *name, const char *title) : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname) { fName = name ? name : "MCT1SupercutsCalc"; fTitle = title ? title : "Class to evaluate the Supercuts"; InitParams(); } // -------------------------------------------------------------------------- // Bool_t MCT1SupercutsCalc::PreProcess(MParList *pList) { fHil = (MHillas*)pList->FindObject(fHilName, "MHillas"); if (!fHil) { *fLog << err << fHilName << " [MHillas] not found... aborting." << endl; return kFALSE; } fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); if (!fHilSrc) { *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; return kFALSE; } fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); if (!fMcEvt) { *fLog << err << "MMcEvt not found... aborting." << endl; return kFALSE; } MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); if (!cam) { *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; return kFALSE; } fMm2Deg = cam->GetConvMm2Deg(); fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName); if (!fHadronness) return kFALSE; return kTRUE; } // -------------------------------------------------------------------------- // // Calculation of upper and lower limits // Double_t MCT1SupercutsCalc::CtsMCut(Double_t *a, Double_t ls, Double_t ct, Double_t ls2, Double_t dd2) { // define cut-function // // dNOMLOGSIZE = 4.1 (=log(60.0) // dNOMCOSZA = 1.0 // // a: array of cut parameters // ls: log(SIZE) - dNOMLOGSIZE // ls2: ls^2 // ct: Cos(ZA.) - dNOMCOSZA // dd2: DIST^2 const Double_t limit = a[0] + a[1] * dd2 + a[2] * ct + ls * (a[3] + a[4] * dd2 + a[5] * ct) + ls2 * (a[6] + a[7] * dd2); //*fLog << "MCT1SupercutsCalc::CtsMCut; *a = " // << *a << ", " << *(a+1) << ", " << *(a+2) << ", " // << *(a+3) << ", " << *(a+4) << ", " << *(a+5) << ", " // << *(a+6) << ", " << *(a+7) << endl; //*fLog << "MCT1SupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls // << ", " << ls2 << ", " << ct << ", " << dd2 << ", " // << limit << endl; return limit; } // --------------------------------------------------------------------------- // // Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich) // optimized for mkn 421 2001 data // // set hadronness to 0.25 if cuts are fullfilled // 0.75 otherwise // Bool_t MCT1SupercutsCalc::Process() { const Double_t dNOMLOGSIZE = 4.1; const Double_t dNOMCOSZA = 1.0; const Double_t newdist = fHilSrc->GetDist() * fMm2Deg; const Double_t xbar = fHil->GetMeanX(); const Double_t ybar = fHil->GetMeanY(); const Double_t dist = sqrt(xbar*xbar + ybar*ybar) * fMm2Deg; const Double_t dd2 = dist * dist; const Double_t dmls = log(fHil->GetSize()) - dNOMLOGSIZE; const Double_t dmls2 = dmls * dmls; const Double_t dmcza = cos(fMcEvt->GetTelescopeTheta()) - dNOMCOSZA; const Double_t length = fHil->GetLength() * fMm2Deg; const Double_t width = fHil->GetWidth() * fMm2Deg; //*fLog << "MCT1SupercutsCalc::Process; dmls, dmcza, dmls2, dd2 = " // << dmls << ", " << dmcza << ", " << dmls2 << ", " // << dd2 << endl; //*fLog << "MCT1SupercutsCalc::Process; newdist, dist, length, width = " // << newdist << ", " << dist << ", " << length << ", " // << width << endl; if (newdist < 1.05 && newdist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) && newdist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) && dist < 1.05 && length < CtsMCut (fLengthUp, dmls, dmcza, dmls2, dd2) && length > CtsMCut (fLengthLo, dmls, dmcza, dmls2, dd2) && width < CtsMCut (fWidthUp, dmls, dmcza, dmls2, dd2) && width > CtsMCut (fWidthLo, dmls, dmcza, dmls2, dd2) && //asym < CtsMCut (asymup, dmls, dmcza, dmls2, dd2) && //asym > CtsMCut (asymlow, dmls, dmcza, dmls2, dd2) && dist < CtsMCut (fDistUp, dmls, dmcza, dmls2, dd2) && dist > CtsMCut (fDistLo, dmls, dmcza, dmls2, dd2) ) fHadronness->SetHadronness(0.25); else fHadronness->SetHadronness(0.75); return kTRUE; }