| 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 | // MHDisp //
|
|---|
| 29 | // //
|
|---|
| 30 | // container holding the histograms for Disp //
|
|---|
| 31 | // also computes the minimization parameter of the Disp optimization //
|
|---|
| 32 | // //
|
|---|
| 33 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 34 | #include "MHDisp.h"
|
|---|
| 35 |
|
|---|
| 36 | #include <math.h>
|
|---|
| 37 |
|
|---|
| 38 | #include <TH1.h>
|
|---|
| 39 | #include <TH2.h>
|
|---|
| 40 | #include <TProfile.h>
|
|---|
| 41 | #include <TArrayI.h>
|
|---|
| 42 | #include <TPad.h>
|
|---|
| 43 | #include <TCanvas.h>
|
|---|
| 44 | #include <TStyle.h>
|
|---|
| 45 |
|
|---|
| 46 | #include "MGeomCam.h"
|
|---|
| 47 | #include "MSrcPosCam.h"
|
|---|
| 48 | #include "MHillas.h"
|
|---|
| 49 | #include "MHillasExt.h"
|
|---|
| 50 | #include "MNewImagePar.h"
|
|---|
| 51 | #include "MMcEvt.hxx"
|
|---|
| 52 | #include "MPointingPos.h"
|
|---|
| 53 | #include "MImageParDisp.h"
|
|---|
| 54 |
|
|---|
| 55 | #include "MHMatrix.h"
|
|---|
| 56 | #include "MParameters.h"
|
|---|
| 57 |
|
|---|
| 58 | #include "MParList.h"
|
|---|
| 59 |
|
|---|
| 60 | #include "MLog.h"
|
|---|
| 61 | #include "MLogManip.h"
|
|---|
| 62 |
|
|---|
| 63 | ClassImp(MHDisp);
|
|---|
| 64 |
|
|---|
| 65 | using namespace std;
|
|---|
| 66 |
|
|---|
| 67 | enum dispVar_t {kX,kY,kMeanX,kMeanY,kDelta,kSize,kM3Long,kAsym,
|
|---|
| 68 | kEnergy,kImpact,kLongitmax,kZd,kAz,kTotVar}; // enum variables for the
|
|---|
| 69 | // matrix columns mapping
|
|---|
| 70 |
|
|---|
| 71 | // --------------------------------------------------------------------------
|
|---|
| 72 | //
|
|---|
| 73 | // Constructor
|
|---|
| 74 | //
|
|---|
| 75 | MHDisp::MHDisp(const char *imagepardispname,
|
|---|
| 76 | const char *name, const char *title)
|
|---|
| 77 | : fImageParDispName(imagepardispname)
|
|---|
| 78 | {
|
|---|
| 79 | fName = name ? name : "MHDisp";
|
|---|
| 80 | fTitle = title ? title : "Histograms for Disp";
|
|---|
| 81 |
|
|---|
| 82 | fSelectedPos = 1; // default MHDisp flag (selects Correct Disp source position solution)
|
|---|
| 83 |
|
|---|
| 84 | fMatrix = NULL;
|
|---|
| 85 |
|
|---|
| 86 | // initialize mapping array dimension to the number of columns of fMatrix
|
|---|
| 87 | fMap.Set(kTotVar);
|
|---|
| 88 |
|
|---|
| 89 | //--------------------------------------------------
|
|---|
| 90 | // creating the Disp related histograms
|
|---|
| 91 |
|
|---|
| 92 | fHistEnergy = new TH1F("fHistEnergy",
|
|---|
| 93 | "log10(Energy)", 50, 1., 3.);
|
|---|
| 94 | fHistEnergy->SetDirectory(NULL);
|
|---|
| 95 | fHistEnergy->UseCurrentStyle();
|
|---|
| 96 | fHistEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 97 | fHistEnergy->SetYTitle("# events");
|
|---|
| 98 |
|
|---|
| 99 | fHistSize = new TH1F("fHistSize",
|
|---|
| 100 | "log10(Size)", 50, 2., 4.);
|
|---|
| 101 | fHistSize->SetDirectory(NULL);
|
|---|
| 102 | fHistSize->UseCurrentStyle();
|
|---|
| 103 | fHistSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 104 | fHistSize->SetYTitle("# events");
|
|---|
| 105 |
|
|---|
| 106 | fHistcosZA = new TH1F("fHistcosZA",
|
|---|
| 107 | "cos(Zenith Angle)", 10, 0., 1.);
|
|---|
| 108 | fHistcosZA->SetDirectory(NULL);
|
|---|
| 109 | fHistcosZA->UseCurrentStyle();
|
|---|
| 110 | fHistcosZA->SetXTitle("cos(Theta)");
|
|---|
| 111 | fHistcosZA->SetYTitle("# events");
|
|---|
| 112 |
|
|---|
| 113 | fSkymapXY = new TH2F("fSkymapXY",
|
|---|
| 114 | "Disp estimated source positions Skymap", 71, -2., 2., 71, -2., 2.);
|
|---|
| 115 | fSkymapXY->SetDirectory(NULL);
|
|---|
| 116 | fSkymapXY->UseCurrentStyle();
|
|---|
| 117 | fSkymapXY->SetXTitle("X Disp source position [deg]");
|
|---|
| 118 | fSkymapXY->SetYTitle("Y Disp source position [deg]");
|
|---|
| 119 |
|
|---|
| 120 | fHistMinPar = new TH1F("fHistMinPar" ,
|
|---|
| 121 | "Distance^2 between Disp and real srcpos", 100, 0., 2.);
|
|---|
| 122 | fHistMinPar->SetDirectory(NULL);
|
|---|
| 123 | fHistMinPar->UseCurrentStyle();
|
|---|
| 124 | fHistMinPar->SetXTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
|
|---|
| 125 | fHistMinPar->SetYTitle("# events");
|
|---|
| 126 |
|
|---|
| 127 | fHistDuDv = new TH2F("fHistDuDv",
|
|---|
| 128 | "Du vs. Dv (distances between Disp and real srcpos)",
|
|---|
| 129 | 100, -2., 2., 100, -2., 2.);
|
|---|
| 130 | fHistDuDv->SetDirectory(NULL);
|
|---|
| 131 | fHistDuDv->UseCurrentStyle();
|
|---|
| 132 | fHistDuDv->SetXTitle("Dv = transveral distance [deg]");
|
|---|
| 133 | fHistDuDv->SetYTitle("Du = longitudinal distance [deg]");
|
|---|
| 134 |
|
|---|
| 135 | fHistMinParEnergy = new TH2F("fHistMinParEnergy",
|
|---|
| 136 | "Minimization parameter vs. Energy", 50, 1., 3., 100, 0., 2.);
|
|---|
| 137 | fHistMinParEnergy->SetDirectory(NULL);
|
|---|
| 138 | fHistMinParEnergy->UseCurrentStyle();
|
|---|
| 139 | fHistMinParEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 140 | fHistMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
|
|---|
| 141 |
|
|---|
| 142 | fHistMinParSize = new TH2F("fHistMinParSize",
|
|---|
| 143 | "Minimization parameter vs. Size", 50, 2., 4., 100, 0., 2.);
|
|---|
| 144 | fHistMinParSize->SetDirectory(NULL);
|
|---|
| 145 | fHistMinParSize->UseCurrentStyle();
|
|---|
| 146 | fHistMinParSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 147 | fHistMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
|
|---|
| 148 |
|
|---|
| 149 | fHistDuEnergy = new TH2F("fHistDuEnergy",
|
|---|
| 150 | "Du vs. Energy", 50, 1., 3., 100, -2., 2.);
|
|---|
| 151 | fHistDuEnergy->SetDirectory(NULL);
|
|---|
| 152 | fHistDuEnergy->UseCurrentStyle();
|
|---|
| 153 | fHistDuEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 154 | fHistDuEnergy->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]");
|
|---|
| 155 |
|
|---|
| 156 | fHistDuSize = new TH2F("fHistDuSize",
|
|---|
| 157 | "Du vs. Size", 50, 2., 4., 100, -2., 2.);
|
|---|
| 158 | fHistDuSize->SetDirectory(NULL);
|
|---|
| 159 | fHistDuSize->UseCurrentStyle();
|
|---|
| 160 | fHistDuSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 161 | fHistDuSize->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]");
|
|---|
| 162 |
|
|---|
| 163 | fHistDvEnergy = new TH2F("fHistDvEnergy",
|
|---|
| 164 | "Dv vs. Energy", 50, 1., 3., 100, -2., 2.);
|
|---|
| 165 | fHistDvEnergy->SetDirectory(NULL);
|
|---|
| 166 | fHistDvEnergy->UseCurrentStyle();
|
|---|
| 167 | fHistDvEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 168 | fHistDvEnergy->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]");
|
|---|
| 169 |
|
|---|
| 170 | fHistDvSize = new TH2F("fHistDvSize",
|
|---|
| 171 | "Dv vs. Size", 50, 2., 4., 100, -2., 2.);
|
|---|
| 172 | fHistDvSize->SetDirectory(NULL);
|
|---|
| 173 | fHistDvSize->UseCurrentStyle();
|
|---|
| 174 | fHistDvSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 175 | fHistDvSize->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]");
|
|---|
| 176 |
|
|---|
| 177 | fEvCorrAssign = new TProfile("fEvCorrAssign",
|
|---|
| 178 | "Fraction events srcpos well assign vs. Size", 50, 2., 4., 0., 1.);
|
|---|
| 179 | fEvCorrAssign->SetDirectory(NULL);
|
|---|
| 180 | fEvCorrAssign->SetStats(0);
|
|---|
| 181 | fEvCorrAssign->SetXTitle("log10(Size) [#phot]");
|
|---|
| 182 | fEvCorrAssign->SetYTitle("Fraction of events with srcpos well assign");
|
|---|
| 183 | }
|
|---|
| 184 |
|
|---|
| 185 |
|
|---|
| 186 | // --------------------------------------------------------------------------
|
|---|
| 187 | //
|
|---|
| 188 | // Destructor
|
|---|
| 189 | //
|
|---|
| 190 | MHDisp::~MHDisp()
|
|---|
| 191 | {
|
|---|
| 192 | delete fHistEnergy;
|
|---|
| 193 | delete fHistSize;
|
|---|
| 194 | delete fHistcosZA;
|
|---|
| 195 | delete fSkymapXY;
|
|---|
| 196 | delete fHistMinPar;
|
|---|
| 197 | delete fHistDuDv;
|
|---|
| 198 | delete fHistMinParEnergy;
|
|---|
| 199 | delete fHistMinParSize;
|
|---|
| 200 | delete fHistDuEnergy;
|
|---|
| 201 | delete fHistDuSize;
|
|---|
| 202 | delete fHistDvEnergy;
|
|---|
| 203 | delete fHistDvSize;
|
|---|
| 204 | delete fEvCorrAssign;
|
|---|
| 205 | }
|
|---|
| 206 |
|
|---|
| 207 | // --------------------------------------------------------------------------
|
|---|
| 208 | //
|
|---|
| 209 | // Set the pointers to the containers
|
|---|
| 210 | //
|
|---|
| 211 | //
|
|---|
| 212 | Bool_t MHDisp::SetupFill(const MParList *pList)
|
|---|
| 213 | {
|
|---|
| 214 | // reset all histograms and Minimization parameter computing variables
|
|---|
| 215 | // before each new eventloop
|
|---|
| 216 | fNumEv = 0;
|
|---|
| 217 | fSumMinPar = 0.;
|
|---|
| 218 | fMinPar = (MParameterD*)const_cast<MParList*>(pList)->FindCreateObj("MParameterD", "MinimizationParameter");
|
|---|
| 219 | if (!fMinPar)
|
|---|
| 220 | {
|
|---|
| 221 | *fLog << err << "MParameterD (MinimizationParameter) not found and could not be created... aborting."
|
|---|
| 222 | << endl;
|
|---|
| 223 | return kFALSE;
|
|---|
| 224 | }
|
|---|
| 225 | fMinPar->SetVal(0);
|
|---|
| 226 |
|
|---|
| 227 | fHistEnergy->Reset();
|
|---|
| 228 | fHistSize->Reset();
|
|---|
| 229 | fHistcosZA->Reset();
|
|---|
| 230 | fSkymapXY->Reset();
|
|---|
| 231 | fHistMinPar->Reset();
|
|---|
| 232 | fHistDuDv->Reset();
|
|---|
| 233 | fHistMinParEnergy->Reset();
|
|---|
| 234 | fHistMinParSize->Reset();
|
|---|
| 235 | fHistDuEnergy->Reset();
|
|---|
| 236 | fHistDuSize->Reset();
|
|---|
| 237 | fHistDvEnergy->Reset();
|
|---|
| 238 | fHistDvSize->Reset();
|
|---|
| 239 | fEvCorrAssign->Reset();
|
|---|
| 240 |
|
|---|
| 241 |
|
|---|
| 242 | // look for the defined camera geometry to get mm to deg conversion factor
|
|---|
| 243 | MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
|
|---|
| 244 | if (!cam)
|
|---|
| 245 | {
|
|---|
| 246 | *fLog << err << "MGeomCam (Camera Geometry) not found... aborting."
|
|---|
| 247 | << endl;
|
|---|
| 248 | return kFALSE;
|
|---|
| 249 | }
|
|---|
| 250 | fMm2Deg = cam->GetConvMm2Deg();
|
|---|
| 251 |
|
|---|
| 252 |
|
|---|
| 253 | // look for Disp related containers
|
|---|
| 254 | fImageParDisp = (MImageParDisp*)pList->FindObject(fImageParDispName,
|
|---|
| 255 | "MImageParDisp");
|
|---|
| 256 | if (!fImageParDisp)
|
|---|
| 257 | {
|
|---|
| 258 | *fLog << err << fImageParDispName
|
|---|
| 259 | << " [MImageParDisp] not found... aborting." << endl;
|
|---|
| 260 | return kFALSE;
|
|---|
| 261 | }
|
|---|
| 262 |
|
|---|
| 263 |
|
|---|
| 264 | //-----------------------------------------------------------
|
|---|
| 265 | // if fMatrix exists, skip preprocessing and just read events from matrix;
|
|---|
| 266 | // if doesn't exist, read variables values from containers, so look for them
|
|---|
| 267 | if (fMatrix)
|
|---|
| 268 | return kTRUE;
|
|---|
| 269 |
|
|---|
| 270 | fSrcPos = (MSrcPosCam*)pList->FindObject("MSrcPosCam");
|
|---|
| 271 | if (!fSrcPos)
|
|---|
| 272 | {
|
|---|
| 273 | *fLog << err << "MSrcPosCam not found... aborting." << endl;
|
|---|
| 274 | return kFALSE;
|
|---|
| 275 | }
|
|---|
| 276 |
|
|---|
| 277 | fHil = (MHillas*)pList->FindObject("MHillas");
|
|---|
| 278 | if (!fHil)
|
|---|
| 279 | {
|
|---|
| 280 | *fLog << err << "MHillas not found... aborting." << endl;
|
|---|
| 281 | return kFALSE;
|
|---|
| 282 | }
|
|---|
| 283 |
|
|---|
| 284 | fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
|
|---|
| 285 | if (!fHilExt)
|
|---|
| 286 | {
|
|---|
| 287 | *fLog << err << "MHillasExt not found... aborting." << endl;
|
|---|
| 288 | return kFALSE;
|
|---|
| 289 | }
|
|---|
| 290 |
|
|---|
| 291 | fNewPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
|
|---|
| 292 | if (!fNewPar)
|
|---|
| 293 | {
|
|---|
| 294 | *fLog << err << "MNewImagePar not found... aborting." << endl;
|
|---|
| 295 | return kFALSE;
|
|---|
| 296 | }
|
|---|
| 297 |
|
|---|
| 298 | fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
|
|---|
| 299 | if (!fMcEvt)
|
|---|
| 300 | {
|
|---|
| 301 | *fLog << err << "MMcEvt not found... This is not a MC file,"
|
|---|
| 302 | << " you are not trying to optimize Disp, just calculating it."
|
|---|
| 303 | << endl;
|
|---|
| 304 | // return kFALSE;
|
|---|
| 305 | }
|
|---|
| 306 |
|
|---|
| 307 | fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
|
|---|
| 308 | if (!fPointing)
|
|---|
| 309 | {
|
|---|
| 310 | *fLog << err << "MPointingPos not found... aborting." << endl;
|
|---|
| 311 | return kFALSE;
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | //------------------------------------------
|
|---|
| 315 |
|
|---|
| 316 | return kTRUE;
|
|---|
| 317 | }
|
|---|
| 318 |
|
|---|
| 319 |
|
|---|
| 320 | // --------------------------------------------------------------------------
|
|---|
| 321 | //
|
|---|
| 322 | // Set which selection algorithm for the Disp estimated source position
|
|---|
| 323 | // solutions we want to follow when filling the histograms:
|
|---|
| 324 | // * iflag == 1 --> Correct source position, the closest to the real srcpos
|
|---|
| 325 | // (only applicable to MC or well known source real data)
|
|---|
| 326 | // * iflag == 2 --> Wrong source position, the furthest to the real srcpos
|
|---|
| 327 | // (same applicability than case 1)
|
|---|
| 328 | // * iflag == 3 --> Source position assigned as M3Long sign indicates
|
|---|
| 329 | // * iflag == 4 --> Source position assigned as Asym sign indicates
|
|---|
| 330 | //
|
|---|
| 331 | void MHDisp::SetSelectedPos(Int_t iflag)
|
|---|
| 332 | {
|
|---|
| 333 | fSelectedPos = iflag;
|
|---|
| 334 | }
|
|---|
| 335 |
|
|---|
| 336 |
|
|---|
| 337 | // --------------------------------------------------------------------------
|
|---|
| 338 | //
|
|---|
| 339 | // Sets the Matrix and the array of mapped values for each Matrix column
|
|---|
| 340 | //
|
|---|
| 341 | void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map)
|
|---|
| 342 | {
|
|---|
| 343 | fMatrix = matrix;
|
|---|
| 344 | fMap = map;
|
|---|
| 345 | }
|
|---|
| 346 |
|
|---|
| 347 |
|
|---|
| 348 | // --------------------------------------------------------------------------
|
|---|
| 349 | //
|
|---|
| 350 | // Returns the Matrix and the mapped value for each Matrix column
|
|---|
| 351 | //
|
|---|
| 352 | void MHDisp::GetMatrixMap(MHMatrix* &matrix, TArrayI &map)
|
|---|
| 353 | {
|
|---|
| 354 | map = fMap;
|
|---|
| 355 | matrix = fMatrix;
|
|---|
| 356 | }
|
|---|
| 357 |
|
|---|
| 358 |
|
|---|
| 359 | // --------------------------------------------------------------------------
|
|---|
| 360 | //
|
|---|
| 361 | // You can use this function if you want to use a MHMatrix instead of the
|
|---|
| 362 | // given containers for the Disp optimization. This function adds all
|
|---|
| 363 | // necessary columns to the given matrix. Afterwards, you should fill
|
|---|
| 364 | // the matrix with the corresponding data (eg from a file by using
|
|---|
| 365 | // MHMatrix::Fill). Then, if you loop through the matrix (eg using
|
|---|
| 366 | // MMatrixLoop), MHDisp::Fill will take the values from the matrix
|
|---|
| 367 | // instead of the containers.
|
|---|
| 368 | //
|
|---|
| 369 | void MHDisp::InitMapping(MHMatrix *mat)
|
|---|
| 370 | {
|
|---|
| 371 | if (fMatrix)
|
|---|
| 372 | return;
|
|---|
| 373 |
|
|---|
| 374 | fMatrix = mat;
|
|---|
| 375 |
|
|---|
| 376 | fMap[kX] = fMatrix->AddColumn("MSrcPosCam.fX");
|
|---|
| 377 | fMap[kY] = fMatrix->AddColumn("MSrcPosCam.fY");
|
|---|
| 378 | fMap[kMeanX] = fMatrix->AddColumn("MHillas.fMeanX");
|
|---|
| 379 | fMap[kMeanY] = fMatrix->AddColumn("MHillas.fMeanY");
|
|---|
| 380 | fMap[kDelta] = fMatrix->AddColumn("MHillas.fDelta");
|
|---|
| 381 |
|
|---|
| 382 | fMap[kSize] = fMatrix->AddColumn("MHillas.fSize");
|
|---|
| 383 |
|
|---|
| 384 | fMap[kM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long");
|
|---|
| 385 | fMap[kAsym] = fMatrix->AddColumn("MHillasExt.fAsym");
|
|---|
| 386 |
|
|---|
| 387 | fMap[kEnergy] = fMatrix->AddColumn("MMcEvt.fEnergy");
|
|---|
| 388 | fMap[kImpact] = fMatrix->AddColumn("MMcEvt.fImpact");
|
|---|
| 389 | fMap[kLongitmax] = fMatrix->AddColumn("MMcEvt.fLongitmax");
|
|---|
| 390 |
|
|---|
| 391 | fMap[kZd] = fMatrix->AddColumn("MPointingPos.fZd");
|
|---|
| 392 | fMap[kAz] = fMatrix->AddColumn("MPointingPos.fAz");
|
|---|
| 393 | }
|
|---|
| 394 |
|
|---|
| 395 |
|
|---|
| 396 | // --------------------------------------------------------------------------
|
|---|
| 397 | //
|
|---|
| 398 | // Returns a mapped value from the Matrix
|
|---|
| 399 | //
|
|---|
| 400 | Double_t MHDisp::GetVal(Int_t i) const
|
|---|
| 401 | {
|
|---|
| 402 | Double_t val = (*fMatrix)[fMap[i]];
|
|---|
| 403 |
|
|---|
| 404 | //*fLog << "MHDisp::GetVal; i, fMatrix, fMap, val = "
|
|---|
| 405 | // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl;
|
|---|
| 406 |
|
|---|
| 407 | return val;
|
|---|
| 408 | }
|
|---|
| 409 |
|
|---|
| 410 |
|
|---|
| 411 | // --------------------------------------------------------------------------
|
|---|
| 412 | //
|
|---|
| 413 | // Fill the histograms and sum up the Minimization paramter outcome
|
|---|
| 414 | // of each processed event
|
|---|
| 415 | //
|
|---|
| 416 | Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
|
|---|
| 417 | {
|
|---|
| 418 | Double_t energy = 0.;
|
|---|
| 419 | Double_t impact = 0.;
|
|---|
| 420 | Double_t xmax = 0.;
|
|---|
| 421 |
|
|---|
| 422 | if ( fMatrix || (!fMatrix && fMcEvt) )
|
|---|
| 423 | {
|
|---|
| 424 | energy = fMatrix ? GetVal(kEnergy) : fMcEvt->GetEnergy();
|
|---|
| 425 | impact = fMatrix ? GetVal(kImpact) : fMcEvt->GetImpact();
|
|---|
| 426 | xmax = fMatrix ? GetVal(kLongitmax) : fMcEvt->GetLongitmax();
|
|---|
| 427 | }
|
|---|
| 428 |
|
|---|
| 429 | Double_t theta = fMatrix ? GetVal(kZd) : fPointing->GetZd();
|
|---|
| 430 | // Double_t phi = fMatrix ? GetVal(kAz) : fPointing->GetAz();
|
|---|
| 431 |
|
|---|
| 432 | Double_t xsrcpos0 = fMatrix ? GetVal(kX) : fSrcPos->GetX();
|
|---|
| 433 | Double_t ysrcpos0 = fMatrix ? GetVal(kY) : fSrcPos->GetY();
|
|---|
| 434 | Double_t xmean0 = fMatrix ? GetVal(kMeanX) : fHil->GetMeanX();
|
|---|
| 435 | Double_t ymean0 = fMatrix ? GetVal(kMeanY) : fHil->GetMeanY();
|
|---|
| 436 | Double_t delta = fMatrix ? GetVal(kDelta) : fHil->GetDelta();
|
|---|
| 437 |
|
|---|
| 438 | Double_t size = fMatrix ? GetVal(kSize) : fHil->GetSize();
|
|---|
| 439 |
|
|---|
| 440 | Double_t m3long = fMatrix ? GetVal(kM3Long) : fHilExt->GetM3Long();
|
|---|
| 441 | Double_t asym = fMatrix ? GetVal(kAsym) : fHilExt->GetAsym();
|
|---|
| 442 |
|
|---|
| 443 | //------------------------------------------
|
|---|
| 444 | // convert to deg
|
|---|
| 445 | Double_t xsrcpos = xsrcpos0 * fMm2Deg;
|
|---|
| 446 | Double_t ysrcpos = ysrcpos0 * fMm2Deg;
|
|---|
| 447 | Double_t xmean = xmean0 * fMm2Deg;
|
|---|
| 448 | Double_t ymean = ymean0 * fMm2Deg;
|
|---|
| 449 |
|
|---|
| 450 | // calculate cosinus of the angle between d and a vectors
|
|---|
| 451 | Double_t a = (xmean-xsrcpos)*cos(delta) + (ymean-ysrcpos)*sin(delta);
|
|---|
| 452 | Double_t b = sqrt( (xmean-xsrcpos)*(xmean-xsrcpos) + (ymean-ysrcpos)*(ymean-ysrcpos) );
|
|---|
| 453 | Double_t cosda = a/b;
|
|---|
| 454 |
|
|---|
| 455 | // sign of cosda
|
|---|
| 456 | Int_t s0 = cosda>0 ? 1 : -1;
|
|---|
| 457 |
|
|---|
| 458 | // get the sign of M3Long and Asym variables
|
|---|
| 459 | Int_t sm3 = m3long>0 ? 1 : -1;
|
|---|
| 460 | Int_t sa = asym>0 ? 1 : -1;
|
|---|
| 461 |
|
|---|
| 462 | // set the sign "s" that will select one Disp source position for each
|
|---|
| 463 | // shower, according to each of the possible algorithms for solution selection:
|
|---|
| 464 | // SelectedPos = 1 means choose the right source position
|
|---|
| 465 | // 2 wrong
|
|---|
| 466 | // 3 the position according to M3Long
|
|---|
| 467 | // 4 the position according to Asym
|
|---|
| 468 | Int_t s = s0;
|
|---|
| 469 | if (fSelectedPos == 1)
|
|---|
| 470 | s = s0;
|
|---|
| 471 | else if (fSelectedPos == 2)
|
|---|
| 472 | s = -s0;
|
|---|
| 473 | else if (fSelectedPos == 3)
|
|---|
| 474 | s = sm3;
|
|---|
| 475 | else if (fSelectedPos == 4)
|
|---|
| 476 | s = sa;
|
|---|
| 477 | else
|
|---|
| 478 | *fLog << "Illegal value for Disp srcpos selection algorithm: "
|
|---|
| 479 | << " fSelectedPos = " << fSelectedPos << endl;
|
|---|
| 480 |
|
|---|
| 481 | // count the number of events where source position has been correctly assigned
|
|---|
| 482 | if (s == s0)
|
|---|
| 483 | fEvCorrAssign->Fill(log10(size), 1.);
|
|---|
| 484 | else
|
|---|
| 485 | fEvCorrAssign->Fill(log10(size), 0.);
|
|---|
| 486 |
|
|---|
| 487 | // get estimated Disp value
|
|---|
| 488 | Double_t disp = fImageParDisp->GetDisp();
|
|---|
| 489 |
|
|---|
| 490 | //------------------------------------------
|
|---|
| 491 | // Disp estimated source position
|
|---|
| 492 | Double_t xdisp = xmean - s*cos(delta)*disp;
|
|---|
| 493 | Double_t ydisp = ymean - s*sin(delta)*disp;
|
|---|
| 494 | fSkymapXY->Fill(xdisp,ydisp);
|
|---|
| 495 |
|
|---|
| 496 | // Distance between estimated Disp and real source position
|
|---|
| 497 | Double_t d2 = (xdisp-xsrcpos)*(xdisp-xsrcpos) + (ydisp-ysrcpos)*(ydisp-ysrcpos);
|
|---|
| 498 | fHistMinPar->Fill(d2);
|
|---|
| 499 |
|
|---|
| 500 | // Longitudinal and transversal distances between Disp and real source positon
|
|---|
| 501 | Double_t Du = -s*( (xdisp-xsrcpos)*cos(delta) + (ydisp-ysrcpos)*sin(delta));
|
|---|
| 502 | Double_t Dv = -s*(-(xdisp-xsrcpos)*sin(delta) + (ydisp-ysrcpos)*cos(delta));
|
|---|
| 503 | fHistDuDv->Fill(Dv,Du);
|
|---|
| 504 |
|
|---|
| 505 | // Fill Energy, Size and ZA distributions
|
|---|
| 506 | if (fMatrix || (!fMatrix && fMcEvt))
|
|---|
| 507 | fHistEnergy->Fill(log10(energy));
|
|---|
| 508 | fHistSize->Fill(log10(size));
|
|---|
| 509 | fHistcosZA->Fill(cos(theta/kRad2Deg));
|
|---|
| 510 |
|
|---|
| 511 | // to check the size and energy dependence of the optimization
|
|---|
| 512 | if (fMatrix || (!fMatrix && fMcEvt))
|
|---|
| 513 | {
|
|---|
| 514 | fHistMinParEnergy->Fill(log10(energy),d2);
|
|---|
| 515 | fHistDuEnergy->Fill(log10(energy),Du);
|
|---|
| 516 | fHistDvEnergy->Fill(log10(energy),Dv);
|
|---|
| 517 | }
|
|---|
| 518 | fHistMinParSize->Fill(log10(size),d2);
|
|---|
| 519 | fHistDuSize->Fill(log10(size),Du);
|
|---|
| 520 | fHistDvSize->Fill(log10(size),Dv);
|
|---|
| 521 |
|
|---|
| 522 | // variables for the Minimization parameter calculation (= d^2)
|
|---|
| 523 | fNumEv += 1;
|
|---|
| 524 | fSumMinPar += d2;
|
|---|
| 525 |
|
|---|
| 526 | return kTRUE;
|
|---|
| 527 | }
|
|---|
| 528 |
|
|---|
| 529 |
|
|---|
| 530 | // --------------------------------------------------------------------------
|
|---|
| 531 | //
|
|---|
| 532 | // Calculates the final Minimization parameter of the Disp optimization
|
|---|
| 533 | //
|
|---|
| 534 | Bool_t MHDisp::Finalize()
|
|---|
| 535 | {
|
|---|
| 536 | fMinPar->SetVal(fSumMinPar/fNumEv);
|
|---|
| 537 | *fLog << endl;
|
|---|
| 538 | *fLog << "MHDisp::Finalize: SumMinPar, NumEv = " << fSumMinPar << ", " << fNumEv << endl;
|
|---|
| 539 | *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl;
|
|---|
| 540 |
|
|---|
| 541 | fMinPar->SetVal(fHistMinPar->GetMean());
|
|---|
| 542 | *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl;
|
|---|
| 543 |
|
|---|
| 544 | SetReadyToSave();
|
|---|
| 545 |
|
|---|
| 546 | return kTRUE;
|
|---|
| 547 | }
|
|---|
| 548 |
|
|---|
| 549 |
|
|---|
| 550 | // --------------------------------------------------------------------------
|
|---|
| 551 | //
|
|---|
| 552 | // Creates a new canvas and draws the Disp related histograms into it.
|
|---|
| 553 | // Be careful: The histograms belongs to this object and won't get deleted
|
|---|
| 554 | // together with the canvas.
|
|---|
| 555 | //
|
|---|
| 556 | void MHDisp::Draw(Option_t *opt)
|
|---|
| 557 | {
|
|---|
| 558 | gStyle->SetPalette(1);
|
|---|
| 559 |
|
|---|
| 560 | TString title = GetName();
|
|---|
| 561 | title += ": ";
|
|---|
| 562 | title += GetTitle();
|
|---|
| 563 |
|
|---|
| 564 | TCanvas *pad = new TCanvas(GetName(),title,0,0,900,1500);
|
|---|
| 565 | pad->SetBorderMode(0);
|
|---|
| 566 | pad->Divide(3,5);
|
|---|
| 567 |
|
|---|
| 568 | pad->cd(1);
|
|---|
| 569 | gPad->SetBorderMode(0);
|
|---|
| 570 | fHistcosZA->SetTitleOffset(1.2,"Y");
|
|---|
| 571 | fHistcosZA->DrawClone(opt);
|
|---|
| 572 | fHistcosZA->SetBit(kCanDelete);
|
|---|
| 573 | gPad->Modified();
|
|---|
| 574 |
|
|---|
| 575 | pad->cd(2);
|
|---|
| 576 | gPad->SetBorderMode(0);
|
|---|
| 577 | fHistEnergy->SetTitleOffset(1.2,"Y");
|
|---|
| 578 | fHistEnergy->DrawClone(opt);
|
|---|
| 579 | fHistEnergy->SetBit(kCanDelete);
|
|---|
| 580 | gPad->Modified();
|
|---|
| 581 |
|
|---|
| 582 | pad->cd(3);
|
|---|
| 583 | gPad->SetBorderMode(0);
|
|---|
| 584 | fHistSize->SetTitleOffset(1.2,"Y");
|
|---|
| 585 | fHistSize->DrawClone(opt);
|
|---|
| 586 | fHistSize->SetBit(kCanDelete);
|
|---|
| 587 | gPad->Modified();
|
|---|
| 588 |
|
|---|
| 589 | pad->cd(4);
|
|---|
| 590 | gPad->SetBorderMode(0);
|
|---|
| 591 | fHistMinPar->SetTitleOffset(1.2,"Y");
|
|---|
| 592 | fHistMinPar->DrawClone(opt);
|
|---|
| 593 | fHistMinPar->SetBit(kCanDelete);
|
|---|
| 594 | gPad->Modified();
|
|---|
| 595 |
|
|---|
| 596 | TProfile *profMinParEnergy = fHistMinParEnergy->ProfileX("Minimization parameter vs. Energy",0,99999,"s");
|
|---|
| 597 | profMinParEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 598 | profMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
|
|---|
| 599 | pad->cd(5);
|
|---|
| 600 | gPad->SetBorderMode(0);
|
|---|
| 601 | profMinParEnergy->SetTitleOffset(1.2,"Y");
|
|---|
| 602 | profMinParEnergy->SetStats(0);
|
|---|
| 603 | profMinParEnergy->DrawClone(opt);
|
|---|
| 604 | profMinParEnergy->SetBit(kCanDelete);
|
|---|
| 605 | gPad->Modified();
|
|---|
| 606 |
|
|---|
| 607 | TProfile *profMinParSize = fHistMinParSize->ProfileX("Minimization parameter vs. Size",0,99999,"s");
|
|---|
| 608 | profMinParSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 609 | profMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
|
|---|
| 610 | pad->cd(6);
|
|---|
| 611 | gPad->SetBorderMode(0);
|
|---|
| 612 | profMinParSize->SetTitleOffset(1.2,"Y");
|
|---|
| 613 | profMinParSize->SetStats(0);
|
|---|
| 614 | profMinParSize->DrawClone(opt);
|
|---|
| 615 | profMinParSize->SetBit(kCanDelete);
|
|---|
| 616 | gPad->Modified();
|
|---|
| 617 |
|
|---|
| 618 | pad->cd(7);
|
|---|
| 619 | gPad->SetBorderMode(0);
|
|---|
| 620 | fSkymapXY->SetTitleOffset(1.2,"Y");
|
|---|
| 621 | fSkymapXY->DrawClone("COLZ");
|
|---|
| 622 | fSkymapXY->SetBit(kCanDelete);
|
|---|
| 623 | gPad->Modified();
|
|---|
| 624 |
|
|---|
| 625 | pad->cd(8);
|
|---|
| 626 | gPad->SetBorderMode(0);
|
|---|
| 627 | fEvCorrAssign->SetTitleOffset(1.2,"Y");
|
|---|
| 628 | fEvCorrAssign->DrawClone(opt);
|
|---|
| 629 | fEvCorrAssign->SetBit(kCanDelete);
|
|---|
| 630 | gPad->Modified();
|
|---|
| 631 |
|
|---|
| 632 | pad->cd(9);
|
|---|
| 633 | gPad->SetBorderMode(0);
|
|---|
| 634 | fHistDuDv->SetTitleOffset(1.2,"Y");
|
|---|
| 635 | fHistDuDv->DrawClone("COLZ");
|
|---|
| 636 | fHistDuDv->SetBit(kCanDelete);
|
|---|
| 637 | gPad->Modified();
|
|---|
| 638 |
|
|---|
| 639 | TH1F *histDu = (TH1F*)fHistDuDv->ProjectionY("histDu");
|
|---|
| 640 | histDu->SetTitle("Longitudinal distance Du");
|
|---|
| 641 | histDu->SetXTitle("Du = longitudinal distance [deg]");
|
|---|
| 642 | pad->cd(10);
|
|---|
| 643 | gPad->SetBorderMode(0);
|
|---|
| 644 | histDu->SetTitleOffset(1.2,"Y");
|
|---|
| 645 | histDu->DrawClone(opt);
|
|---|
| 646 | histDu->SetBit(kCanDelete);
|
|---|
| 647 | gPad->Modified();
|
|---|
| 648 |
|
|---|
| 649 | TProfile *profDuEnergy = fHistDuEnergy->ProfileX("Du vs. Energy",0,99999,"s");
|
|---|
| 650 | profDuEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 651 | profDuEnergy->SetYTitle("Du = longitudinal distance [deg]");
|
|---|
| 652 | pad->cd(11);
|
|---|
| 653 | gPad->SetBorderMode(0);
|
|---|
| 654 | profDuEnergy->SetTitleOffset(1.2,"Y");
|
|---|
| 655 | profDuEnergy->SetStats(0);
|
|---|
| 656 | profDuEnergy->DrawClone(opt);
|
|---|
| 657 | profDuEnergy->SetBit(kCanDelete);
|
|---|
| 658 | gPad->Modified();
|
|---|
| 659 |
|
|---|
| 660 | TProfile *profDuSize = fHistDuSize->ProfileX("Du vs. Size",0,99999,"s");
|
|---|
| 661 | profDuSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 662 | profDuSize->SetYTitle("Du = longitudinal distance [deg]");
|
|---|
| 663 | pad->cd(12);
|
|---|
| 664 | gPad->SetBorderMode(0);
|
|---|
| 665 | profDuSize->SetTitleOffset(1.2,"Y");
|
|---|
| 666 | profDuSize->SetStats(0);
|
|---|
| 667 | profDuSize->DrawClone(opt);
|
|---|
| 668 | profDuSize->SetBit(kCanDelete);
|
|---|
| 669 | gPad->Modified();
|
|---|
| 670 |
|
|---|
| 671 | TH1F *histDv = (TH1F*)fHistDuDv->ProjectionX("histDv");
|
|---|
| 672 | histDv->SetTitle("Transversal distance Dv");
|
|---|
| 673 | histDv->SetXTitle("Dv = transversal distance [deg]");
|
|---|
| 674 | pad->cd(13);
|
|---|
| 675 | gPad->SetBorderMode(0);
|
|---|
| 676 | histDv->SetTitleOffset(1.2,"Y");
|
|---|
| 677 | histDv->DrawClone(opt);
|
|---|
| 678 | histDv->SetBit(kCanDelete);
|
|---|
| 679 | gPad->Modified();
|
|---|
| 680 |
|
|---|
| 681 | TProfile *profDvEnergy = fHistDvEnergy->ProfileX("Dv vs. Energy",0,99999,"s");
|
|---|
| 682 | profDvEnergy->SetXTitle("log10(Energy) [GeV]");
|
|---|
| 683 | profDvEnergy->SetYTitle("Dv = transversal distance [deg]");
|
|---|
| 684 | pad->cd(14);
|
|---|
| 685 | gPad->SetBorderMode(0);
|
|---|
| 686 | profDvEnergy->SetTitleOffset(1.2,"Y");
|
|---|
| 687 | profDvEnergy->SetStats(0);
|
|---|
| 688 | profDvEnergy->DrawClone(opt);
|
|---|
| 689 | profDvEnergy->SetBit(kCanDelete);
|
|---|
| 690 | gPad->Modified();
|
|---|
| 691 |
|
|---|
| 692 | TProfile *profDvSize = fHistDvSize->ProfileX("Dv vs. Size",0,99999,"s");
|
|---|
| 693 | profDvSize->SetXTitle("log10(Size) [#phot]");
|
|---|
| 694 | profDvSize->SetYTitle("Dv = transversal distance [deg]");
|
|---|
| 695 | pad->cd(15);
|
|---|
| 696 | gPad->SetBorderMode(0);
|
|---|
| 697 | profDvSize->SetTitleOffset(1.2,"Y");
|
|---|
| 698 | profDvSize->SetStats(0);
|
|---|
| 699 | profDvSize->DrawClone(opt);
|
|---|
| 700 | profDvSize->SetBit(kCanDelete);
|
|---|
| 701 | gPad->Modified();
|
|---|
| 702 | }
|
|---|
| 703 |
|
|---|
| 704 |
|
|---|
| 705 |
|
|---|
| 706 |
|
|---|
| 707 |
|
|---|
| 708 |
|
|---|
| 709 |
|
|---|
| 710 |
|
|---|
| 711 |
|
|---|
| 712 |
|
|---|
| 713 |
|
|---|
| 714 |
|
|---|
| 715 |
|
|---|
| 716 |
|
|---|
| 717 |
|
|---|