| 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): Thomas Bretz 12/2000 <mailto:tbretz@uni-sw.gwdg.de>
|
|---|
| 19 | ! Author(s): Harald Kornmayer 1/2001
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2001
|
|---|
| 22 | !
|
|---|
| 23 | ! Modified 4/7/2002 Abelardo Moralejo: now the dimension of fTrigger is
|
|---|
| 24 | ! set dinamically, to allow an arbitrary large number of trigger
|
|---|
| 25 | ! conditions to be processed.
|
|---|
| 26 | !
|
|---|
| 27 | !
|
|---|
| 28 | \* ======================================================================== */
|
|---|
| 29 |
|
|---|
| 30 | #include "MMcTriggerRateCalc.h"
|
|---|
| 31 |
|
|---|
| 32 | #include "MLog.h"
|
|---|
| 33 | #include "MLogManip.h"
|
|---|
| 34 |
|
|---|
| 35 | #include "MParList.h"
|
|---|
| 36 | #include "MHMcRate.h"
|
|---|
| 37 |
|
|---|
| 38 | #include "MMcEvt.hxx"
|
|---|
| 39 | #include "MMcTrig.hxx"
|
|---|
| 40 | #include "MMcRunHeader.hxx"
|
|---|
| 41 | #include "MMcCorsikaRunHeader.h"
|
|---|
| 42 |
|
|---|
| 43 | ClassImp(MMcTriggerRateCalc);
|
|---|
| 44 |
|
|---|
| 45 | void MMcTriggerRateCalc::Init(int dim, float *trigbg, float simbg,
|
|---|
| 46 | const char *name, const char *title)
|
|---|
| 47 | {
|
|---|
| 48 | fName = name ? name : "MMcTriggerRateCalc";
|
|---|
| 49 | fTitle = title ? title : "Task to calc the trigger rate ";
|
|---|
| 50 |
|
|---|
| 51 | fMcTrig = NULL;
|
|---|
| 52 | fMcRate = NULL;
|
|---|
| 53 |
|
|---|
| 54 | fTrigNSB = trigbg;
|
|---|
| 55 | fSimNSB = simbg;
|
|---|
| 56 |
|
|---|
| 57 | fPartId = -1;
|
|---|
| 58 |
|
|---|
| 59 | fShowers = 0;
|
|---|
| 60 | fAnalShow = 0;
|
|---|
| 61 |
|
|---|
| 62 | fFirst = dim>0 ? 1 : -dim;
|
|---|
| 63 | fLast = dim>0 ? dim : -dim;
|
|---|
| 64 |
|
|---|
| 65 | fNum = fLast-fFirst+1;
|
|---|
| 66 | fTrigger = new float[fNum];
|
|---|
| 67 |
|
|---|
| 68 | for (UInt_t i=0;i<fNum;i++)
|
|---|
| 69 | fTrigger[i]=0;
|
|---|
| 70 |
|
|---|
| 71 | AddToBranchList("MMcEvt.fPartId");
|
|---|
| 72 | AddToBranchList("MMcEvt.fImpact");
|
|---|
| 73 | AddToBranchList("MMcEvt.fEnergy");
|
|---|
| 74 | AddToBranchList("MMcEvt.fPhi");
|
|---|
| 75 | AddToBranchList("MMcEvt.fTheta");
|
|---|
| 76 | AddToBranchList("MMcEvt.fPhotElfromShower");
|
|---|
| 77 | AddToBranchList("MMcTrig", "fNumFirstLevel", fFirst, fLast);
|
|---|
| 78 |
|
|---|
| 79 | }
|
|---|
| 80 |
|
|---|
| 81 | // ReInit: read run header to get some info later:
|
|---|
| 82 |
|
|---|
| 83 | Bool_t MMcTriggerRateCalc::ReInit(MParList *pList)
|
|---|
| 84 | {
|
|---|
| 85 | fMcRunHeader = (MMcRunHeader*) pList->FindObject("MMcRunHeader");
|
|---|
| 86 | if (!fMcRunHeader)
|
|---|
| 87 | {
|
|---|
| 88 | *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
|
|---|
| 89 | return kFALSE;
|
|---|
| 90 | }
|
|---|
| 91 |
|
|---|
| 92 | fMcCorRunHeader = (MMcCorsikaRunHeader*) pList->FindObject("MMcCorsikaRunHeader");
|
|---|
| 93 | if (!fMcCorRunHeader)
|
|---|
| 94 | {
|
|---|
| 95 | *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl;
|
|---|
| 96 | return kFALSE;
|
|---|
| 97 | }
|
|---|
| 98 |
|
|---|
| 99 | if (fMcRunHeader->GetAllEvtsTriggered())
|
|---|
| 100 | fShowers += fMcRunHeader->GetNumSimulatedShowers();
|
|---|
| 101 |
|
|---|
| 102 | if (fMcRunHeader->GetAllEvtsTriggered())
|
|---|
| 103 | {
|
|---|
| 104 | *fLog << endl << all << endl <<
|
|---|
| 105 | "WARNING: the input data file contains only the" << endl <<
|
|---|
| 106 | "events that triggered. I will assume the standard value" << endl <<
|
|---|
| 107 | "for maximum impact parameter (400 m)" <<endl;
|
|---|
| 108 |
|
|---|
| 109 |
|
|---|
| 110 | if (fTrigNSB[0] > 0)
|
|---|
| 111 | *fLog << endl << all <<
|
|---|
| 112 | "WARNING: NSB rate can be overestimated by up to 5%." << endl <<
|
|---|
| 113 | "For a precise estimate of the total rate including NSB" << endl <<
|
|---|
| 114 | "accidental triggers I need a file containing all event headers."
|
|---|
| 115 | << endl;
|
|---|
| 116 | else
|
|---|
| 117 | *fLog << endl << all <<
|
|---|
| 118 | "WARNING: calculating only shower rate (no NSB accidental triggers)" << endl;
|
|---|
| 119 | }
|
|---|
| 120 |
|
|---|
| 121 | *fLog << endl << all <<
|
|---|
| 122 | "WARNING: I will assume the standard maximum off axis angle" << endl <<
|
|---|
| 123 | "(5 degrees) for the original showers" << endl;
|
|---|
| 124 |
|
|---|
| 125 | for (UInt_t i=0; i<fNum; i++)
|
|---|
| 126 | {
|
|---|
| 127 | if (fMcRunHeader->GetAllEvtsTriggered())
|
|---|
| 128 | {
|
|---|
| 129 | GetRate(i)->SetImpactMin(0.);
|
|---|
| 130 | GetRate(i)->SetImpactMax(40000.); // in cm
|
|---|
| 131 | }
|
|---|
| 132 | GetRate(i)->SetSolidAngle(2.390941702e-2); // sr
|
|---|
| 133 |
|
|---|
| 134 | //
|
|---|
| 135 | // Set limits of energy range, coverting from GeV to TeV:
|
|---|
| 136 | //
|
|---|
| 137 | GetRate(i)->SetEnergyMin(fMcCorRunHeader->GetELowLim()/1000.);
|
|---|
| 138 | GetRate(i)->SetEnergyMax(fMcCorRunHeader->GetEUppLim()/1000.);
|
|---|
| 139 |
|
|---|
| 140 | }
|
|---|
| 141 |
|
|---|
| 142 | return kTRUE;
|
|---|
| 143 | }
|
|---|
| 144 |
|
|---|
| 145 | // --------------------------------------------------------------------------
|
|---|
| 146 | //
|
|---|
| 147 | // overloaded constructor I
|
|---|
| 148 | //
|
|---|
| 149 | // dim: fDimension
|
|---|
| 150 | // *trigbg: number of NSB triggers for
|
|---|
| 151 | // a given trigger condition.
|
|---|
| 152 | // simbg: Number of simulated events for the NSB
|
|---|
| 153 | // rate: rate of incident showers
|
|---|
| 154 | //
|
|---|
| 155 |
|
|---|
| 156 | MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim,
|
|---|
| 157 | float *trigbg, float simbg,
|
|---|
| 158 | const char *name, const char *title)
|
|---|
| 159 | {
|
|---|
| 160 | Init(dim, trigbg, simbg, name, title);
|
|---|
| 161 | }
|
|---|
| 162 |
|
|---|
| 163 |
|
|---|
| 164 | // --------------------------------------------------------------------------
|
|---|
| 165 | //
|
|---|
| 166 | // overloaded constructor II
|
|---|
| 167 | //
|
|---|
| 168 | // dim: fDimension
|
|---|
| 169 | // *trigbg: number of NSB triggers for
|
|---|
| 170 | // a given trigger condition.
|
|---|
| 171 | // simbg: Number of simulated events for the NSB
|
|---|
| 172 | //
|
|---|
| 173 | MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, float *trigbg,float simbg,
|
|---|
| 174 | const char *name, const char *title)
|
|---|
| 175 | {
|
|---|
| 176 | Init(dim, trigbg, simbg, name, title);
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 | MMcTriggerRateCalc::~MMcTriggerRateCalc()
|
|---|
| 180 | {
|
|---|
| 181 | if (fMcTrig)
|
|---|
| 182 | delete fMcTrig;
|
|---|
| 183 |
|
|---|
| 184 | if (fMcRate)
|
|---|
| 185 | delete fMcRate;
|
|---|
| 186 | }
|
|---|
| 187 |
|
|---|
| 188 |
|
|---|
| 189 | // --------------------------------------------------------------------------
|
|---|
| 190 | //
|
|---|
| 191 | // The PreProcess connects the raw data with this task. It checks if the
|
|---|
| 192 | // input containers exist, if not a kFalse flag is returned. It also checks
|
|---|
| 193 | // if the output contaniers exist, if not they are created.
|
|---|
| 194 | // This task can read either Montecarlo files with multiple trigger
|
|---|
| 195 | // options, either Montecarlo files with a single trigger option.
|
|---|
| 196 | //
|
|---|
| 197 | Bool_t MMcTriggerRateCalc::PreProcess (MParList *pList)
|
|---|
| 198 | {
|
|---|
| 199 | // connect the raw data with this task
|
|---|
| 200 |
|
|---|
| 201 | fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
|
|---|
| 202 | if (!fMcEvt)
|
|---|
| 203 | {
|
|---|
| 204 | *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
|
|---|
| 205 | return kFALSE;
|
|---|
| 206 | }
|
|---|
| 207 |
|
|---|
| 208 | UInt_t num;
|
|---|
| 209 |
|
|---|
| 210 | fMcTrig = new TObjArray(pList->FindObjectList("MMcTrig", fFirst, fLast));
|
|---|
| 211 | num = fMcTrig->GetEntriesFast();
|
|---|
| 212 | if (num != fNum)
|
|---|
| 213 | {
|
|---|
| 214 | *fLog << err << dbginf << fNum << " MMcTrig objects requested, ";
|
|---|
| 215 | *fLog << num << " are available... aborting." << endl;
|
|---|
| 216 | return kFALSE;
|
|---|
| 217 | }
|
|---|
| 218 |
|
|---|
| 219 | fMcRate = new TObjArray(pList->FindObjectList("MHMcRate", fFirst, fLast));
|
|---|
| 220 | num = fMcRate->GetEntriesFast();
|
|---|
| 221 | if (num != fNum)
|
|---|
| 222 | {
|
|---|
| 223 | *fLog << err << dbginf << fNum << " MHMcRate objects requested, ";
|
|---|
| 224 | *fLog << num << " are available... aborting." << endl;
|
|---|
| 225 | return kFALSE;
|
|---|
| 226 | }
|
|---|
| 227 |
|
|---|
| 228 | for (UInt_t i=0;i<fNum;i++)
|
|---|
| 229 | {
|
|---|
| 230 | MHMcRate &rate = *GetRate(i);
|
|---|
| 231 |
|
|---|
| 232 | if (fTrigNSB)
|
|---|
| 233 | rate.SetBackground(fTrigNSB[i], fSimNSB);
|
|---|
| 234 | else
|
|---|
| 235 | rate.SetBackground(0., fSimNSB);
|
|---|
| 236 | }
|
|---|
| 237 |
|
|---|
| 238 | return kTRUE;
|
|---|
| 239 | }
|
|---|
| 240 |
|
|---|
| 241 | // --------------------------------------------------------------------------
|
|---|
| 242 | //
|
|---|
| 243 | // The Process-function counts the number of simulated showers, the
|
|---|
| 244 | // number of analised showers and the number of triggers. It also updates
|
|---|
| 245 | // the limits for theta, phi, energy and impact parameter in the
|
|---|
| 246 | // MHMcRate container.
|
|---|
| 247 | //
|
|---|
| 248 | Bool_t MMcTriggerRateCalc::Process()
|
|---|
| 249 | {
|
|---|
| 250 | //
|
|---|
| 251 | // Counting analysed showers (except in the case in which the file
|
|---|
| 252 | // contains only events that triggered, since then we do not know
|
|---|
| 253 | // how many showers were analysed).
|
|---|
| 254 | //
|
|---|
| 255 |
|
|---|
| 256 | if ( ! fMcRunHeader->GetAllEvtsTriggered() &&
|
|---|
| 257 | fMcEvt->GetPhotElfromShower() )
|
|---|
| 258 | fAnalShow++;
|
|---|
| 259 |
|
|---|
| 260 | //
|
|---|
| 261 | // Set PartId and check it is the same for all events
|
|---|
| 262 | //
|
|---|
| 263 | if (fPartId < 0)
|
|---|
| 264 | fPartId = fMcEvt->GetPartId();
|
|---|
| 265 | else if (fPartId != fMcEvt->GetPartId())
|
|---|
| 266 | {
|
|---|
| 267 | *fLog << err << dbginf << "Incident particle type seems to change";
|
|---|
| 268 | *fLog << "from " << fPartId << " to " << fMcEvt->GetPartId() << endl;
|
|---|
| 269 | *fLog << "in input data files... aborting." << endl;
|
|---|
| 270 | return kFALSE;
|
|---|
| 271 | }
|
|---|
| 272 |
|
|---|
| 273 | //
|
|---|
| 274 | // Getting angles, energy and impact parameter to set boundaries
|
|---|
| 275 | //
|
|---|
| 276 | const Float_t theta = fMcEvt->GetTheta();
|
|---|
| 277 | const Float_t phi = fMcEvt->GetPhi();
|
|---|
| 278 | const Float_t param = fMcEvt->GetImpact();
|
|---|
| 279 | const Float_t ener = fMcEvt->GetEnergy()/1000.0;
|
|---|
| 280 |
|
|---|
| 281 | //
|
|---|
| 282 | // Counting number of triggers
|
|---|
| 283 | //
|
|---|
| 284 | for (UInt_t i=0; i<fNum; i++)
|
|---|
| 285 | {
|
|---|
| 286 | if (GetTrig(i)->GetFirstLevel())
|
|---|
| 287 | fTrigger[i] ++;
|
|---|
| 288 |
|
|---|
| 289 | if ( ! fMcRunHeader->GetAllEvtsTriggered())
|
|---|
| 290 | GetRate(i)->UpdateBoundaries(ener, theta, phi, param);
|
|---|
| 291 | }
|
|---|
| 292 |
|
|---|
| 293 | return kTRUE;
|
|---|
| 294 | }
|
|---|
| 295 |
|
|---|
| 296 | // --------------------------------------------------------------------------
|
|---|
| 297 | //
|
|---|
| 298 | // The PostProcess-function calculates and shows the trigger rate
|
|---|
| 299 | //
|
|---|
| 300 | Bool_t MMcTriggerRateCalc::PostProcess()
|
|---|
| 301 | {
|
|---|
| 302 | for (UInt_t i=0; i<fNum; i++)
|
|---|
| 303 | {
|
|---|
| 304 | MHMcRate &rate = *GetRate(i);
|
|---|
| 305 |
|
|---|
| 306 | rate.SetParticle(fPartId);
|
|---|
| 307 | switch (fPartId)
|
|---|
| 308 | {
|
|---|
| 309 | case kPROTON:
|
|---|
| 310 | if (fMcCorRunHeader->GetSlopeSpec() != -2.75)
|
|---|
| 311 | {
|
|---|
| 312 | *fLog << err << dbginf <<
|
|---|
| 313 | "Spectrum slope as read from input file ("<<
|
|---|
| 314 | fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for protons" << endl << "... aborting." << endl;
|
|---|
| 315 | return kFALSE;
|
|---|
| 316 | }
|
|---|
| 317 | rate.SetFlux(0.1091, 2.75);
|
|---|
| 318 | break;
|
|---|
| 319 | case kHELIUM:
|
|---|
| 320 | if (fMcCorRunHeader->GetSlopeSpec() != -2.62)
|
|---|
| 321 | {
|
|---|
| 322 | *fLog << err << dbginf <<
|
|---|
| 323 | "Spectrum slope as read from input file ("<<
|
|---|
| 324 | fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for Helium" << endl << "... aborting." << endl;
|
|---|
| 325 | return kFALSE;
|
|---|
| 326 | }
|
|---|
| 327 | rate.SetFlux(0.0660, 2.62);
|
|---|
| 328 | break;
|
|---|
| 329 | default:
|
|---|
| 330 | *fLog << err << dbginf << "Unknown incident flux parameters for ";
|
|---|
| 331 | *fLog << fPartId<< " particle Id ... aborting." << endl;
|
|---|
| 332 | return kFALSE;
|
|---|
| 333 | }
|
|---|
| 334 | }
|
|---|
| 335 |
|
|---|
| 336 | //
|
|---|
| 337 | // Computing trigger rate
|
|---|
| 338 | //
|
|---|
| 339 | for (UInt_t i=0; i<fNum; i++)
|
|---|
| 340 | GetRate(i)->CalcRate(fTrigger[i], fAnalShow, fShowers);
|
|---|
| 341 |
|
|---|
| 342 | return kTRUE;
|
|---|
| 343 | }
|
|---|
