- Timestamp:
- 12/10/09 10:08:55 (15 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r9524 r9525 18 18 19 19 -*-*- END OF LINE -*-*- 20 21 2009/12/10 Thomas Bretz 22 23 * ceres.rc: 24 - added line for new option 25 26 * mhflux/MAlphaFitter.cc: 27 - added "Sensitivity" to output 28 - prevent zero or negative divisor 29 30 * msim/MSimMMCS.h: 31 - removed pointers from i/o 32 33 * msimcamera/MSimCamera.h: 34 - fixed a typo in a comment 35 36 * msimcamera/MSimRandomPhotons.[h,cc]: 37 - changed to allow simulation of a nsb spectrum 38 39 20 40 21 41 2009/12/01 Thomas Bretz -
trunk/MagicSoft/Mars/ceres.rc
r9518 r9525 152 152 # Dark Counts per APD: ~4MHz 153 153 #MSimRandomPhotons.FrequencyFixed: 0.004 154 # NSB photon rate [ph/m^2/nm/sr/ns] 154 # Spectrum to use as NSB spectrum (see Benn et al. for example) 155 #MSimRandomPhotons.FileNameNSB: resmc/night-sky-la-palma.txt 156 # If an empty file name is given a constant is used instead: 157 # NSB photon rate [ph/m^2/nm/sr/ns] 155 158 #MSimRandomPhotons.FrequencyNSB: 5.8 156 159 157 160 # FIXME: With a class describing the cones we could give NSB as 158 161 # per sr and cm^2 162 159 163 160 164 # ------------------------------------------------------------------------- -
trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc
r9367 r9525 598 598 *fLog << " - Signal Events " << fEventsSignal << endl; 599 599 *fLog << " - Background Events " << fEventsBackground << endl; 600 *fLog << " - E/sqrt(B>=Alpha) " << fEventsExcess/TMath::Sqrt(TMath::Max(fEventsBackground,fScaleFactor)) << endl; 600 601 *fLog << " - Chi^2/ndf (Signal) " << fChiSqSignal << endl; 601 602 *fLog << " - Chi^2/ndf (Background) " << fChiSqBg << endl; 602 603 *fLog << " - Signal integrated up to " << fIntegralMax << "°" << endl; 603 *fLog << " - Scale Factor (Off)" << fScaleFactor << endl;604 *fLog << " - Off Scale Alpha (Off) " << fScaleFactor << endl; 604 605 } 605 606 } … … 858 859 return GetGausSigma(); 859 860 case kWeakSource: 860 return GetEventsBackground()<1 ? -GetEventsExcess() : -GetEventsExcess()/TMath::Sqrt(GetEventsBackground()); 861 if (GetEventsExcess()<1) 862 return 0; 863 return -GetEventsExcess()/TMath::Sqrt(TMath::Max(GetEventsBackground(), GetEventsBackground())); 861 864 case kWeakSourceLogExcess: 862 return GetEventsBackground()<1 ? -GetEventsExcess() : -GetEventsExcess()/TMath::Sqrt(GetEventsBackground())*TMath::Log10(GetEventsExcess()); 865 if (GetEventsExcess()<1) 866 return 0; 867 return -GetEventsExcess()/TMath::Sqrt(TMath::Max(GetEventsBackground(), GetEventsBackground()))*TMath::Log10(GetEventsExcess()); 863 868 } 864 869 return 0; -
trunk/MagicSoft/Mars/mjobs/MJSimulation.cc
r9482 r9525 368 368 MParSpline splinemirror("MirrorReflectivity"); 369 369 MParSpline splinecones("ConesAngularAcceptance"); 370 MParSpline splinecones2("ConesTransmission"); 370 371 plist.AddToList(&splinepde); 371 372 plist.AddToList(&splinemirror); 372 373 plist.AddToList(&splinecones); 374 plist.AddToList(&splinecones2); 373 375 374 376 const TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())"; … … 385 387 MSimAbsorption absmir("SimMirrorReflectivity"); 386 388 MSimAbsorption cones("SimConesAngularAcceptance"); 389 MSimAbsorption cones2("SimConesTransmission"); 387 390 absapd.SetParName("PhotonDetectionEfficiency"); 388 391 absmir.SetParName("MirrorReflectivity"); 389 392 cones.SetParName("ConesAngularAcceptance"); 390 393 cones.SetUseTheta(); 391 394 cones2.SetParName("ConesTransmission"); 395 392 396 MSimPointingPos pointing; 393 397 … … 666 670 } 667 671 tasks.AddToList(&cones); 672 tasks.AddToList(&cones2); 668 673 //if (header.IsPointRun()) 669 674 // tasks.AddToList(&fillP); … … 834 839 835 840 TCanvas &d = fDisplay->AddTab("Info2"); 836 d.Divide(2, 2);841 d.Divide(2,3); 837 842 838 843 d.cd(1); … … 844 849 splinepde.DrawClone()->SetBit(kCanDelete); 845 850 851 d.cd(3); 852 gPad->SetBorderMode(0); 853 gPad->SetFrameBorderMode(0); 854 gPad->SetGridx(); 855 gPad->SetGridy(); 856 gROOT->SetSelectedPad(0); 857 splinecones2.DrawClone()->SetBit(kCanDelete); 858 859 d.cd(5); 860 gPad->SetBorderMode(0); 861 gPad->SetFrameBorderMode(0); 862 gPad->SetGridx(); 863 gPad->SetGridy(); 864 gROOT->SetSelectedPad(0); 865 splinemirror.DrawClone()->SetBit(kCanDelete); 866 846 867 d.cd(2); 847 868 gPad->SetBorderMode(0); … … 850 871 gPad->SetGridy(); 851 872 gROOT->SetSelectedPad(0); 852 splinemirror.DrawClone()->SetBit(kCanDelete);853 854 d.cd(3);855 gPad->SetBorderMode(0);856 gPad->SetFrameBorderMode(0);857 gPad->SetGridx();858 gPad->SetGridy();859 gROOT->SetSelectedPad(0);860 873 splinecones.DrawClone()->SetBit(kCanDelete); 861 862 d.cd(4); 874 //splinecones2.DrawClone("same")->SetBit(kCanDelete); 875 876 d.cd(6); 863 877 gPad->SetBorderMode(0); 864 878 gPad->SetFrameBorderMode(0); … … 869 883 //all->SetTitle("Combined acceptance"); 870 884 all->SetBit(kCanDelete); 871 all->Multiply(*splinemirror.GetSpline()); 885 if (splinemirror.GetSpline()) 886 all->Multiply(*splinemirror.GetSpline()); 887 if (splinecones2.GetSpline()) 888 all->Multiply(*splinecones2.GetSpline()); 872 889 } 873 890 -
trunk/MagicSoft/Mars/msim/MSimMMCS.h
r9359 r9525 17 17 { 18 18 private: 19 MMcEvtBasic *fMcEvtBasic; 20 MMcEvt *fMcEvt; 19 MMcEvtBasic *fMcEvtBasic; //! 20 MMcEvt *fMcEvt; //! 21 21 MPointingPos *fPointingTel; //! telescope poiting position in local (telescope) coordinate system 22 MCorsikaEvtHeader *fEvtHeader; 23 MCorsikaRunHeader *fRunHeader; 24 MMcRunHeader *fMcRunHeader; 22 MCorsikaEvtHeader *fEvtHeader; //! 23 MCorsikaRunHeader *fRunHeader; //! 24 MMcRunHeader *fMcRunHeader; //! 25 25 26 26 // MTask -
trunk/MagicSoft/Mars/msimcamera/MSimCamera.h
r9425 r9525 28 28 MMcEvt *fMcEvt; //! For information stored in MMcEvt 29 29 30 const MSpline3 *fSpline; // Pu sle Shape30 const MSpline3 *fSpline; // Pulse Shape 31 31 32 32 Bool_t fBaselineGain; // Should the gain be applied to baseline and electronic noise? -
trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc
r9521 r9525 52 52 // R = R0 * Ai * f 53 53 // 54 // R0 is the night sky background rate as given in Eckart's paper. Ai the 55 // area of the cones acceptance window, f is given as: 54 // R0 is the night sky background rate as given in Eckart's paper (divided 55 // by the wavelength window). Ai the area of the cones acceptance window, 56 // f is given as: 56 57 // 57 58 // f = nm * Min(Ar, sr*d^2) … … 68 69 // sr is the effective solid angle corresponding to the integral of the 69 70 // cones angular acceptance 71 // 72 // Alternatively, the night-sky background rate can be calculated from 73 // a spectrum as given in Fig. 1 (but versus Nanometers) in 74 // 75 // Chris R. Benn & Sara L. Ellison La Palma Night-Sky Brightness 76 // 77 // After proper conversion of the units, the rate of the pixel 0 78 // is then calculated by 79 // 80 // rate = f * nsb 81 // 82 // With nsb 83 // 84 // nsb = Integral(nsb spectrum * combines efficiencies) 85 // 86 // and f can be either 87 // 88 // Eff. angular acceptance Cones (e.g. 20deg) * Cone-Area (mm^2) 89 // f = sr * A0 90 // 91 // or 92 // 93 // Mirror-Area * Field of view of cones (deg^2) 94 // f = Ar * A0; 70 95 // 71 96 // … … 118 143 MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title) 119 144 : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0), 120 fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam") 145 fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"), 146 fFileNameNSB("resmc/night-sky-la-palma.txt") 121 147 { 122 148 fName = name ? name : "MSimRandomPhotons"; … … 169 195 }*/ 170 196 171 fRunHeader = 0; 172 if (fSimulateWavelength) 173 { 174 fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); 175 if (!fRunHeader) 197 fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); 198 if (fSimulateWavelength && !fRunHeader) 199 { 200 *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; 201 return kFALSE; 202 } 203 204 MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector"); 205 if (!r) 206 { 207 *fLog << err << "Reflector [MReflector] not found... aborting." << endl; 208 return kFALSE; 209 } 210 211 const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline"); 212 const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesTransmission", "MParSpline"); 213 const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity", "MParSpline"); 214 const MParSpline *s4 = (MParSpline*)pList->FindObject("ConesAngularAcceptance", "MParSpline"); 215 216 // Multiply all relevant efficiencies to get the total tarnsmission 217 MParSpline *eff = (MParSpline*)s1->Clone(); 218 eff->Multiply(*s2->GetSpline()); 219 eff->Multiply(*s3->GetSpline()); 220 221 // Effectively transmitted wavelength band 222 const Double_t nm = eff && eff->GetSpline() ? eff->GetSpline()->Integral() : 1; 223 224 // Angular acceptance of the cones 225 const Double_t sr = s4 && s4->GetSpline() ? s4->GetSpline()->IntegralSolidAngle() : 1; 226 227 { 228 const Double_t d2 = fGeom->GetCameraDist()*fGeom->GetCameraDist(); 229 const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad(); 230 const Double_t f1 = TMath::Min(r->GetA()/1e4, sr*d2) * conv*conv; 231 232 // Rate in GHz / mm^2 233 fScale = fFreqNSB * nm * f1; // [GHz/mm^2] efficiency * m^2 *rad^2 *mm^2 234 235 const Double_t freq0 = fScale*(*fGeom)[0].GetA()*1000; 236 237 *fLog << inf << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", freq0) << "MHz" << endl; 238 239 // FIXME: Scale with the number of pixels 240 if (freq0>1000) 176 241 { 177 *fLog << err << " MCorsikaRunHeader not found... aborting." << endl;242 *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl; 178 243 return kFALSE; 179 244 } 180 245 } 181 246 182 MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector"); 183 if (!r) 184 { 185 *fLog << err << "Reflector [MReflector] not found... aborting." << endl; 186 return kFALSE; 187 } 188 189 *fLog << warn << "FIXME: A check is needed that the PDE is 0 at both ends!" << endl; 190 191 const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline"); 192 const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance", "MParSpline"); 193 const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity", "MParSpline"); 194 const MParSpline *s5 = (MParSpline*)pList->FindObject("ConesTransmission", "MParSpline"); 195 196 const Double_t d2 = fGeom->GetCameraDist()*fGeom->GetCameraDist(); 197 // const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1; 198 const Double_t sr = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1; 199 // const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1; 200 const Double_t Ar = r->GetA()/1e4; 201 202 // Conversion factor to convert pixel area to steradians (because it 203 // is a rather small area we can assume it is flat) 204 const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad(); 205 206 // Multiply all relevant efficiencies 207 MParSpline *s4 = (MParSpline*)s1->Clone(); 208 s4->Multiply(*s3->GetSpline()); 209 s4->Multiply(*s5->GetSpline()); 210 211 const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1; 212 213 delete s4; 214 215 // /100 to convert the pixel area from mm^2 to cm^2 216 fScale = nm * TMath::Min(Ar, sr*d2) * conv*conv; 217 218 *fLog << inf; 219 *fLog << "Effective cone acceptance: " << Form("%.2f", sr*d2) << "m^2" << endl; 220 *fLog << "Reflector area: " << Form("%.2f", Ar) << "m^2" << endl; 221 *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl; 222 // *fLog << "Eff. wavelength band (PDE): " << Form("%.1f", pde) << "nm" << endl; 223 // *fLog << "Eff. wavelength band (Mirror): " << Form("%.1f", mir) << "nm" << endl; 224 *fLog << "Eff. wavelength band (MIR+Cone+PDE): " << Form("%.1f", nm) << "nm" << endl; 225 *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl; 226 //*fLog << "Effective angular acceptance: " << sr << "sr" << endl; 227 //*fLog << "Resulting NSB frequency: " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl; 228 *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl; 247 if (fFileNameNSB.IsNull()) 248 { 249 delete eff; 250 return kTRUE; 251 } 229 252 230 253 // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader"); … … 237 260 // ampl = MMcFadcHeader->GetAmplitud() 238 261 // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio) 262 263 // Conversion of the y-axis 264 // ------------------------ 265 // Double_t ff = 1; // myJy / arcsec^2 per nm 266 // ff *= 1e-6; // Jy / arcsec^2 per nm 267 // ff *= 3600*3600; // Jy / deg^2 268 // ff *= 1./TMath::DegToRad()/TMath::DegToRad(); // Jy/sr = 1e-26J/s/m^2/Hz/sr 269 // ff *= 1e-26; // J/s/m^2/Hz/sr per nm 270 271 const Double_t arcsec2rad = TMath::DegToRad()/3600.; 272 const Double_t f = 1e-32 / (arcsec2rad*arcsec2rad); 273 274 // Read night sky background flux from file 275 MParSpline flux; 276 if (!flux.ReadFile(fFileNameNSB)) 277 return kFALSE; 278 279 const Int_t min = TMath::FloorNint(flux.GetXmin()); 280 const Int_t max = TMath::CeilNint( flux.GetXmax()); 281 282 if (fRunHeader) 283 { 284 if (min>fRunHeader->GetWavelengthMin()) 285 { 286 *fLog << warn << "WARNING - Minimum wavelength of night sky background flux ("; 287 *fLog << min << "nm) from " << fFileNameNSB; 288 *fLog << " exceeds minimum wavelength simulated "; 289 *fLog << fRunHeader->GetWavelengthMin() << "nm." << endl; 290 } 291 if (max<fRunHeader->GetWavelengthMax()) 292 { 293 *fLog << warn << "WARNING - Maximum wavelength of night sky background flux ("; 294 *fLog << max << "nm) from " << fFileNameNSB; 295 *fLog << " undershoots maximum wavelength simulated "; 296 *fLog << fRunHeader->GetWavelengthMax() << "nm." << endl; 297 } 298 } 299 300 MParSpline nsb; 301 302 // Normalization to our units, 303 // conversion from energy flux to photon flux 304 nsb.SetFunction(Form("%.12e/(x*TMath::H())", f), max-min, min, max); 305 306 // multiply night sky background flux with normalization 307 nsb.Multiply(*flux.GetSpline()); 308 309 // Multiply with the total transmission 310 nsb.Multiply(*eff->GetSpline()); 311 312 // Check if the photon flux is zero at both ends 313 if (nsb.GetSpline()->Eval(min)>1e-5) 314 { 315 *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << min << "nm is not zero... abort." << endl; 316 return kFALSE; 317 } 318 if (nsb.GetSpline()->Eval(max)>1e-5) 319 { 320 *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << max << "nm is not zero... abort." << endl; 321 return kFALSE; 322 } 323 324 if (fRunHeader) 325 { 326 if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMin())>1e-5) 327 *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMin() << "nm is not zero... abort." << endl; 328 if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMax())>1e-5) 329 *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMax() << "nm is not zero... abort." << endl; 330 } 331 332 // Conversion from m to radians 333 const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad()*1e3; 334 335 // Angular acceptance of the cones 336 //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr 337 // Absolute reflector area 338 const Double_t Ar = r->GetA()/1e4; // m^2 339 // Size of the cone's entrance window 340 const Double_t A0 = (*fGeom)[0].GetA()*1e-6; // m^2 341 342 // Rate * m^2 * Solid Angle 343 // ------------------------- 344 345 // Angular acceptance Cones (e.g. 20deg) * Cone-Area 346 const Double_t f1 = A0 * sr; // m^2 sr 347 348 // Mirror-Area * Field of view of cones (e.g. 0.1deg) 349 const Double_t f2 = Ar * A0*conv*conv; // m^2 sr 350 351 // FIXME: Calculate the reflectivity of the bottom by replacing 352 // MirrorReflectivity by bottom reflectivity and reflect 353 // and use it to reflect the difference between f1 and f2 354 // if any. 355 356 // Total NSB rate in MHz per m^2 and sr 357 const Double_t rate = nsb.GetSpline()->Integral() * 1e-6; 358 359 *fLog << inf; 360 361 // Resulting rates as if Razmick's constant had been used 362 // *fLog << 1.75e6/(600-300) * f1 * eff->GetSpline()->Integral() << " MHz" << endl; 363 // *fLog << 1.75e6/(600-300) * f2 * eff->GetSpline()->Integral() << " MHz" << endl; 364 365 *fLog << "Conversion factor Fnu: " << f << endl; 366 *fLog << "Total reflective area: " << Form("%.2f", Ar) << " m" << UTF8::kSquare << endl; 367 *fLog << "Acceptance area of cone 0: " << Form("%.2f", A0*1e6) << " mm" << UTF8::kSquare << " = "; 368 *fLog << A0*conv*conv << " sr" << endl; 369 *fLog << "Cones angular acceptance: " << sr << " sr" << endl; 370 *fLog << "ConeArea*MirrorAngle (f1): " << f1 << " m^2 sr" << endl; 371 *fLog << "MirrorArea*ConeAngle (f2): " << f2 << " m^2 sr" << endl; 372 *fLog << "Effective. transmission: " << Form("%.1f", nm) << " nm" << endl; 373 *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f1): " << Form("%.2f", rate * f1) << " MHz" << endl; 374 *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f2): " << Form("%.2f", rate * f2) << " MHz" << endl; 375 *fLog << "Using f1." << endl; 376 377 // Scale the rate per mm^2 and to GHz 378 fScale = rate * f1 / (*fGeom)[0].GetA() / 1000; 379 380 // FIXME: Scale with the number of pixels 381 if (rate*f1>1000) 382 { 383 *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl; 384 return kFALSE; 385 } 386 387 delete eff; 239 388 240 389 return kTRUE; … … 274 423 { 275 424 // Scale the rate with the pixel size. 276 const Double_t rate = fFreqFixed +fFreqNSB*(*fGeom)[idx].GetA()*fScale;425 const Double_t rate = fFreqFixed + fScale*(*fGeom)[idx].GetA(); 277 426 278 427 (*fRates)[idx].SetPedestal(rate); … … 305 454 // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant) FIXME: Reset? 306 455 307 if (f RunHeader)456 if (fSimulateWavelength) 308 457 { 309 458 const Float_t wmin = fRunHeader->GetWavelengthMin(); … … 353 502 } 354 503 504 if (IsEnvDefined(env, prefix, "FileNameNSB", print)) 505 { 506 rc = kTRUE; 507 fFileNameNSB = GetEnvValue(env, prefix, "FileNameNSB", fFileNameNSB); 508 } 509 510 if (IsEnvDefined(env, prefix, "SimulateCherenkovSpectrum", print)) 511 { 512 rc = kTRUE; 513 fSimulateWavelength = GetEnvValue(env, prefix, "SimulateCherenkovSpectrum", fSimulateWavelength); 514 } 515 355 516 return rc; 356 517 } -
trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h
r9425 r9525 33 33 34 34 TString fNameGeomCam; 35 TString fFileNameNSB; 35 36 36 37 // MTask
Note:
See TracChangeset
for help on using the changeset viewer.