Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1843) +++ /trunk/MagicSoft/Mars/Changelog (revision 1844) @@ -37,4 +37,7 @@ * mtools/Makefile, mtools/ToolsLinkDef.h: - added MChisqEval + + * manalysis/MEnergyEstParam.[h,cc]: + - slight changes Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc (revision 1843) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc (revision 1844) @@ -52,16 +52,35 @@ void MEnergyEstParam::InitCoefficients() { - fA[0] = 39.667402; // [cm] - fA[1] = 143.081912; // [cm] - fA[2] = -395.501677; // [cm] - fA[3] = 0.006193; - - fB[0] = 45.037867; // [GeV] - fB[1] = 0.087222; // [GeV] - fB[2] = -0.208065; // [GeV/cm] - fB[3] = 78.164942; // [GeV] - fB[4] = -159.361283; // [GeV] - fB[5] = -0.130455; // [GeV/cm] - fB[6] = 0.051139; + fA.Set(5); + fA[0] = -2539; // [cm] + fA[1] = 900; // [cm] + fA[2] = 17.5; // [cm] + fA[3] = 0.006; + fA[4] = 58.3; + /* + fA[0] = 39.667402; // [cm] + fA[1] = 143.081912; // [cm] + fA[2] = -395.501677; // [cm] + fA[3] = 0.006193; + fA[4] = 0; + */ + + fB.Set(7); + fB[0] = -8.69; // [GeV] + fB[1] = -0.069; // [GeV] + fB[2] = 0.000655; // [GeV] + fB[3] = 0.0326; // [GeV] + fB[4] = 0.000225; // [GeV] + fB[5] = 4.13e-8; // [GeV] + fB[6] = 0.05; + /* + fB[0] = 45.037867; // [GeV] + fB[1] = 0.087222; // [GeV] + fB[2] = -0.208065; // [GeV/cm] + fB[3] = 78.164942; // [GeV] + fB[4] = -159.361283; // [GeV] + fB[5] = -0.130455; // [GeV/cm] + fB[6] = 0.051139; + */ } @@ -73,5 +92,5 @@ // MEnergyEstParam::MEnergyEstParam(const char *hillas, const char *name, const char *title) - : fMatrix(NULL), fA(4), fB(7) + : fMatrix(NULL) { fName = name ? name : "MEnergyEstParam"; @@ -88,8 +107,8 @@ InitCoefficients(); - AddToBranchList("MMcEvt.fTheta"); - AddToBranchList(fHillasName+"fSize"); - AddToBranchList(fHillasName+"fWidth"); - AddToBranchList(fHillasName+"fLength"); + AddToBranchList("MMcEvt.fTelescopeTheta"); + AddToBranchList(fHillasName+".fSize"); + AddToBranchList(fHillasName+".fWidth"); + AddToBranchList(fHillasName+".fLength"); } @@ -159,7 +178,7 @@ // // Set the four coefficients for the estimation of the impact parameter. -// Argument must ba a TArrayD of size 4. -// -void MEnergyEstParam::SetCoeffA(TArrayD arr) +// Argument must ba a TArrayD of size 5. +// +void MEnergyEstParam::SetCoeffA(const TArrayD &arr) { if (arr.GetSize() != fA.GetSize()) @@ -177,5 +196,5 @@ // Argument must ba a TArrayD of size 7. // -void MEnergyEstParam::SetCoeffB(TArrayD arr) +void MEnergyEstParam::SetCoeffB(const TArrayD &arr) { if (arr.GetSize() != fB.GetSize()) @@ -186,4 +205,21 @@ fB = arr; +} + +// -------------------------------------------------------------------------- +// +// Set the four coefficients for the estimation of impact and energy. +// Argument must ba a TArrayD of size 12. +// +void MEnergyEstParam::SetCoeff(const TArrayD &arr) +{ + if (arr.GetSize() != fA.GetSize()+fB.GetSize()) + { + *fLog << err << dbginf << "Error - Wrong number of coefficients!" << endl; + return; + } + + fA = TArrayD(fA.GetSize(), arr.GetArray()); + fB = TArrayD(fB.GetSize(), arr.GetArray() + fA.GetSize()); } @@ -213,5 +249,5 @@ fMatrix = mat; - fMap[0] = fMatrix->AddColumn("MMcEvt.fTheta"); + fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta"); fMap[1] = fMatrix->AddColumn(fHillasName+".fWidth"); fMap[2] = fMatrix->AddColumn(fHillasName+".fLength"); @@ -245,18 +281,24 @@ Bool_t MEnergyEstParam::Process() { - const Double_t theta = fMatrix ? GetVal(0) : fMc->GetTheta(); + const Double_t theta = fMatrix ? GetVal(0) : fMc->GetTelescopeTheta(); const Double_t width = fMatrix ? GetVal(1) : fHillas->GetWidth(); const Double_t length = fMatrix ? GetVal(2) : fHillas->GetLength(); const Double_t size = fMatrix ? GetVal(3) : fHillas->GetSize(); - const Double_t k = 1/cos(theta); + const Double_t k = 1./cos(theta); const Double_t k2 = k*k; - const Double_t i0 = k * (1+fA[3]*k)/(1+fA[3]); + const Double_t i0 = k * (fA[3]*k+1)/(fA[3]+1); const Double_t i1 = fA[0] + fA[2]*width; - const Double_t e0 = k2 * (1+fB[6]*k2)/(1+fB[6]); - const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*length; - const Double_t e2 = fB[2] + fB[5]*length; + const Double_t e0 = k2 * (fB[6]*k2+1)/(fB[6]+1); + + /* MY PARAM */ + const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*width; + const Double_t e2 = fB[2] + fB[5]*size*width; + + /* MARCOS */ + //const Double_t e1 = fB[0] + fB[1]*size + fB[3]*length + fB[4]*size*length; + //const Double_t e2 = fB[2] + fB[5]*length; TIter NextH(fHillasSrc); @@ -273,9 +315,31 @@ const Double_t dist = fMatrix ? GetVal(col++) : ((MHillasSrc*)NextH())->GetDist(); - const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm] - const Double_t er = -e0 * (e1 + e2*ir); // [GeV] + + /* MARCOS */ + //const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm] + /*const*/// Double_t er = e0 * (e1 + e2*ir); // [GeV] + + /* MY PARAM */ + // if (width==0) return kCONTINUE; + const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm] + Double_t er = e0 * (e1 + e2*ir); // [GeV] + + /* MKA */ + //Double_t er = e0 * (fB[0] + fB[1]*size/width + fB[2]*ir /*+ d*leakage*/); + + if (width==0) + return kCONTINUE; + + if (TMath::IsNaN(ir)) + *fLog << all << theta << " " << width << " " << length << " " << size << " " << dist << endl; + if (TMath::IsNaN(er)) + { + *fLog << all << theta << " " << width << " " << length << " " << size << " " << dist << endl; + er = 0; + } est->SetEnergy(er); est->SetImpact(ir); + est->SetReadyToSave(); } Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 1843) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 1844) @@ -41,5 +41,5 @@ public: - MEnergyEstParam(const char *hil="MHillas", const char *name=NULL, const char *title=NULL); + MEnergyEstParam(const char *hil="Hillas", const char *name=NULL, const char *title=NULL); ~MEnergyEstParam(); @@ -47,10 +47,13 @@ Bool_t Process(); - void Add(const TString hillas, const TString energy); + void Add(const TString hillas, const TString energy="MEnergyEst"); void InitMapping(MHMatrix *mat); - void SetCoeffA(TArrayD arr); - void SetCoeffB(TArrayD arr); + Int_t GetNumCoeff() const { return fA.GetSize()+fB.GetSize(); } + + void SetCoeff(const TArrayD &arr); + void SetCoeffA(const TArrayD &arr); + void SetCoeffB(const TArrayD &arr); ClassDef(MEnergyEstParam, 0) // Task to estimate the energy