/* This is an example setup for an analysis of multiplicities and inclusive quantities in hadron decays. It has to include Main_FullDecay.C, and a FullDecay has to be started with the following commands (in the directory containing the provided Run.dat): - ./FullDecay DECAYER= EVENTS=100000 ANALYSIS=1 */ #include "Main_FullDecay.C" #include "ATOOLS/Org/Shell_Tools.H" #ifdef USING__ROOT static TFile* rootfile; static map multiplicities; static map energies; static map directenergies; typedef map::iterator FTMapIt; static map int_multiplicities; typedef map::iterator FIMapIt; #endif void InitialiseAnalysis() { #ifdef USING__ROOT std::string adir = "Analysis"; ATOOLS::MakeDir(adir,493); rootfile = new TFile(string(adir+"/Multiplicities_"+ mother_flav.ShellName()+".root").c_str(), "RECREATE"); PRINT_INFO("Writing to: "< flavours; flavours.push_back(Flavour(kf_pi_plus)); flavours.push_back(Flavour(kf_pi_plus).Bar()); flavours.push_back(Flavour(kf_K_plus)); flavours.push_back(Flavour(kf_K_plus).Bar()); flavours.push_back(Flavour(kf_K_L)); flavours.push_back(Flavour(kf_e)); flavours.push_back(Flavour(kf_e).Bar()); flavours.push_back(Flavour(kf_nue)); flavours.push_back(Flavour(kf_nue).Bar()); flavours.push_back(Flavour(kf_numu)); flavours.push_back(Flavour(kf_numu).Bar()); flavours.push_back(Flavour(kf_mu)); flavours.push_back(Flavour(kf_mu).Bar()); flavours.push_back(Flavour(kf_tau)); flavours.push_back(Flavour(kf_tau).Bar()); flavours.push_back(Flavour(kf_nutau)); flavours.push_back(Flavour(kf_nutau).Bar()); flavours.push_back(Flavour(kf_n)); flavours.push_back(Flavour(kf_n).Bar()); flavours.push_back(Flavour(kf_p_plus)); flavours.push_back(Flavour(kf_p_plus).Bar()); flavours.push_back(Flavour(kf_photon)); for(list::iterator it=flavours.begin(); it!=flavours.end(); it++) { multiplicities[*it] = makeTH1D(it->ShellName()+"_multi", "", 9, -0.5, 8.5, it->RootName()+" multiplicity","Events"); energies[*it] = makeTH1D(it->ShellName()+"_energy", "", 25, 0.0, mother_flav.HadMass()/2.0, "E_{"+it->RootName()+"}", "#frac{1}{#Gamma} #frac{d#Gamma}{dE}"); directenergies[*it] = makeTH1D(it->ShellName()+"_directenergy", "", 25, 0.0, mother_flav.HadMass()/2.0, "E_{"+it->RootName()+"}", "#frac{1}{#Gamma} #frac{d#Gamma}{dE}"); int_multiplicities[*it] = 0; } #endif } void AnalyseEvent(Blob_List* blobs) { #ifdef USING__ROOT int outgoing=1; Particle_List outparts = blobs->ExtractParticles(part_status::active, outgoing); for(Particle_List::iterator part=outparts.begin(); part!=outparts.end(); part++) { if(multiplicities.find((*part)->Flav())!=multiplicities.end()) { int_multiplicities[(*part)->Flav()]++; energies[(*part)->Flav()]->Fill((*part)->Momentum()[0]); if((*part)->ProductionBlob()->InParticle(0)->Flav()==mother_flav) { directenergies[(*part)->Flav()]->Fill((*part)->Momentum()[0]); } } else { PRINT_INFO("stable but unanalysed outparticle encountered: " <<**part); } } for(FIMapIt it=int_multiplicities.begin(); it!=int_multiplicities.end(); it++) { Flavour flav = it->first; multiplicities[flav]->Fill(it->second); it->second = 0; } #endif } void FinishAnalysis() { #ifdef USING__ROOT for(FTMapIt it=multiplicities.begin(); it!=multiplicities.end(); it++) { it->second->Write(); } for(FTMapIt it=energies.begin(); it!=energies.end(); it++) { it->second->Write(); } for(FTMapIt it=directenergies.begin(); it!=directenergies.end(); it++) { it->second->Write(); } #endif }