| 1 | /* ======================================================================== *\
|
|---|
| 2 | !
|
|---|
| 3 | ! *
|
|---|
| 4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
|---|
| 5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
|---|
| 6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
|---|
| 7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
|---|
| 8 | ! *
|
|---|
| 9 | ! * Permission to use, copy, modify and distribute this software and its
|
|---|
| 10 | ! * documentation for any purpose is hereby granted without fee,
|
|---|
| 11 | ! * provided that the above copyright notice appear in all copies and
|
|---|
| 12 | ! * that both that copyright notice and this permission notice appear
|
|---|
| 13 | ! * in supporting documentation. It is provided "as is" without express
|
|---|
| 14 | ! * or implied warranty.
|
|---|
| 15 | ! *
|
|---|
| 16 | !
|
|---|
| 17 | !
|
|---|
| 18 | ! Author(s): Eva Domingo , 12/2004 <mailto:domingo@ifae.es>
|
|---|
| 19 | ! Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2005
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 27 | // //
|
|---|
| 28 | // MDispCalc //
|
|---|
| 29 | // //
|
|---|
| 30 | // This task calculates Disp with a given parameterization //
|
|---|
| 31 | // (parameters are stored in the MDispParameters container) //
|
|---|
| 32 | // Also the matrix of variables to be used in the Disp //
|
|---|
| 33 | // parameterization is defined //
|
|---|
| 34 | // //
|
|---|
| 35 | // //
|
|---|
| 36 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 37 | #include "MDispCalc.h"
|
|---|
| 38 |
|
|---|
| 39 | #include <math.h>
|
|---|
| 40 | #include <TArray.h>
|
|---|
| 41 |
|
|---|
| 42 | #include "MGeomCam.h"
|
|---|
| 43 | #include "MSrcPosCam.h"
|
|---|
| 44 | #include "MHillas.h"
|
|---|
| 45 | #include "MHillasExt.h"
|
|---|
| 46 | #include "MNewImagePar.h"
|
|---|
| 47 | #include "MMcEvt.hxx"
|
|---|
| 48 | #include "MPointingPos.h"
|
|---|
| 49 | #include "MImageParDisp.h"
|
|---|
| 50 | #include "MDispParameters.h"
|
|---|
| 51 |
|
|---|
| 52 | #include "MHMatrix.h"
|
|---|
| 53 |
|
|---|
| 54 | #include "MParList.h"
|
|---|
| 55 |
|
|---|
| 56 | #include "MLog.h"
|
|---|
| 57 | #include "MLogManip.h"
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 | ClassImp(MDispCalc);
|
|---|
| 61 |
|
|---|
| 62 | using namespace std;
|
|---|
| 63 |
|
|---|
| 64 | static const Int_t nPars=5; // number of parameters used in the Disp expression
|
|---|
| 65 |
|
|---|
| 66 | enum dispVar_t {kSize,kWidth,kLength,kConc,kLeakage1,kLeakage2,kTotVar}; // enum variables for the
|
|---|
| 67 | // matrix columns mapping
|
|---|
| 68 |
|
|---|
| 69 | // --------------------------------------------------------------------------
|
|---|
| 70 | //
|
|---|
| 71 | // constructor
|
|---|
| 72 | //
|
|---|
| 73 | MDispCalc::MDispCalc(const char *imagepardispname,
|
|---|
| 74 | const char *dispparametersname,
|
|---|
| 75 | const char *name, const char *title)
|
|---|
| 76 | : fImageParDispName(imagepardispname),
|
|---|
| 77 | fDispParametersName(dispparametersname),
|
|---|
| 78 | fLogSize0(3.), fLength0(0.1), fWidth0(0.05),
|
|---|
| 79 | fConc0(0.35), fLeakage10(0.05), fLeakage20(0.1)
|
|---|
| 80 | {
|
|---|
| 81 | fName = name ? name : "MDispCalc";
|
|---|
| 82 | fTitle = title ? title : "Class to calculate Disp";
|
|---|
| 83 |
|
|---|
| 84 | // object of MDispParamteres that will hold the Disp parameters
|
|---|
| 85 | fDispParameters = new MDispParameters(fDispParametersName);
|
|---|
| 86 | // initialize the size of the Disp parameters and parameters stepsizes arrays
|
|---|
| 87 | fDispParameters->GetParameters().Set(nPars);
|
|---|
| 88 | fDispParameters->GetStepsizes().Set(nPars);
|
|---|
| 89 | // set Disp parameters to their default values
|
|---|
| 90 | fDispParameters->InitParameters();
|
|---|
| 91 |
|
|---|
| 92 | fMatrix = NULL;
|
|---|
| 93 |
|
|---|
| 94 | // initialize mapping array dimension to the number of columns of fMatrix
|
|---|
| 95 | fMap.Set(kTotVar);
|
|---|
| 96 | }
|
|---|
| 97 |
|
|---|
| 98 |
|
|---|
| 99 | // --------------------------------------------------------------------------
|
|---|
| 100 | //
|
|---|
| 101 | // Default destructor.
|
|---|
| 102 | //
|
|---|
| 103 | MDispCalc::~MDispCalc()
|
|---|
| 104 | {
|
|---|
| 105 | delete fDispParameters;
|
|---|
| 106 | }
|
|---|
| 107 |
|
|---|
| 108 |
|
|---|
| 109 | // --------------------------------------------------------------------------
|
|---|
| 110 | //
|
|---|
| 111 | Int_t MDispCalc::PreProcess(MParList *pList)
|
|---|
| 112 | {
|
|---|
| 113 | // look for the defined camera geometry to get mm to deg conversion factor
|
|---|
| 114 | MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
|
|---|
| 115 | if (!cam)
|
|---|
| 116 | {
|
|---|
| 117 | *fLog << err << "MGeomCam (Camera Geometry) not found... aborting."
|
|---|
| 118 | << endl;
|
|---|
| 119 | return kFALSE;
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | fMm2Deg = cam->GetConvMm2Deg();
|
|---|
| 123 |
|
|---|
| 124 |
|
|---|
| 125 | // look for Disp related containers
|
|---|
| 126 | fImageParDisp = (MImageParDisp*)pList->FindCreateObj("MImageParDisp",
|
|---|
| 127 | fImageParDispName);
|
|---|
| 128 | if (!fImageParDisp)
|
|---|
| 129 | {
|
|---|
| 130 | *fLog << err << fImageParDispName
|
|---|
| 131 | << " [MImageParDisp] not found... aborting." << endl;
|
|---|
| 132 | return kFALSE;
|
|---|
| 133 | }
|
|---|
| 134 |
|
|---|
| 135 |
|
|---|
| 136 | //-----------------------------------------------------------
|
|---|
| 137 | // if fMatrix exists, skip preprocessing and just read events from matrix;
|
|---|
| 138 | // if doesn't exist, read variables values from containers, so look for them
|
|---|
| 139 | if (fMatrix)
|
|---|
| 140 | return kTRUE;
|
|---|
| 141 |
|
|---|
| 142 | fSrcPos = (MSrcPosCam*)pList->FindObject("MSrcPosCam");
|
|---|
| 143 | if (!fSrcPos)
|
|---|
| 144 | {
|
|---|
| 145 | *fLog << err << "MSrcPosCam not found... aborting." << endl;
|
|---|
| 146 | return kFALSE;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | fHil = (MHillas*)pList->FindObject("MHillas");
|
|---|
| 150 | if (!fHil)
|
|---|
| 151 | {
|
|---|
| 152 | *fLog << err << "MHillas not found... aborting." << endl;
|
|---|
| 153 | return kFALSE;
|
|---|
| 154 | }
|
|---|
| 155 |
|
|---|
| 156 | fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
|
|---|
| 157 | if (!fHilExt)
|
|---|
| 158 | {
|
|---|
| 159 | *fLog << err << "MHillasExt not found... aborting." << endl;
|
|---|
| 160 | return kFALSE;
|
|---|
| 161 | }
|
|---|
| 162 |
|
|---|
| 163 | fNewPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
|
|---|
| 164 | if (!fNewPar)
|
|---|
| 165 | {
|
|---|
| 166 | *fLog << err << "MNewImagePar not found... aborting." << endl;
|
|---|
| 167 | return kFALSE;
|
|---|
| 168 | }
|
|---|
| 169 |
|
|---|
| 170 | fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
|
|---|
| 171 | if (!fMcEvt)
|
|---|
| 172 | {
|
|---|
| 173 | *fLog << err << "MMcEvt not found... This is not a MC file,"
|
|---|
| 174 | << " you are not trying to optimize Disp, just calculating it."
|
|---|
| 175 | << endl;
|
|---|
| 176 | // return kFALSE;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
|
|---|
| 180 | if (!fPointing)
|
|---|
| 181 | {
|
|---|
| 182 | *fLog << err << "MPointingPos not found... aborting." << endl;
|
|---|
| 183 | return kFALSE;
|
|---|
| 184 | }
|
|---|
| 185 |
|
|---|
| 186 | //------------------------------------------
|
|---|
| 187 |
|
|---|
| 188 | return kTRUE;
|
|---|
| 189 | }
|
|---|
| 190 |
|
|---|
| 191 |
|
|---|
| 192 | // --------------------------------------------------------------------------
|
|---|
| 193 | //
|
|---|
| 194 | // You can use this function if you want to use a MHMatrix instead of the
|
|---|
| 195 | // given containers for the Disp parameterization. This function adds all
|
|---|
| 196 | // necessary columns to the given matrix. Afterwards, you should fill
|
|---|
| 197 | // the matrix with the corresponding data (eg from a file by using
|
|---|
| 198 | // MHMatrix::Fill). Then, if you loop through the matrix (eg using
|
|---|
| 199 | // MMatrixLoop), MDispCalc::Calc will take the values from the matrix
|
|---|
| 200 | // instead of the containers.
|
|---|
| 201 | //
|
|---|
| 202 | void MDispCalc::InitMapping(MHMatrix *mat)
|
|---|
| 203 | {
|
|---|
| 204 | if (fMatrix)
|
|---|
| 205 | return;
|
|---|
| 206 |
|
|---|
| 207 | fMatrix = mat;
|
|---|
| 208 |
|
|---|
| 209 | fMap[kSize] = fMatrix->AddColumn("MHillas.fSize");
|
|---|
| 210 | fMap[kWidth] = fMatrix->AddColumn("MHillas.fWidth");
|
|---|
| 211 | fMap[kLength] = fMatrix->AddColumn("MHillas.fLength");
|
|---|
| 212 |
|
|---|
| 213 | fMap[kConc] = fMatrix->AddColumn("MNewImagePar.fConc");
|
|---|
| 214 | fMap[kLeakage1] = fMatrix->AddColumn("MNewImagePar.fLeakage1");
|
|---|
| 215 | fMap[kLeakage2] = fMatrix->AddColumn("MNewImagePar.fLeakage2");
|
|---|
| 216 | }
|
|---|
| 217 |
|
|---|
| 218 |
|
|---|
| 219 | // --------------------------------------------------------------------------
|
|---|
| 220 | //
|
|---|
| 221 | // Returns a mapped value from the Matrix
|
|---|
| 222 | //
|
|---|
| 223 | Double_t MDispCalc::GetVal(Int_t i) const
|
|---|
| 224 | {
|
|---|
| 225 | Double_t val = (*fMatrix)[fMap[i]];
|
|---|
| 226 |
|
|---|
| 227 | // *fLog << "MDispCalc::GetVal; i, fMatrix, fMap, val = "
|
|---|
| 228 | // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl;
|
|---|
| 229 |
|
|---|
| 230 | return val;
|
|---|
| 231 | }
|
|---|
| 232 |
|
|---|
| 233 |
|
|---|
| 234 | // ---------------------------------------------------------------------------
|
|---|
| 235 | //
|
|---|
| 236 | // Calculate Disp
|
|---|
| 237 | //
|
|---|
| 238 | //
|
|---|
| 239 | Int_t MDispCalc::Process()
|
|---|
| 240 | {
|
|---|
| 241 | // get variables needed to compute disp from the matrix or the containers
|
|---|
| 242 | const Double_t size = fMatrix ? GetVal(kSize) : fHil->GetSize();
|
|---|
| 243 | const Double_t width0 = fMatrix ? GetVal(kWidth) : fHil->GetWidth();
|
|---|
| 244 | const Double_t length0 = fMatrix ? GetVal(kLength) : fHil->GetLength();
|
|---|
| 245 | const Double_t conc = fMatrix ? GetVal(kConc) : fNewPar->GetConc();
|
|---|
| 246 | const Double_t leakage1 = fMatrix ? GetVal(kLeakage1) : fNewPar->GetLeakage1();
|
|---|
| 247 | const Double_t leakage2 = fMatrix ? GetVal(kLeakage2) : fNewPar->GetLeakage2();
|
|---|
| 248 |
|
|---|
| 249 | // convert to deg
|
|---|
| 250 | const Double_t width = width0 * fMm2Deg;
|
|---|
| 251 | const Double_t length = length0 * fMm2Deg;
|
|---|
| 252 |
|
|---|
| 253 | // create an array to pass the variables to MDispCalc::Calc
|
|---|
| 254 | TArrayD imagevars(kTotVar);
|
|---|
| 255 | imagevars[kSize] = log10(size);
|
|---|
| 256 | imagevars[kWidth] = width;
|
|---|
| 257 | imagevars[kLength] = length;
|
|---|
| 258 | imagevars[kConc] = conc;
|
|---|
| 259 | imagevars[kLeakage1] = leakage1;
|
|---|
| 260 | imagevars[kLeakage2] = leakage2;
|
|---|
| 261 |
|
|---|
| 262 | // compute Disp
|
|---|
| 263 | Double_t disp = Calc(imagevars);
|
|---|
| 264 |
|
|---|
| 265 | // save value into MImageParDisp container
|
|---|
| 266 | fImageParDisp->SetDisp(disp);
|
|---|
| 267 | fImageParDisp->SetReadyToSave();
|
|---|
| 268 |
|
|---|
| 269 | return kTRUE;
|
|---|
| 270 | }
|
|---|
| 271 |
|
|---|
| 272 |
|
|---|
| 273 | // --------------------------------------------------------------------------
|
|---|
| 274 | //
|
|---|
| 275 | // Set the Disp parameterization and Calculate Disp
|
|---|
| 276 | //
|
|---|
| 277 | //
|
|---|
| 278 | Double_t MDispCalc::Calc(TArrayD &imagevar)
|
|---|
| 279 | {
|
|---|
| 280 | // get parameters
|
|---|
| 281 | const TArrayD& p = fDispParameters->GetParameters();
|
|---|
| 282 |
|
|---|
| 283 | // get image variables to be used in the Disp parameterization
|
|---|
| 284 | Double_t logsize = imagevar[0];
|
|---|
| 285 | Double_t width = imagevar[1];
|
|---|
| 286 | Double_t length = imagevar[2];
|
|---|
| 287 | // Double_t conc = imagevar[3];
|
|---|
| 288 | // Double_t leakage1 = imagevar[4];
|
|---|
| 289 | // Double_t leakage2 = imagevar[5];
|
|---|
| 290 |
|
|---|
| 291 | // Disp parameterization to be optimized
|
|---|
| 292 | // Note: fLogSize0, fLength0... variables are introduced to uncorrelate as much
|
|---|
| 293 | // as possible the parameters in the Disp expression, with the purpose of
|
|---|
| 294 | // helping the minimization algorithm to converge. They are set approx.
|
|---|
| 295 | // to their distribution mean value in the MDisp constructor, but can be
|
|---|
| 296 | // changed using the corresponding set function.
|
|---|
| 297 | //
|
|---|
| 298 |
|
|---|
| 299 | // Double_t disp = p[0] + p[1]*(logsize-fLogSize0)
|
|---|
| 300 | // + (p[2] + p[3]*(logsize-fLogSize0))*width/length;
|
|---|
| 301 |
|
|---|
| 302 | /*
|
|---|
| 303 | Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0)
|
|---|
| 304 | + (p[2] + p[3]*(logsize-fLogSize0))*width/length;
|
|---|
| 305 | */
|
|---|
| 306 |
|
|---|
| 307 | // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0)
|
|---|
| 308 | // + (p[2] + p[3]*(logsize-fLogSize0))*length/width;
|
|---|
| 309 |
|
|---|
| 310 | // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0)
|
|---|
| 311 | // + (p[2] + p[3]*(logsize-fLogSize0) + p[5]*(length-fLength0))*width/length;
|
|---|
| 312 |
|
|---|
| 313 | // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(width-fWidth0)
|
|---|
| 314 | // + (p[2] + p[3]*(logsize-fLogSize0))*width/length; // + p[5]*(width-fWidth0))*width/length;
|
|---|
| 315 |
|
|---|
| 316 | // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(conc-fConc0)
|
|---|
| 317 | // + (p[2] + p[3]*(logsize-fLogSize0))*width/length; // + p[5]*(conc-fConc0))*width/length;
|
|---|
| 318 |
|
|---|
| 319 | // Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0)
|
|---|
| 320 | // + p[2]*pow(logsize-fLogSize0,2)
|
|---|
| 321 | // + p[3]*pow(logsize-fLogSize0,3)
|
|---|
| 322 | // + p[4]*pow(logsize-fLogSize0,4) )
|
|---|
| 323 | // *(1-width/length);
|
|---|
| 324 |
|
|---|
| 325 | Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0)
|
|---|
| 326 | + p[2]*pow(logsize-fLogSize0,2)
|
|---|
| 327 | + p[3]*pow(logsize-fLogSize0,3)
|
|---|
| 328 | + p[4]*pow(logsize-fLogSize0,4) )
|
|---|
| 329 | *( 1/(1+width/length) );
|
|---|
| 330 |
|
|---|
| 331 | /*
|
|---|
| 332 | Double_t disp = ( p[0] + p[1]*(logsize)
|
|---|
| 333 | + p[2]*pow(logsize,2)
|
|---|
| 334 | + p[3]*pow(logsize,3)
|
|---|
| 335 | + p[4]*pow(logsize,4) )
|
|---|
| 336 | *( 1/(1+width/length) );
|
|---|
| 337 | */
|
|---|
| 338 |
|
|---|
| 339 | return disp;
|
|---|
| 340 | }
|
|---|
| 341 |
|
|---|
| 342 | //===========================================================================
|
|---|
| 343 |
|
|---|
| 344 |
|
|---|
| 345 |
|
|---|
| 346 |
|
|---|
