/* This is an example setup for an analysis of photon multiplicities and radiative energy loss as well as an angular radiation pattern analysis. It has to include Main_FullDecay.C, and a FullDecay has to be started with the following specialities: - ./PhotonAnalysis DECAYER= EVENTS=100000 ANALYSIS=1 */ #include "Main_FullDecay.C" #include "ATOOLS/Org/Shell_Tools.H" #ifdef USING__ROOT static TFile* rootfile; static TH1D * photonmultiplicity; static TH1D * decayframeenergy; static TH1D * multipoleframeangles; #endif void InitialiseAnalysis() { #ifdef USING__ROOT std::string adir = ""; ATOOLS::MakeDir("PhotonAnalysisDirectory/"+mother_flav.ShellName()+"_decays",true,493); rootfile = new TFile(string("PhotonAnalysisDirectory/"+ mother_flav.ShellName()+"_decays/"+ mother_flav.ShellName()+"__"+ "DAUGHTERS"+".root").c_str(), "RECREATE"); photonmultiplicity = makeTH1D("photon_multiplicity","", 10, -0.5, 10.5, Flavour(kf_photon).RootName()+" multiplicity","Events"); decayframeenergy = makeTH1D("decayframeenergy","", 1000, 0., mother_flav.HadMass(), "total energy radiated in decay frame","Events"); multipoleframeangles = makeTH1D("multipoleframeangles","", 1000, 0., M_PI, "angular radiation pattern","Events"); #endif } void AnalyseEvent(Blob_List* blobs) { #ifdef USING__ROOT // int outgoing = 1; // int incoming = -1; // Particle_List outparts = blobs->ExtractParticles(part_status::active, outgoing); /////////////////////////////////////////////////////////////////////////////////// // analyse primary decay blob, ignore subsequent decays // /////////////////////////////////////////////////////////////////////////////////// Blob * primarydecayblob = blobs->FindFirst(btp::Hadron_Decay); // msg_Out()<<"primary decay blob:"<NOutP(); i++) { if ((primarydecayblob->OutParticle(i)->Flav().IsPhoton() == true) && (primarydecayblob->OutParticle(i)->Info() == 'S')) { photmult++; photener = photener + primarydecayblob->OutParticle(i)->Momentum()[0]; } } photonmultiplicity->Fill(photmult); if (photener != 0.) decayframeenergy->Fill(photener); // multipole rest frame angles Vec4D multipolesum = Vec4D(0.,0.,0.,0.); Vec4D axis = Vec4D(0.,0.,0.,1.); std::list multipole; std::list newphot; for (int i=0; iNOutP(); i++) { if (primarydecayblob->OutParticle(i)->Flav().Charge() != 0.) { multipolesum = multipolesum + primarydecayblob->OutParticle(i)->Momentum(); multipole.push_back(primarydecayblob->OutParticle(i)->Momentum()); } } if (primarydecayblob->InParticle(0)->Flav().Charge() != 0) { multipolesum = multipolesum + primarydecayblob->InParticle(0)->Momentum(); multipole.push_front(primarydecayblob->InParticle(0)->Momentum()); } Poincare boost(multipolesum); Poincare rotate; // charged initial state: rotate such that initial state at theta = 0 std::list::iterator heaviest = multipole.begin(); if (mother_flav.Charge() != 0.) { Vec4D inmom = *multipole.begin(); boost.Boost(inmom); rotate = Poincare(inmom,axis); } // neutral initial state: rotate such that heaviest charged final state at theta = 0 else { for (std::list::iterator iter=multipole.begin(); iter!=multipole.end(); iter++) { if (abs((iter->Abs2() - heaviest->Abs2())/(iter->Abs2() + heaviest->Abs2())) > 1E-6) { heaviest = iter; } } boost.Boost(*heaviest); rotate = Poincare(*heaviest,axis); } for (int i=0; iNOutP(); i++) { if (primarydecayblob->OutParticle(i)->Flav().IsPhoton() == true) { Vec4D mom = primarydecayblob->OutParticle(i)->Momentum(); boost.Boost(mom); rotate.Rotate(mom); double theta = acos((Vec3D(mom)*Vec3D(axis))/(Vec3D(mom).Abs()*Vec3D(axis).Abs())); multipoleframeangles->Fill(theta); } } /////////////////////////////////////////////////////////////////////////////////// // inclusive analysis of whole decay chain // /////////////////////////////////////////////////////////////////////////////////// // to be done .. #endif } void FinishAnalysis() { #ifdef USING__ROOT photonmultiplicity->Write(); decayframeenergy->Write(); multipoleframeangles->Write(); #endif }