Changeset 1844

Timestamp:

03/19/03 16:59:38 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• manalysis/MEnergyEstParam.cc (modified) (9 diffs)
• manalysis/MEnergyEstParam.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -37 +37 @@
     * mtools/Makefile, mtools/ToolsLinkDef.h:
       - added MChisqEval
+
+    * manalysis/MEnergyEstParam.[h,cc]:
+      - slight changes
 
 

trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc

@@ -52 +52 @@
 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 +92 @@
 //
 MEnergyEstParam::MEnergyEstParam(const char *hillas, const char *name, const char *title)
-   : fMatrix(NULL), fA(4), fB(7)
+   : fMatrix(NULL)
 {
     fName  = name  ? name  : "MEnergyEstParam";
@@ -88 +107 @@
     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 +178 @@
 //
 // 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 +196 @@
 // 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 +205 @@
 
     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 +249 @@
     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 +281 @@
 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 +315 @@
         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();
     }
 

trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h

@@ -41 +41 @@
 
 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 +47 @@
     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