#include "MdEventPlotter.h" #include #include #include #include #include #include #include #include #include #include // #include #include #include #include #include #include // #include // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // #include #include #include #include // #include #include // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using std::cout; using std::cerr; using std::endl; using namespace utl; using namespace fwk; namespace MdEventPlotterAG { MdEventPlotter::MdEventPlotter() : fRunNumber(0), fEventNumber(0) { } MdEventPlotter::~MdEventPlotter() { } VModule::ResultFlag MdEventPlotter::Init() { fFile1 = new ofstream("MdLDFs.dat"); *fFile1 << "NEvent, Distance [m], Muon Density [1/m^2], Error [1/m^2]" << endl; // Configuration reading. utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdEventPlotter"); //Output directory LoadConfig(config, "outdir", fOutDir, "./figs_EVENTPLOTTER" ); std::stringstream cmd; cmd << "test -d " << fOutDir << " || mkdir -p " << fOutDir; system( cmd.str().c_str() ); LoadConfig(config, "run", fRun, "run2" ); LoadConfig(config, "plotExtension", fPlotExtension1, ".png" ); LoadConfig(config, "plotExtension2", fPlotExtension2, ".eps" ); LoadConfig(config, "plotExtension3", fPlotExtension3, ".root" ); InitializeStyle(); fStyle->cd(); return eSuccess; } VModule::ResultFlag MdEventPlotter::Run(evt::Event& event) { INFO( "\n\n**************\nSTARTING EVENT ANALYSIS\n**************\n"); //Plot event pattern on ground PlotEvent( event ); //Plot the reconstructed MLDF PlotMLDF( event ); fEventNumber++; return eSuccess; } VModule::ResultFlag MdEventPlotter::PlotMLDF(evt::Event& event) { INFO( "\n\n-----\nPLOTTING LDF\n"); fStyle->cd(); const mdet::MDetector& mdDe = det::Detector::GetInstance().GetMDetector(); typedef mdet::MDetector::CounterConstIterator Itc; const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem(); if (! (event.HasSimShower() || event.HasRecShower()) ) { ERROR("\n---\nCurrent event does not have a simulated or reconstructed shower.\n---\n"); return eFailure; } const evt::ShowerSimData& simShower = event.GetSimShower(); // const Point& corePos = simShower.GetPosition(); // const TimeStamp& coreTime = simShower.GetTimeStamp(); const double showerE = simShower.GetEnergy(); const double showerZenith = simShower.GetZenith(); const double showerAzimuth = simShower.GetAzimuth(); const int showerNumber = simShower.GetShowerNumber(); fRunNumber = showerNumber; std::string showerPrimary, showerModel; if ( simShower.GetPrimaryParticle() == Particle::eProton ) showerPrimary = "pr"; if ( simShower.GetPrimaryParticle() == Particle::eIron ) showerPrimary = "fe"; //showerModel = "qgsjetII"; std::ostringstream os; os.str(""); os << " EVENT NUMBER -> " << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << fEventNumber << endl; os << " SHOWER PRIMARY -> " << showerPrimary << endl; os << " SHOWER ENERGY -> " << std::fixed << std::setprecision(2) << log10(showerE/eV) << endl; os << " SHOWER ZENITH ANGLE -> " << std::setfill('0') << std::setw(2) << std::fixed << std::setprecision(1) << showerZenith/utl::deg << endl; os << " SHOWER AZIMUTH ANGLE -> " << std::setfill('0') << std::setw(2) << std::fixed << std::setprecision(1) << showerAzimuth/utl::deg << endl; os << " SHOWER MODEL -> " << showerModel << endl; os << " SHOWER NUMBER -> " << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << showerNumber << endl; cout << os.str() << endl; //TODO:MLDF from sim shower at a given depth //// if( simShower.HasGroundParticles () ) //// { //// const utl::ShowerParticleList & plist = simShower.GetGroundParticles (); //// for ( ShowerParticleIterator pIt = plist.ParticlesBegin(), pEnd = plist.ParticlesEnd(); pIt != pEnd; ++pIt ){ //// const Particle & part = (*pIt); //// if( part.GetType() != Particle::eMuon || part.GetType() != Particle::eAntiMuon ) continue; //// } //// } // nothing to do if there is no RecShower if (!(event.GetRecShower().HasSRecShower()) ) { ERROR("Current event does not have a SD recontructed shower."); return eFailure; } const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower() ; const utl::Point& rCore = sRecShower.GetCorePosition(); const double rEnergy = sRecShower.GetEnergy(); const utl::Vector& rAxis = sRecShower.GetAxis(); const utl::UTMPoint utmcore( rCore, utl::ReferenceEllipsoid::GetWGS84()); cout << "Rec Core Position= (" << std::fixed << std::setprecision(2) << utmcore.GetEasting()/utl::km << "," << std::fixed << std::setprecision(2) << utmcore.GetNorthing()/utl::km << ") [km]\n" << "Rec Energy = " << std::fixed << std::setprecision(3) << rEnergy/utl::EeV << " [EeV]\n" << "Rec Energy = " << std::fixed << std::setprecision(2) << log10(rEnergy/utl::eV) << " [log(eV)]\n"; if ( !(event.GetRecShower().HasMRecShower()) ) { ERROR("\n---\nCurrent event does not have a MD reconstructed shower\n---\n"); return eFailure; } const evt::ShowerMRecData& mrecShower = event.GetRecShower().GetMRecShower(); // const Point& baryCenter = mrecShower.GetBarycenter(); TGraphErrors grMuons; TGraphAsymmErrors grSilent; TObjArray arSilent; TObjArray liSilent; TObjArray arSaturated; TObjArray liSaturated; TGraph grMLDFFit; if (event.HasMEvent()) { mevt::MEvent& me = event.GetMEvent(); unsigned int slim = me.GetSilentLimit(); TLegend leg; leg.SetNColumns( 1 ); int idx = 0; // int idxSi = 0; for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic ) { if ( ic->GetRecStatus() == mevt::RecStatus::eRejected ) { // JMF 03-May-2013 // JMF 03-May-2013 cout << "Counter " << ic->GetId() << " REJECTED, IT IS NOT IN THE ARRAY "; // JMF 03-May-2013 cout << " --\n"; // JMF 03-May-2013 continue; // JMF 03-May-2013 } // JMF 03-May-2013 if(!ic->HasRecData()) { cout << "Counter " << ic->GetId() << " DOES NOT HAVE REC DATA "; // JMF 03-May-2013 cout << " --\n"; // JMF 03-May-2013 //continue; } //if( fEventNumber == 0 ) fAreaDet = CounterProjectedArea( mdDe.GetCounter( ic->GetId() ), (*ic) , rCore, rAxis, siteCS )/utl::m2; const mdet::Counter& mdCo = mdDe.GetCounter( ic->GetId() ); const utl::Point& r_i = mdCo.GetPosition(); double counterRho = cross( rAxis, r_i - rCore ).GetMag() ; unsigned int nMuons = 0; double MuonDensity = 0; double MuonDensityError = 0; double rerr = 25*utl::m; fAreaDet = CounterProjectedArea( mdCo, (*ic) , rCore, rAxis, siteCS )/utl::m2; if( ic->HasRecData() ) fWindowSize = ic->GetRecData().GetWindowSize(); mevt::RecStatus::Status status = ic->GetRecStatus(); if ( status == mevt::RecStatus::eCandidate ) { cout << "CANDIDATE COUNTER\t"; fCandidates.push_back( counterRho ); nMuons = ic->GetRecData().GetNumberOfCorrectedMuons(); MuonDensity = nMuons/fAreaDet; MuonDensityError = TMath::Sqrt(nMuons)/fAreaDet; // MuonDensity = nMuons/fAreaDet -> Var(MuonDensityError) = Var(nMuons)/fAreaDet^2 grMuons.SetPoint( idx, counterRho/utl::m, MuonDensity ); grMuons.SetPointError( idx, 0, MuonDensityError ); ++idx; // All LDFs of candidate counters are stored for writting them in 'Finish' function LDFCandidates.push_back(boost::make_tuple( fEventNumber, counterRho, MuonDensity , MuonDensityError ) ); } else if ( status == mevt::RecStatus::eSilent ) { cout << "SILENT COUNTER\t"; if( ic->HasRecData() )//silent could not have rec data but are still silent! { nMuons = ic->GetRecData().GetNumberOfCorrectedMuons(); MuonDensity = nMuons/fAreaDet; } fSilents.push_back( counterRho ); TArrow *aSilent = new TArrow ( counterRho/utl::m,// x1 slim/fAreaDet, // y1 counterRho/utl::m,// x2 slim/fAreaDet*(1 - 0.5), // y2 0.02,">" ); arSilent.Add( aSilent ); TLine *lSilent = new TLine ( counterRho/utl::m+rerr/utl::m,// x1 slim/fAreaDet, // y1 counterRho/utl::m-rerr/utl::m,// x2 slim/fAreaDet // y2 ); liSilent.Add( lSilent ); // grSilent.SetPoint( idxSi, counterRho/utl::m, slim ); // grSilent.SetPointError( idxSi, 0, 0, slim*(1 - 0.5), 0 ); // ++idxSi; // All LDFs of silent counters are stored for writting them in 'Finish' function LDFSilents.push_back(boost::make_tuple( fEventNumber, counterRho, slim/fAreaDet , 0 ) ); } else if ( status == mevt::RecStatus::eSaturated ) { cout << "SATURATED COUNTER\t"; nMuons = ic->GetRecData().GetNumberOfCorrectedMuons(); MuonDensity = nMuons/fAreaDet; fSaturated.push_back( counterRho ); TArrow *aSaturated = new TArrow ( counterRho/utl::m,// x1 MuonDensity, // y1 counterRho/utl::m,// x2 MuonDensity*(1 + 0.5), // y2 0.02,">" ); arSaturated.Add( aSaturated ); TLine *lSaturated = new TLine ( counterRho/utl::m+rerr/utl::m,// x1 MuonDensity, // y1 counterRho/utl::m-rerr/utl::m,// x2 MuonDensity // y2 ); liSaturated.Add( lSaturated ); // grSaturated.SetPoint( idxS, counterRho/utl::m, nMuons ); // grSaturated.SetPointError( idxS, 0, 0, nMuons*(1 + 0.5), 0 ); // ++idxS; // All LDFs of saturated counters are stored for writting them in 'Finish' function LDFSaturated.push_back(boost::make_tuple( fEventNumber, counterRho, MuonDensity, 0 ) ); } else if ( status == mevt::RecStatus::eRejected ) { cout << "Rejected COUNTER\t"; if( ic->HasRecData() )//silent could not have rec data but are still silent! { nMuons = ic->GetRecData().GetNumberOfCorrectedMuons(); MuonDensity = nMuons/fAreaDet; } fRejected.push_back( counterRho ); } cout << "Counter Id: " << ic->GetId() << " "; cout << "nMuons: " << nMuons << " "; cout << "Area: "<< fAreaDet << " m^2 "; cout << "MuonDensity: "<< MuonDensity << " 1/m^2 "; cout << "Rho: " << counterRho/utl::m << " m "; cout << "(Rec. status "<< status << ") "; cout << endl; } // MLDF from Fit if (mrecShower.HasMLDF()) { const TabulatedFunctionErrors& ldf = mrecShower.GetMLDF(); const unsigned int nLDFPoints = ldf.GetNPoints(); for (unsigned int i = 0; i < nLDFPoints; ++i) grMLDFFit.SetPoint(i, ldf.GetX(i)/utl::m, ldf.GetY(i) ); } else { INFO("\nMLDF Not FOUND in this SHOWER\n"); return eFailure; } double beta = mrecShower.GetBeta(); double betaErr = mrecShower.GetBetaError(); double nMuRef = mrecShower.GetNMuRef(); double nMuRefErr = mrecShower.GetNMuRefError(); double refDist = mrecShower.GetReferenceDistance(); std::stringstream st; st.str(""); st << "#font[102]{"; st << "#color[2]{#scale[0.6]{#rho_{#mu}(r) = #rho_{#mu}("<< int(refDist/utl::m)<<") #times f(r;#beta)}}"; st << "}"; //TLatex * latexMLDF = new TLatex( refDist/utl::m-100, nMuRef*(10), st.str().c_str() ); TLatex * latexMLDF = new TLatex( 500, 100, st.str().c_str() ); st.str(""); st << "#font[102]{"; st << "#color[2]{#scale[0.6]{"; st << "#rho_{#mu}("<< int(refDist/utl::m)<<") = " << std::fixed << std::setprecision(2) << nMuRef << "#pm" << nMuRefErr; st << ", #beta = " << std::fixed << std::setprecision(1) << beta << "#pm"; if( betaErr == 0.5 ) st << std::fixed << std::setprecision(1) << 0.0 << " (fixed)"; else st << std::fixed << std::setprecision(1) << betaErr; st <<"}}"; st << "}"; //TLatex * latexNMu = new TLatex( refDist/utl::m-100, nMuRef*(5), st.str().c_str() ); TLatex * latexNMu = new TLatex( 500, 50, st.str().c_str() ); cout << "\nRho_mu("<< int(refDist/utl::m)<<") = " << std::fixed << std::setprecision(2) << nMuRef << "+/-" << nMuRefErr; cout << ", Beta = " << std::fixed << std::setprecision(1) << beta << "+/-"; if( betaErr == 0.5 ) cout << std::fixed << std::setprecision(1) << 0.0 << " (fixed)"; else cout << std::fixed << std::setprecision(1) << betaErr << "\n"; std::stringstream canvasName; canvasName << "MEvent_MLDF_" << "EventNum_" << std::setfill('0') << std::setw(4) << std::fixed << fEventNumber << "_RecE_" << std::fixed << std::setprecision(2) << log10(rEnergy/utl::eV); TCanvas canvas(canvasName.str().c_str() ,"MEventLDF", 200, 10, 600, 600); TPad pad1("pad1", "", 0, 0, 1, 1); pad1.Draw(); pad1.cd(); pad1.SetLogy(); st.str(""); st << "#font[102]{#color[4]{#scale[0.7]{" << "Shw " << std::setfill('0') << std::setw(4) << std::fixed << fEventNumber << ", " << "E=" << std::fixed << std::setprecision(3) << showerE/utl::EeV << " EeV, " << "#theta="<< std::fixed << std::setprecision(1) << showerZenith/utl::deg << "#circ, " << "W="<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns." << "}}}"; grMuons.SetTitle( st.str().c_str() ); grMuons.SetMarkerColor(kBlack); grMuons.SetMarkerStyle(kFullCircle); grMuons.SetMarkerSize(0.8); grMLDFFit.SetTitle( st.str().c_str() ); grMLDFFit.SetMarkerColor(kRed); grMLDFFit.SetMarkerStyle(kFullCircle); grMLDFFit.SetMarkerSize(0.5); grMLDFFit.SetLineColor(kRed); grMLDFFit.SetLineWidth( 2 ); grMLDFFit.SetLineStyle(kSolid); grMLDFFit.GetXaxis()->SetTitle( "r [m]"); grMLDFFit.GetYaxis()->SetTitle( "#rho_{#mu} [m^{-2}]"); grMLDFFit.GetXaxis()->CenterTitle( ); grMLDFFit.GetYaxis()->CenterTitle( ); grMLDFFit.GetXaxis()->SetRangeUser( 0, 2500 ); grMLDFFit.GetYaxis()->SetRangeUser( 0.01, 1000 ); grMLDFFit.GetXaxis()->SetNdivisions(505, kTRUE); /* the axis will be optimized around 5 primary and 5 secondary divisions. */ // grMLDFFit.GetXaxis()->SetNdivisions(504, kFALSE); // the axis will be forced to use exactly 4 primary and 5 secondary divisions. /* n = n1 + 100*n2 + 10000*n3, where n1 is the number of primary divisions, n2 is the number of second order divisions and n3 is the number of third order divisions. */ //grMLDFFit.GetXaxis()->SetRangeUser( 1e-3, 5e2 ); /* grSilent.SetMarkerColor(kBlue); grSilent.SetMarkerStyle(2); grSilent.SetMarkerSize(1.5); grSilent.SetLineColor(kBlue); */ // grMLDFFit.SetMinimum( 3e-2 ); // grMLDFFit.SetMaximum( 15e3 ); grMLDFFit.Draw("AL"); grMuons.Draw("P"); // grSilent.Draw("P>"); os.str(""); os << "LDF Fit"; leg.AddEntry( &grMLDFFit, os.str().c_str(),"L"); for (int i = 0; i < arSilent.GetEntriesFast(); ++i ) { TArrow* arrow = (TArrow*)arSilent.UncheckedAt(i); arrow->SetAngle(45); arrow->SetLineWidth(1); arrow->SetLineColor(kBlue); arrow->Draw(); TLine* line = (TLine*)liSilent.UncheckedAt(i); line->SetLineWidth(2); line->SetLineColor(kBlue); line->Draw(); } for (int i = 0; i < arSaturated.GetEntriesFast(); ++i ) { TArrow* arrow = (TArrow*)arSaturated.UncheckedAt(i); arrow->SetAngle(45); arrow->SetLineWidth(1); arrow->SetLineColor(kBlue); arrow->Draw(); TLine* line = (TLine*)liSaturated.UncheckedAt(i); line->SetLineWidth(2); line->SetLineColor(kBlue); line->Draw(); } latexNMu->Draw(); latexMLDF->Draw(); /* os.str(""); os << "/home/juan/Offline/Examples/MdEventPlotter_3_offline_file"; // os << "/" << showerModel; // os << "/" << showerPrimary; // os << "/" << fRun; // os << "/" << std::fixed << std::setprecision(2) << log10(showerE/utl::eV); // os << "/a"<< std::setfill('0') << std::setw(2) << std::fixed << std::setprecision(0) << showerZenith/utl::deg; os << "/"; os << "Mean_MLDF_" << showerModel << "_"<< showerPrimary; os << "_LogE_" << std::fixed << std::setprecision(2) << log10(showerE/utl::eV) << "eV"; os << "_Zenith_" << std::fixed << std::setprecision(0) << showerZenith/utl::deg << "deg"; os << ".root"; std::cout << "RETRIEVING SIM SHOWER FILE " << os.str() << endl; TFile iRootFile( os.str().c_str() ); iRootFile.cd(); TList * list =iRootFile.GetListOfKeys(); for( int l=1; l <= (list->GetEntries()-2)/2; ++l ){ os.str(""); os << "Underground MLDF " << l; if( l == showerNumber ) { cout << "TH1D " << os.str() << endl; //iRootFile.ls(); TH1D* Udensity = (TH1D*) iRootFile.Get( os.str().c_str() ); if( ! Udensity ) { cerr << "Udensity " << l << " not retrieved\n"; return eFailure; } //double AreaDet = 30*cos( showerZenith );//m2 Udensity->Scale(fAreaDet); os.str(""); Udensity->SetTitle( os.str().c_str() ); Udensity->SetMarkerColor(kGray+2); Udensity->SetMarkerStyle(kFullCircle); Udensity->SetMarkerSize(0.5); Udensity->SetLineColor(kGray+2); Udensity->SetLineWidth( 1 ); Udensity->SetLineStyle(kSolid); Udensity->GetXaxis()->SetTitle( "r (m)"); Udensity->GetYaxis()->SetTitle( "#rho_{#mu}(r) (m^{-2})"); Udensity->GetYaxis()->SetRangeUser( 1e-3, 5e2 ); Udensity->GetXaxis()->CenterTitle( ); Udensity->GetYaxis()->CenterTitle( ); Udensity->SetStats(0); Udensity->Draw("SAME"); os.str(""); os << "Injected"; leg.AddEntry( Udensity, os.str().c_str(),"L"); fNmuSizes.push_back(boost::make_tuple( Udensity->GetBinContent( Udensity->FindBin(refDist) ) , nMuRef )); break; } } */ leg.SetEntrySeparation( 0.8 ); leg.SetX1NDC( 0.5 ); leg.SetX2NDC( 1 - 1.1*fStyle->GetPadRightMargin() ); leg.SetY1NDC( 1 - 1.6*fStyle->GetPadTopMargin() ); leg.SetY2NDC( 1 - 1.1*fStyle->GetPadTopMargin() ); leg.SetBorderSize( 0 ); leg.SetTextFont( fStyle->GetTitleFont() ); leg.SetTextSize( fStyle->GetTitleSize()-0.005 ); leg.SetFillColor(canvas.GetFillColor()); leg.Draw(); // canvasName << fPlotExtension1; canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension1).c_str() ); // JMF 14-Nov-2012 canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension2).c_str() ); // JMF 14-Nov-2012 canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension3).c_str() ); // JMF 14-Nov-2012 } return eSuccess; } VModule::ResultFlag MdEventPlotter::PlotEvent(evt::Event& event) { fStyle->cd(); const mdet::MDetector& mdDe = det::Detector::GetInstance().GetMDetector(); typedef mdet::MDetector::CounterConstIterator Itc; if (!event.HasSimShower()) { ERROR("Current event does not have a simulated shower."); return eFailure; } const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower() ; // const utl::Point& rCore = sRecShower.GetCorePosition(); const double rEnergy = sRecShower.GetEnergy(); const evt::ShowerSimData& simShower = event.GetSimShower(); const Point& corePos = simShower.GetPosition(); // const TimeStamp& coreTime = simShower.GetTimeStamp(); const double showerE = simShower.GetEnergy(); const double showerZenith = simShower.GetZenith(); const double showerAzimuth = simShower.GetAzimuth(); const int showerNumber = simShower.GetShowerNumber(); fRunNumber = showerNumber; std::string showerPrimary, showerModel; if ( simShower.GetPrimaryParticle() == Particle::eProton ) showerPrimary = "pr"; if ( simShower.GetPrimaryParticle() == Particle::eIron ) showerPrimary = "fe"; //showerModel = "qgsjetII"; cout << " EVENT NUMBER -> " << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << fEventNumber << endl; cout << " SHOWER PRIMARY -> " << showerPrimary << endl; cout << " SHOWER ENERGY -> " << std::fixed << std::setprecision(2) << log10(showerE/eV) << endl; cout << " SHOWER ZENITH ANGLE -> " << std::setfill('0') << std::setw(2) << std::fixed << std::setprecision(1) << showerZenith/utl::deg << endl; cout << " SHOWER AZIMUTH ANGLE -> " << std::setfill('0') << std::setw(2) << std::fixed << std::setprecision(1) << showerAzimuth/utl::deg << endl; cout << " SHOWER MODEL -> " << showerModel << endl; cout << " SHOWER NUMBER -> " << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << showerNumber << endl << endl; INFO( "\nLOOPING OVER MEVENT DATA"); if (event.HasMEvent()) { mevt::MEvent& me = event.GetMEvent(); TGraph grCounters; TGraph grBarycenter; // TGraph grCountersWithSig; unsigned int ncounters = me.GetNumberOfCounters(); TObjArray LatexArray(ncounters); TObjArray grCountersWithSigArray; if ( !(event.GetRecShower().HasMRecShower()) ) { ERROR("\n---\nCurrent event does not have a MD reconstructed shower\n---\n"); return eFailure; } unsigned int pt, ptwS; pt= ptwS = 0 ; for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd(); ic != ec; ++ic ) { if ( ic->GetRecStatus() == mevt::RecStatus::eRejected ) { // JMF 03-May-2013 // JMF 03-May-2013 cout << "Counter " << ic->GetId() << " REJECTED, IT IS NOT IN THE ARRAY "; // JMF 03-May-2013 cout << " --\n"; // JMF 03-May-2013 continue; // JMF 03-May-2013 } // JMF 03-May-2013 //if(!ic->HasRecData()) { //cout << "Counter " << ic->GetId() << " DOES NOT HAVE REC DATA "; // JMF 03-May-2013 //cout << " --\n"; // JMF 03-May-2013 //} const mdet::Counter& mdCo = mdDe.GetCounter( ic->GetId() ); // GetPosition actually is from Station obj which is a friend of Counter obj const Point countPos = mdCo.GetPosition(); const UTMPoint utm( countPos, ReferenceEllipsoid::GetWGS84()); if( ic->HasRecData() ) fWindowSize = ic->GetRecData().GetWindowSize(); grCounters.SetPoint( pt++, utm.GetEasting() /utl::km, utm.GetNorthing()/utl::km ); //bool hastrace = false; unsigned int nChPerCounter = 0; unsigned int nPartPerCounterInj = 0; unsigned int nPartPerCounterRec = 0; for (mevt::Counter::ModuleIterator im = ic->ModulesBegin(), em = ic->ModulesEnd(); im != em; ++im) { // const mdet::Module& mdMo = mdCo.GetModule( im->GetId() ); for (mevt::Module::ChannelIterator ch = im->ChannelsBegin(), ech = im->ChannelsEnd(); ch != ech; ++ch) { // if ( ch->HasTrace()) { const utl::TraceB& trace = ch->GetTrace(); bool atleastonesample = false; for (utl::TraceB::SizeType s = 0; s < trace.GetSize(); ++s) { //cout << bool(trace[s]) << " "; if( (atleastonesample = bool(trace[s]) ) ) break; } if( atleastonesample ) { //hastrace = true; ++nChPerCounter; const unsigned int scId = mdMo.ChannelToScintillatorId( ch->GetId() ); if( im->HasScintillator( scId ) ) { //cout << "RETRIEVING SCINTILLATOR " << scId << " SIM DATA\n"; //JMF unsigned int injp , recp; injp = recp = 0; const mevt::Scintillator& sc = im->GetScintillator( scId ); if (sc.HasSimData()) { const mevt::ScintillatorSimData& ssd = sc.GetSimData(); //injp = ssd.GetNumberOfParticles(); injp = ssd.GetNumberOfInjectedMuons(); nPartPerCounterInj += injp; } if (sc.HasRecData()) { const mevt::ScintillatorRecData& srd = sc.GetRecData(); recp = srd.GetNumberOfMuons(); nPartPerCounterRec += recp; } /* cout << "Number of Part\n" << injp << " injected -- " << recp << " reconstructed " << "\n"; cout << "Counter Id: " << ic->GetId() << endl; cout << "Module Id: " << im->GetId() << endl; cout << "Channel Id: " << ch->GetId() << endl; */ //JMF } //cout << "Northing " << utm.GetNorthing()/utl::km << " "; //cout << "Easting " << utm.GetEasting() /utl::km << " "; //cout << endl; } } } }//end modules loop //Add Counter with signal (hastrace = kTrue) or at leat three muon (nPartPerCounterRec >= 3) if( nPartPerCounterRec >= 3 ) { TMarker *gCountersWithSig = new TMarker ( utm.GetEasting() /utl::km, // x utm.GetNorthing()/utl::km, // y kFullCircle ); /////////////// // To plot 'particles injected' change nPartPerCounterRec to // nPartPerCounterInj in SetMarkerSize( ) 10 lines below /////////////// double m, b, SizeMin; //******* unsigned int nPartPerCounterRecMax, nPartPerCounterRecMin, SizeMax; //* nPartPerCounterRecMax = 1000; //* nPartPerCounterRecMin = 3; //* SizeMax = 5; //* This is for scaling the SizeMin = 1; //* size of the markers m = (SizeMax-SizeMin)/(nPartPerCounterRecMax-nPartPerCounterRecMin); //* b = SizeMax - m*nPartPerCounterRecMax; //* gCountersWithSig->SetMarkerSize(m * nPartPerCounterRec + b); //******* gCountersWithSig->SetMarkerColor(kRed); grCountersWithSigArray.Add( gCountersWithSig ); ptwS += 1; std::stringstream ps; ps.str(""); ps << "#font[102]{"; ps << "#color[4]{#scale[0.55]{"<SetTitle( "Easting [km]"); grCounters.GetYaxis()->SetTitle( "Northing [km]"); grCounters.GetXaxis()->CenterTitle( ); grCounters.GetYaxis()->CenterTitle( ); grCounters.Draw("AP"); // Plotting counters with signal for (int i = 0; i < grCountersWithSigArray.GetEntriesFast(); ++i ) { TMarker* marker = (TMarker*)grCountersWithSigArray.UncheckedAt(i); marker->Draw(); } // Plotting text for (int i = 0; i < LatexArray.GetEntriesFast(); ++i ) { TLatex* tex = (TLatex*)LatexArray.UncheckedAt(i); tex->Draw(); } // Plotting barycenter const evt::ShowerMRecData& mrecShower = event.GetRecShower().GetMRecShower(); const Point& baryCenter = mrecShower.GetBarycenter(); const UTMPoint utmBary( baryCenter, ReferenceEllipsoid::GetWGS84()); grBarycenter.SetPoint( 0 , utmBary.GetEasting() /utl::km, utmBary.GetNorthing()/utl::km ); grBarycenter.SetMarkerColor(kBlue); grBarycenter.SetMarkerStyle(kOpenCircle); grBarycenter.SetMarkerSize(1); grBarycenter.Draw("P"); // Plotting core position const UTMPoint utmcore( corePos, ReferenceEllipsoid::GetWGS84()); double l = 0.6; TArrow ar( utmcore.GetEasting()/utl::km ,// x1 utmcore.GetNorthing()/utl::km,// y1 utmcore.GetEasting()/utl::km + l*cos(showerAzimuth/utl::rad), // x2 utmcore.GetNorthing()/utl::km + l*sin(showerAzimuth/utl::rad), // y2 0.02,"<" ); ar.SetAngle(40); ar.SetLineWidth(2); ar.SetLineColor(kBlue); ar.Draw(); // Core positions are stored for plotting them in 'PlotAllCores' function fCores.push_back(boost::make_tuple( utmcore.GetEasting(), utmcore.GetNorthing() )); // canvasName << fPlotExtension1; canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension1 ).c_str() ); canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension2 ).c_str() ); canvas.Print( (fOutDir + "/" + canvasName.str() + fPlotExtension3 ).c_str() ); } } else { INFO("\n++++\nEvent without MEvent: nothing to be done."); } // End event.HasMEvent() return eSuccess; } VModule::ResultFlag MdEventPlotter::Finish() { PlotAllCores(); PlotDistanceDistributions(); PlotNmuDistribution(); ofstream & DensityVsDistance = *fFile1; INFO("\nLateral distribution function:\n" "Rho = Distance Core-Counter (perpendicular distance from axis shower)\n" "Dens = Muon density\n" "DensErr = Muon density error" ); TabularStream tab("|.|.|.|.|"); tab << hline << "Event" << endc << "Rho" << endc << "Dens" << endc << "DensErr" << endr << "Num" << endc << "[m]" << endc << "[1/m^2]" << endc << "[1/m^2]" << endr << hline; /* Printing file with all the LDFs */ for( LDFCounterIt sIt = LDFCandidates.begin(), sEnd = LDFCandidates.end(); sIt != sEnd; ++sIt ) { int fEvNum = (*sIt).get<0>(); double Rho = (*sIt).get<1>(); double MdSig = (*sIt).get<2>(); double MdSigErr = (*sIt).get<3>(); tab << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << fEvNum << endc; tab << std::fixed << std::setprecision(1) << Rho/m << endc; tab << ' ' << std::fixed << std::setprecision(3) << MdSig << endc; tab << ' ' << std::fixed << std::setprecision(3) << MdSigErr << endr; DensityVsDistance << std::setfill('0') << std::setw(4) << std::fixed << std::setprecision(0) << fEvNum << '\t'; DensityVsDistance << std::fixed << std::setprecision(1) << Rho/m << '\t'; DensityVsDistance << std::fixed << std::setprecision(3) << MdSig << '\t'; DensityVsDistance << std::fixed << std::setprecision(3) << MdSigErr << endl; } tab << delr << hline ; cerr << tab; delete fFile1; return eSuccess; } void MdEventPlotter::PlotAllCores() { fStyle->cd(); TGraph grCounters; TGraph grCores; const mdet::MDetector& d = det::Detector::GetInstance().GetMDetector(); typedef mdet::MDetector::CounterConstIterator Itc; int pt; pt = 0; for (Itc ic = d.AllCountersBegin(); ic != d.AllCountersEnd(); ++ic) { //if ( ic->GetRecStatus() == mdet::RecStatus::eRejected ) { // JMF 03-May-2013 // JMF 03-May-2013 //cout << "In Counter " << ic->GetId() << " REJECTED, IT IS NOT IN THE ARRAY "; // JMF 03-May-2013 //cout << " --\n"; // JMF 03-May-2013 //continue; // JMF 03-May-2013 //} // JMF 03-May-2013 const mdet::Counter& mdCo = *ic; const Point countPos = mdCo.GetPosition(); const UTMPoint utm( countPos, ReferenceEllipsoid::GetWGS84()); grCounters.SetPoint( pt++, utm.GetEasting()/utl::km, utm.GetNorthing()/utl::km ); } pt = 0; for(CoordIter coresIt = fCores.begin(), ce = fCores.end(); coresIt != ce ; coresIt++ ){ double e = (*coresIt).get<0>(); double n = (*coresIt).get<1>(); grCores.SetPoint(pt++, e/utl::km, n/utl::km ); } std::stringstream st; st << "#font[102]{#color[4]{#scale[0.65]{" << "Cores distribution, " << "W="<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns." << "}}}"; grCounters.SetTitle( st.str().c_str() ); grCounters.SetMarkerColor(kBlack); grCounters.SetMarkerStyle(kFullCircle); grCounters.SetMarkerSize(0.8); grCounters.GetXaxis()->SetTitle( "Easting (km)"); grCounters.GetYaxis()->SetTitle( "Northing (km)"); grCounters.GetXaxis()->CenterTitle( ); grCounters.GetYaxis()->CenterTitle( ); TCanvas canvas("Cores" ,"Cores", 200, 10, 600, 600); TPad pad1("pad1", "", 0, 0, 1, 1); pad1.Draw(); pad1.cd(); grCounters.Draw("AP"); grCores.SetMarkerColor(kRed); grCores.SetMarkerStyle(kOpenCircle); grCores.SetMarkerSize(0.3); grCores.Draw("P"); st.str(""); st << "CoresDistribution_W"<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns"; // st << fPlotExtension1; canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension1 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension2 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension3 ).c_str() ); return; } void MdEventPlotter::PlotDistanceDistributions() { fStyle->cd(); std::stringstream st; double rmin, rmax; int rbins; rmax = 3500; rmin = 0; double dr = 50; rbins = int((rmax-rmin)/dr); st.str(""); st << "Candidate"; TH1D hCandidate( st.str().c_str(), st.str().c_str(), rbins, rmin, rmax ) ; st.str(""); st << "Silent"; TH1D hSilent( st.str().c_str(), st.str().c_str(), rbins, rmin, rmax ) ; st.str(""); st << "Saturated"; TH1D hSaturated( st.str().c_str(), st.str().c_str(), rbins, rmin, rmax ) ; st.str(""); st << "Rejected"; TH1D hRejected( st.str().c_str(), st.str().c_str(), rbins, rmin, rmax ) ; for( PosIter pIt = fCandidates.begin(), pe = fCandidates.end(); pIt != pe; ++pIt ) hCandidate.Fill( *pIt/utl::m ); for( PosIter pIt = fSilents.begin(), pe = fSilents.end(); pIt != pe; ++pIt ) hSilent.Fill( *pIt/utl::m ); for( PosIter pIt = fRejected.begin(), pe = fRejected.end(); pIt != pe; ++pIt ) hRejected.Fill( *pIt/utl::m ); for( PosIter pIt = fSaturated.begin(), pe = fSaturated.end(); pIt != pe; ++pIt ) hSaturated.Fill( *pIt/utl::m ); st.str(""); st << "#font[102]{#color[4]{#scale[0.7]{" << "Core to Counters distance distribution, " << "W="<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns." << "}}}"; TCanvas canvas("CoreCountersDistance" ,"CoreCountersDistance", 200, 10, 600, 600); TPad pad1("pad1", "", 0, 0, 1, 1); pad1.Draw(); pad1.cd(); THStack hstack; hCandidate.Scale(1./hCandidate.Integral() ); hSilent.Scale(1./hSilent.Integral() ); hRejected.Scale(1./hRejected.Integral() ); hSaturated.Scale(1./hSaturated.Integral() ); hCandidate.SetTitle( st.str().c_str() ); hCandidate.SetLineColor(kBlack); hCandidate.SetLineWidth(2); hCandidate.SetMarkerColor(kBlack); hCandidate.SetMarkerStyle(kFullCircle); hCandidate.SetMarkerSize(0.8); hCandidate.GetXaxis()->SetTitle( "r [m]"); hCandidate.GetYaxis()->SetTitle( "Prob."); hCandidate.GetXaxis()->CenterTitle( ); hCandidate.GetYaxis()->CenterTitle( ); hCandidate.SetStats( kFALSE ); //hCandidate.Draw(); hSilent.SetLineColor(kRed); hSilent.SetLineWidth(2); hSilent.SetMarkerColor(kRed); hSilent.SetMarkerStyle(kOpenCircle); hSilent.SetMarkerSize(0.3); //hSilent.Draw("same"); hRejected.SetLineColor(kGray); hRejected.SetLineWidth(2); hRejected.SetMarkerColor(kGray); hRejected.SetMarkerStyle(kOpenCircle); hRejected.SetMarkerSize(0.3); //hRejected.Draw("same"); hSaturated.SetLineColor(kBlue); hSaturated.SetLineWidth(2); hSaturated.SetMarkerColor(kBlue); hSaturated.SetMarkerStyle(kOpenCircle); hSaturated.SetMarkerSize(0.3); //hSaturated.Draw("same"); hstack.Add( &hCandidate ); hstack.Add( &hSilent ); hstack.Add( &hRejected ); hstack.Add( &hSaturated ); hstack.Draw("nostack"); hstack.GetXaxis()->SetTitle( "r [m]"); hstack.GetYaxis()->SetTitle( "Prob."); hstack.GetXaxis()->CenterTitle( ); hstack.GetYaxis()->CenterTitle( ); TLegend leg; leg.SetNColumns( 1 ); st.str(""); st << " = " << std::fixed << std::setprecision(0) << hCandidate.GetMean() << " m"; leg.AddEntry( &hCandidate, st.str().c_str(),"L"); st.str(""); st << " = " << std::fixed << std::setprecision(0) << hSilent.GetMean() << " m"; leg.AddEntry( &hSilent, st.str().c_str(),"L"); st.str(""); st << " = " << std::fixed << std::setprecision(0) << hRejected.GetMean() << " m"; leg.AddEntry( &hRejected, st.str().c_str(),"L"); st.str(""); st << " = " << std::fixed << std::setprecision(0) << hSaturated.GetMean() << " m"; leg.AddEntry( &hSaturated, st.str().c_str(),"L"); leg.SetEntrySeparation( 0.8 ); leg.SetX1NDC( 0.3 ); leg.SetX2NDC( 1 - 1.05*fStyle->GetPadRightMargin() ); leg.SetY1NDC( 1 - 2.5*fStyle->GetPadTopMargin() ); leg.SetY2NDC( 1 - 1.5*fStyle->GetPadTopMargin() ); leg.SetBorderSize( 0 ); leg.SetTextFont( fStyle->GetTitleFont() ); leg.SetTextSize( fStyle->GetTitleSize()+0.002 ); leg.SetFillColor(canvas.GetFillColor()); leg.Draw(); leg.Draw(); st.str(""); st << "CoreCountersDistanceDistribution_W"<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns"; // st << fPlotExtension1; canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension1 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension2 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension3 ).c_str() ); return; } void MdEventPlotter::PlotNmuDistribution() { fStyle->cd(); std::stringstream st; double nmax = 100; double nmin = -nmax; double dn = 1e0; int nbins = int((nmax-nmin)/dn); st.str(""); st << "NmuSizeDiff"; TH1D hNmuSizeDiff( st.str().c_str(), st.str().c_str(), nbins, nmin, nmax ) ; for(NmuIter sizesIt = fNmuSizes.begin(), ce = fNmuSizes.end(); sizesIt != ce ; sizesIt++ ){ double ninj = (*sizesIt).get<0>(); double nrec = (*sizesIt).get<1>(); hNmuSizeDiff.Fill( (1-nrec/ninj)*100 ); } st.str(""); st << "#font[102]{#color[4]{#scale[0.7]{" << "N_{inj}- N_{rec} distribution, " << "W="<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns." << "}}}"; TCanvas canvas("NmuSizeDist" ,"NmuSizeDist", 200, 10, 600, 600); TPad pad1("pad1", "", 0, 0, 1, 1); pad1.Draw(); pad1.cd(); hNmuSizeDiff.Scale(1./hNmuSizeDiff.Integral() ); hNmuSizeDiff.SetTitle( st.str().c_str() ); hNmuSizeDiff.SetLineColor(kBlack); hNmuSizeDiff.SetLineWidth(2); hNmuSizeDiff.SetMarkerColor(kBlack); hNmuSizeDiff.SetMarkerStyle(kFullCircle); hNmuSizeDiff.SetMarkerSize(0.8); st.str(""); st <<"(N_{inj}- N_{rec})/N_{inj} (%)"; hNmuSizeDiff.GetXaxis()->SetTitle( st.str().c_str() ); hNmuSizeDiff.GetYaxis()->SetTitle( "Prob."); hNmuSizeDiff.GetXaxis()->CenterTitle( ); hNmuSizeDiff.GetYaxis()->CenterTitle( ); hNmuSizeDiff.SetStats( kFALSE ); hNmuSizeDiff.Draw(); TLatex * latex; st.str(""); st << "#font[102]{"; st << "#color[2]{#scale[0.65]{"; st << "#nu= " << std::fixed << std::setprecision(1) << hNmuSizeDiff.GetMean() << "%"; st <<"}}"; st << "}"; latex = new TLatex( 40, hNmuSizeDiff.GetMaximum()*0.9, st.str().c_str() ); latex->Draw(); st.str(""); st << "#font[102]{"; st << "#color[2]{#scale[0.65]{"; st << "#sigma= " << std::fixed << std::setprecision(1) << hNmuSizeDiff.GetRMS() << "%"; st <<"}}"; st << "}"; latex = new TLatex( 40, hNmuSizeDiff.GetMaximum()*0.85, st.str().c_str() ); latex->Draw(); st.str(""); st << "NmuDiffDistribution_W"<< std::fixed << std::setprecision(0) << fWindowSize/utl::ns << "ns"; // st << fPlotExtension1; canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension1 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension2 ).c_str() ); canvas.Print( (fOutDir + "/" + st.str() + fPlotExtension3 ).c_str() ); return; } double MdEventPlotter::CounterProjectedArea(const mdet::Counter& c, const mevt::Counter & evtc, const utl::Point& p, const utl::Vector& v, const utl::CoordinateSystemPtr siteCS) { typedef mdet::Counter::ModuleConstIterator Itm; typedef mdet::Module::ScintillatorConstIterator Its; const utl::Line dir( p, v );//dummy so far double area = 0; for (Itm im = c.ModulesBegin(); im != c.ModulesEnd(); ++im) { const mdet::Module& m = *im; if ( ! evtc.HasModule( m.GetId() ) ) continue; const mevt::Module& evtm = evtc.GetModule( m.GetId() ); if ( evtm.GetRecStatus() == mevt::RecStatus::eRejected ) continue; for (Its is = m.ScintillatorsBegin(); is != m.ScintillatorsEnd(); ++is) { const mdet::Scintillator& s = *is; area += s.ComputeEffectiveArea( dir ); } } double costheta = v.GetCosTheta(siteCS); cout << "Counter " << c.GetId() << " Area " << area/utl::m2 << "--> projected " << area*costheta/utl::m2 << endl; return area*costheta; } void MdEventPlotter::InitializeStyle() { Style_t font=102; Float_t fsize=0.035; Option_t *axis="XYZ"; Float_t margin = 0.15; fStyle = new TStyle; //Canvas style fStyle->SetCanvasBorderMode(0); fStyle->SetCanvasColor(0); //Axis style fStyle->SetTickLength( 0.01, "XYZ" ); fStyle->SetLabelFont( font , axis); fStyle->SetLabelSize( fsize, axis); fStyle->SetLabelOffset( 0.00, axis); fStyle->SetTitleFont( font , axis); fStyle->SetTitleSize( fsize, axis); fStyle->SetTitleOffset( 1.0, axis); fStyle->SetTitleXOffset( 1.0 ); fStyle->SetTitleYOffset( 1.5 ); fStyle->SetTitleAlign( 22 ); fStyle->SetTitleBorderSize( 0 ); fStyle->SetTitleColor( 1, axis ); fStyle->SetTitleFillColor( 0 ); fStyle->SetTitleH( 0 ); fStyle->SetTitleW( 0 ); fStyle->SetTitleX( 0.5 ); fStyle->SetTitleY( 1-margin/2 ); fStyle->SetTitleTextColor( kBlue ); //fStyle->SetTitlePS(const char* pstitle); //fStyle->SetTitleStyle(Style_t style = 1001); fStyle->SetTitleXSize(0.045); fStyle->SetTitleYSize(0.045); fStyle->SetPalette(1,0); fStyle->SetStatColor(0); //Pad style fStyle->SetPadBottomMargin(margin); fStyle->SetPadLeftMargin (margin); fStyle->SetPadRightMargin (margin); fStyle->SetPadTopMargin (margin); fStyle->SetPadBorderMode( 0 ); fStyle->SetPadBorderSize( 0 ); fStyle->SetPadColor( 0 ); fStyle->SetPadGridX(false); fStyle->SetPadGridY(false); fStyle->SetPadTickX(true); fStyle->SetPadTickY(true); //Frame style fStyle->SetFrameBorderMode( 0 ); fStyle->SetFrameBorderSize( 0 ); fStyle->SetFrameFillColor( 0 ); fStyle->SetFrameFillStyle( 1 ); fStyle->SetFrameLineColor( 0 ); fStyle->SetFrameLineStyle( 0 ); fStyle->SetFrameLineWidth( 0 ); } }