/* ======================================================================== *\ ! ! * ! * 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): Eva Domingo , 12/2004 ! Wolfgang Wittek, 12/2004 ! ! Copyright: MAGIC Software Development, 2000-2005 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MDispCalc // // // // This task calculates Disp // // (parameters are taken from the container MDisp) // // Also the matrix of variables to be used in the Disp optimization // // is defined // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MDispCalc.h" #include #include #include "MGeomCam.h" #include "MSrcPosCam.h" #include "MHillas.h" #include "MHillasExt.h" #include "MNewImagePar.h" #include "MMcEvt.hxx" #include "MPointingPos.h" #include "MImageParDisp.h" #include "MDisp.h" #include "MHMatrix.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MDispCalc); using namespace std; // -------------------------------------------------------------------------- // // constructor // MDispCalc::MDispCalc(const char *imagepardispname, const char *dispname, const char *name, const char *title) : fImageParDispName(imagepardispname), fDispName(dispname) { fName = name ? name : "MDispCalc"; fTitle = title ? title : "Class to calculate Disp"; fMatrix = NULL; // initialize mapping array dimension to the number of columns of fMatrix fMap.Set(18); } // -------------------------------------------------------------------------- // Int_t MDispCalc::PreProcess(MParList *pList) { // look for the defined camera geometry to get mm to deg conversion factor MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); if (!cam) { *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; return kFALSE; } fMm2Deg = cam->GetConvMm2Deg(); // look for Disp related containers fImageParDisp = (MImageParDisp*)pList->FindCreateObj("MImageParDisp", fImageParDispName); if (!fImageParDisp) { *fLog << err << fImageParDispName << " [MImageParDisp] not found... aborting." << endl; return kFALSE; } fDisp = (MDisp*)pList->FindObject(fDispName, "MDisp"); if (!fDisp) { *fLog << err << fDispName << " [MDisp] not found... aborting." << endl; return kFALSE; } //----------------------------------------------------------- // if fMatrix exists, skip preprocessing and just read events from matrix; // if doesn't exist, read variables values from containers, so look for them if (fMatrix) return kTRUE; fSrcPos = (MSrcPosCam*)pList->FindObject("MSrcPosCam"); if (!fSrcPos) { *fLog << err << "MSrcPosCam not found... aborting." << endl; return kFALSE; } fHil = (MHillas*)pList->FindObject("MHillas"); if (!fHil) { *fLog << err << "MHillas not found... aborting." << endl; return kFALSE; } fHilExt = (MHillasExt*)pList->FindObject("MHillasExt"); if (!fHilExt) { *fLog << err << "MHillasExt not found... aborting." << endl; return kFALSE; } fNewPar = (MNewImagePar*)pList->FindObject("MNewImagePar"); if (!fNewPar) { *fLog << err << "MNewImagePar not found... aborting." << endl; return kFALSE; } fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); if (!fMcEvt) { *fLog << err << "MMcEvt not found... This is not a MC file," << " you are not trying to optimize Disp, just calculating it." << endl; // return kFALSE; } fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); if (!fPointing) { *fLog << err << "MPointingPos not found... " << "MMcEvt is going to be used to get Theta and Phi." << endl; // return kFALSE; } //------------------------------------------ return kTRUE; } // -------------------------------------------------------------------------- // // You can use this function if you want to use a MHMatrix instead of the // given containers for the Disp optimization. This function adds all // necessary columns to the given matrix. Afterwards, you should fill // the matrix with the corresponding data (eg from a file by using // MHMatrix::Fill). Then, if you loop through the matrix (eg using // MMatrixLoop), MDispCalc::Process will take the values from the matrix // instead of the containers. // void MDispCalc::InitMapping(MHMatrix *mat) { if (fMatrix) return; fMatrix = mat; fMap[0] = fMatrix->AddColumn("MSrcPosCam.fX"); fMap[1] = fMatrix->AddColumn("MSrcPosCam.fY"); fMap[2] = fMatrix->AddColumn("MHillas.fMeanX"); fMap[3] = fMatrix->AddColumn("MHillas.fMeanY"); fMap[4] = fMatrix->AddColumn("MHillas.fDelta"); fMap[5] = fMatrix->AddColumn("MHillas.fSize"); fMap[6] = fMatrix->AddColumn("MHillas.fWidth"); fMap[7] = fMatrix->AddColumn("MHillas.fLength"); fMap[8] = fMatrix->AddColumn("MNewImagePar.fConc"); fMap[9] = fMatrix->AddColumn("MNewImagePar.fLeakage1"); fMap[10] = fMatrix->AddColumn("MNewImagePar.fLeakage2"); fMap[11] = fMatrix->AddColumn("MHillasExt.fM3Long"); fMap[12] = fMatrix->AddColumn("MHillasExt.fAsym"); fMap[13] = fMatrix->AddColumn("MMcEvt.fEnergy"); fMap[14] = fMatrix->AddColumn("MMcEvt.fImpact"); fMap[15] = fMatrix->AddColumn("MMcEvt.fLongitmax"); if (fPointing) { fMap[16] = fMatrix->AddColumn("MPointingPos.fZd"); fMap[17] = fMatrix->AddColumn("MPointingPos.fAz"); } else { fMap[16] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg"); fMap[17] = fMatrix->AddColumn("MMcEvt.fTelescopePhi*kRad2Deg"); } } // -------------------------------------------------------------------------- // // Returns the Matrix and the mapped value for each Matrix column // MHMatrix* MDispCalc::GetMatrixMap(TArrayI &map) { map = fMap; return fMatrix; } // -------------------------------------------------------------------------- // // Returns a mapped value from the Matrix // Double_t MDispCalc::GetVal(Int_t i) const { Double_t val = (*fMatrix)[fMap[i]]; // *fLog << "MDispCalc::GetVal; i, fMatrix, fMap, val = " // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; return val; } // --------------------------------------------------------------------------- // // Calculate Disp // // Int_t MDispCalc::Process() { // get variables needed to compute disp from the matrix or the containers const Double_t size = fMatrix ? GetVal(5) : fHil->GetSize(); const Double_t width0 = fMatrix ? GetVal(6) : fHil->GetWidth(); const Double_t length0 = fMatrix ? GetVal(7) : fHil->GetLength(); const Double_t conc = fMatrix ? GetVal(8) : fNewPar->GetConc(); const Double_t leakage1 = fMatrix ? GetVal(9) : fNewPar->GetLeakage1(); const Double_t leakage2 = fMatrix ? GetVal(10) : fNewPar->GetLeakage2(); // convert to deg const Double_t width = width0 * fMm2Deg; const Double_t length = length0 * fMm2Deg; // create an array to pass the variables to MDisp::Calc Int_t numimagevars = 6; TArrayD imagevars(numimagevars); imagevars[0] = log10(size); imagevars[1] = width; imagevars[2] = length; imagevars[3] = conc; imagevars[4] = leakage1; imagevars[5] = leakage2; // compute Disp Double_t disp = fDisp->Calc(imagevars); // save value into MImageParDisp container fImageParDisp->SetDisp(disp); fImageParDisp->SetReadyToSave(); return kTRUE; } //===========================================================================