//************************************************************************************** // 1event demo script of ROMEO // // Usage: // root -l 'RomeoDemo.C(.....)' // // Arguments: // // maxpe ... input NPE per shot. // beamwidth ... time width of a beam shot. By default // time property of a beam is uniform within the beamwidth. // weightHits ... if True or 1, all PEs within a shot are merged into 1 MCHit with weight. // isAmasimMode ... if True or 1 romeo treats input beam as PhotoElectrons // instead of Photons. IceSim default is True. // isQuickMode ... if True or 1 romeo skips photocathode simulation. // IceSim default is True. // doWaveformSaturation ... if true or 1, it calculates saturation at waveform. // if false or 0, romeo activates RMOPMTCurrentBasedSaturation // which uses average current within outputtimeresol // as an input of saturation curve. (semi-current saturation mode) // outputtimeresol ... any hit within this timewindow will be merged into one // charge object. It's also used to calculate average current // at semi-current saturation mode. // Negative falue is replaced to 2ns. If you set 0, binning is skipped and // RMOPMTCurrentBasedSaturation is not activated. // nbeamshots ... number of beam shots. beam bunch is sequently shot every "beamwidth". // // Beaming examples : // // A) weightHits = kFALSE, maxpe=4, beamwidth=10ns, nbeamshots=4: // // |||||||||||||||| weight 1 / hit // ---- // 10ns // |______________| 4 beam shots // // B) weightHits = kTRUE, maxpe=4, beamwidth=10ns, nbeamshots=4: // // | | | | weight 4 / hit // ---- // 10ns // |______________| 4 beam shots // //************************************************************************************** #include #include #include //TFile foutp; int RomeoDemo(Int_t maxpe = 1000, // input NPE per shot Double_t beamwidth = 1e-8, // 10ns square beam Bool_t weightHits = kFALSE, // will you make weighted hit? Bool_t isAmasimMode=kFALSE, Bool_t isQuickMode=kFALSE, Bool_t doWaveformSaturation=kFALSE, // do saturation at waveform Double_t outputtimeresol = 0, // set positive number if you need charge-bin Int_t nbeamshots = 1) // number of beam shot { gROOT->Reset(); gSystem->Load("libPhysics.so"); gSystem->Load("libHist.so"); gSystem->Load("libGraf3d.so"); gSystem->Load("../lib/libRMONumericalRecipes.so"); gSystem->Load("../lib/libRMOKernel.so"); gSystem->Load("../lib/libRMOOMAnalog.so"); gSystem->Load("../lib/libRMOGlassGel.so"); gSystem->Load("../lib/libRMOPMT.so"); //foutp = new TFile(outname,"RECREATE"); gStyle->SetOptStat(1111); gStyle->SetPalette(1); gStyle->SetCanvasColor(33); gStyle->SetFrameFillColor(18); RMOParameters::SetVerboseLevel(0); //====================================================================== // 1) define data dir // you must chose one data directory which includes your calibration data RMOParameters::SetDataDir("../../resources/data/"); //====================================================================== // 2) create modules // create your modules. // * glass+gel RMOOMAnalog *omanalog = new RMOOMAnalog("StandardDOMIce3"); RMOPhotonGenerator *phgen = new RMOPhotonGenerator("PhotonGeneratorIce3"); RMOWaveLengthSimulator *wlsim = new RMOWaveLengthSimulator("WaveLengthSimIce3"); RMOIsotropicAngleSimulator *angsim = new RMOIsotropicAngleSimulator("IsotropicAngleSimIce3"); RMOQuickPhotonPropagator *phprop = new RMOQuickPhotonPropagator("QuickPropagatorIce3"); // * photo-cathode RMOPMTPhotoCathode *cathode = new RMOPMTPhotoCathode("PhotoCathodeSimIce3"); // old constructor // RMOPMTPhotoCathode *cathode = new RMOPMTPhotoCathode("TA1062.table", // "qe_representative.data", // "qe_representative.data", // "PhotoCathodeIce3TA1062",""); // this is alternative PE generator. // If you use the generator, you should skip all modules listed above. RMOPMTQuickPEGenerator *pegenerator = new RMOPMTQuickPEGenerator("QuickPEGeneratorIce3"); // * charge response and gain RMOPMTChargeResponse *qresponse = new RMOPMTChargeResponse("ChargeResIce3"); // there are two saturation sub-modules. //RMOPMTChargeBasedSaturation *saturation = new RMOPMTChargeBasedSaturation("ChargeBasedSaturationIce3"); RMOPMTCurrentBasedSaturation *saturation = new RMOPMTCurrentBasedSaturation("CurrentBasedSaturationIce3"); // * waveform RMOPMTWaveForm *waveform = new RMOPMTWaveForm("WaveFormIce3"); //====================================================================== // 3) set module parameters. // all configureable parameters should be defined here //-------------------------- // set verbose... // 0 : no output // 1 : default // 2 : function called // 3 : precise // 4 : debug // 5 : trace omanalog->SetVerboseLevel(2); phgen->SetVerboseLevel(2); wlsim->SetVerboseLevel(2); angsim->SetVerboseLevel(2); phprop->SetVerboseLevel(2); cathode->SetVerboseLevel(2); pegenerator->SetVerboseLevel(2); qresponse->SetVerboseLevel(2); saturation->SetVerboseLevel(2); waveform->SetVerboseLevel(2); //-------------------------- // for glass+gel // choose proper wavelength acceptance with current mode. // Note that you must use "withoutQE" wavelength-acceptance // if you kill the AmasimMode. if (isAmasimMode) { RMOWaveLengthAcceptance::SetInputFileName("wavelength_acceptance_withQE_13inch_0deg.data"); RMOAngleAcceptance::SetInputFileName("angle_acceptance_withQE_13inch_420nm.data"); } else { RMOWaveLengthAcceptance::SetInputFileName("wavelength_acceptance_withoutQE_13inch_0deg.data"); //RMOWaveLengthAcceptance::SetInputFileName("wavelength_acceptance_withoutQE_8mm_0deg.data"); RMOAngleAcceptance::SetInputFileName("angle_acceptance_withoutQE_13inch_420nm.data"); //RMOAngleAcceptance::SetInputFileName("angle_acceptance_withoutQE_8mm_420nm.data"); } RMOQuickPhotonPropagator::SetCalibFileNames("QuickPhotonPropagatorMap_13inch_420nm.data"); //RMOQuickPhotonPropagator::SetCalibFileNames("QuickPhotonPropagatorMap_8mm_385nm.data"); //---------------------------- // for photo-cathode // if you want to change QECE manualy set here. default is TA0003 //cathode->Set2DCEFactorAt337nm(0.228897); // CE of TA1062 based on the Chiba calibration //---------------------------- // for charge-response // if you want to change spe parameters manualy set here. default is TA0003 //qresponse->SetSQ0_ov_Q0(3.25837e-01); // Charge Response of TA1062 //qresponse->SetQt_ov_Q0(1.609385e-01); //qresponse->SetPe(2.228744e-01); //qresponse->SetGain(5.8e4); // if icetop low gain // Output time resolution: // its not PMT time resolution but a kind of binning window in order to // generate one charge datum from several PEs within same time window // default is 2ns // If you set 0, it doesn't any binning and // you must switch off RMOCurrentBasedSaturation module. if (outputtimeresol < 0) { cerr << "you put invalid (negative) output time resolution. " << "2ns will be set instead." << endl; outputtimeresol = 2.e-9; } RMOParameters::SetOutputTimeResolution(outputtimeresol); // PMT time resolution: // Default is 4ns. However, recent version of hit-constructor (or hit-maker) // convolutes this value: we need to set 0 if you use romeo with icesim. // For calibration experiment, just leave it as default. RMOParameters::SetPMTTimeResolution(0); //====================================================================== // 4) connect all modules // if you run the script with argument "true", it skips all simulation // of glass+gel and Photo-Cathode simulation. if (isQuickMode) { // add pe generator only omanalog->Add(pegenerator); } else { // simulate photo-cathode omanalog->Add(phgen); omanalog->Add(wlsim); omanalog->Add(angsim); omanalog->Add(phprop); omanalog->Add(cathode); } omanalog->Add(qresponse); if (!doWaveformSaturation) { if (outputtimeresol > 0) { // add saturation module only when you have positive // outputtimeresolution qresponse->Add(saturation); } // switch off waveform saturation (default) waveform->DoSaturationAtWaveform(kFALSE); } else { // set saturation flag on waveform->DoSaturationAtWaveform(kTRUE); } omanalog->Add(waveform); //====================================================================== // 5) Set amasim mode or kill it. //------------------------------------------------------------------ // Set mode flag AFTER added all modules. Default mode is AmasimMode. // AmasimMode = kTRUE : // All incoming photons must reach to the photocathode. // The photon absorption at glass+gel and incoming angle efficiency are // already took into account the number of arrival photons. // ROMEO assigns wavelength and incoming angle for all photons // following analytical model of wavelength spectrum and angle acceptance // without killing any photon. // Other : // All incoming photons are regarded photons at DOM surface. // ROMEO does fast simulation of photon absorption and angle acception // inside the glass_gel modules then some photons are killed by following // wavelength acceptance and angle acceptance. // RMOAttModeFlag::SetAmasimModeAll(isAmasimMode); //====================================================================== // 6) Simulation start. // Each line contains "ExecuteTask" calls specified function for // all module added to omanalog //------------------------------------------------------------------ //---------------------------- // Configure // this function should be called beforhand the Process call // and should be called only once per one run. omanalog->ExecuteTask("Configure"); //---------------------------- // Process // this is the main part of simulation. it can be inside the event loop. // For example of multiple event, see RomeoManyRun.C // create input buffer Int_t ipe = 0; RMOMCPhotonHitBuf hitarray("MCPhotonTimeHit"); hitarray.SetOwner(); // set owner to allow to delete stored components Double_t t0 = 100.e-9; for (Int_t itbin=0; itbinSetTime(time); hit->SetDirection(dir); hit->SetDistanceFromSource(dist); hit->SetSourceType(10); // 10 is SPE pulse hit->SetWeight(maxpe); hitarray.Add(hit); } else { for (ipe=0; ipeExp(tau); //Double_t time = 100.e-9 + gRandom->Uniform(0, 10e-9); //Double_t time = 30.e-9; Double_t time = t0 + itbin*beamwidth + gRandom->Uniform(0, beamwidth); Double_t dist = 50; // 50m //Double_t dist = 100; // 100m TVector3 dir(0, 0, -1); //TVector3 dir(-1, -1, 0); //TVector3 dir(-1, 0, -1); //TVector3 dir(1, 1, -1); //TVector3 dir(1, 1, 0); RMOMCPhotonHit *hit = new RMOMCPhotonHit(); hit->SetTime(time); hit->SetDirection(dir); hit->SetDistanceFromSource(dist); hit->SetSourceType(10); // 10 is SPE pulse hitarray.Add(hit); } } } // set input buffer to first module omanalog->SetBuffer(&hitarray); //**************** // call Process. //**************** omanalog->ExecuteTask("Process"); //**************** // simulation has done! //**************** // output information Int_t ngenpe = maxpe; if (!isQuickMode) { RMOPMTPEBuf &pebuf = (RMOPMTPEBuf&)cathode->GetBuffer(); ngenpe = pebuf.GetEntries(); } RMOPMTWaveFormBuf &wfbuf = (RMOPMTWaveFormBuf&)waveform->GetBuffer(); Double_t peakV = wfbuf.GetMinimum(); Double_t charge = wfbuf.Integral(wfbuf.GetFirstTime(), wfbuf.GetLastTime()) / 43.; // 43 ohm, IceCube default cerr << "***********************************" << endl; cerr << " info of output pulse " << endl; cerr << " ----------------------------------" << endl; cerr << " Input ---------- " << endl; cerr << " * Input ph / pe : " << maxpe << endl; cerr << " From PhotoCathode ---------- " << endl; cerr << " * Npe (survived) : " << ngenpe << endl; cerr << " From ChargeResponse ---------- " << endl; cerr << " * Total charge[C](before saturation) : " << qresponse->GetTotalCharge() << endl; cerr << " * Estimated Npe (before saturation) : " << qresponse->GetTotalNpe() << endl; cerr << " * Total charge[C](after saturation) : " << qresponse->GetSaturatedTotalCharge() << endl; cerr << " * Estimated Npe (after saturation) : " << qresponse->GetSaturatedTotalNpe() << endl; cerr << " From Waveform ---------- " << endl; cerr << " * peak amplitude[V] : " << peakV << endl; cerr << " * Total charge[C] : " << (- charge) << endl; cerr << " * tmin / tmax [s]: " << wfbuf.GetXmin() << " / " << wfbuf.GetXmax() << endl; cerr << " * Estimated Npe : " << (- charge) / (1e7 *1.60217733e-19) << endl; cerr << "***********************************" << endl; wfbuf.DumpInfo(); //---------------------------- // Finish // this function should be called after the Process call // and should be called only once per one run. omanalog->ExecuteTask("Finish"); //---------------------------- // Draw... // if you want to see histograms, call Draw. omanalog->ExecuteTask("Draw"); }