test/utils/Density.C

0001 #define Density_cxx
0002 #include "Density.h"
0003 #include <TH2.h>
0004 #include <TStyle.h>
0005 #include <TCanvas.h>
0006
0007 #include <cmath>
0008
0009 //#define DEBUG
0010
0011 //using namespace std;
0012
0013 void Density::Loop()
0014 {
0015   //   In a ROOT session, you can do:
0016   //      Root > .L Density.C
0017   //      Root > Density t
0018   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0019   //      Root > t.Show();       // Show values of entry 12
0020   //      Root > t.Show(16);     // Read and show values of entry 16
0021   //      Root > t.Loop();       // Loop on all entries
0022   //
0023
0024   //     This is the loop skeleton where:
0025   //    jentry is the global entry number in the chain
0026   //    ientry is the entry number in the current Tree
0027   //  Note that the argument to GetEntry must be:
0028   //    jentry for TChain::GetEntry
0029   //    ientry for TTree::GetEntry and TBranch::GetEntry
0030   //
0031   //       To read only selected branches, Insert statements like:
0032   // METHOD1:
0033   //    fChain->SetBranchStatus("*",0);  // disable all branches
0034   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0035   // METHOD2: replace line
0036   //    fChain->GetEntry(jentry);       //read all branches
0037   //by  b_branchname->GetEntry(ientry); //read only this branch
0038   if (fChain == 0) return;
0039
0040   Long64_t nentries = fChain->GetEntriesFast();
0041
0042   cout << " Entries " << nentries << endl;
0043
0044   Long64_t nbytes = 0, nb = 0;
0045   for (Long64_t jentry=0; jentry<nentries; jentry++) {
0046     // load tree variables
0047     Long64_t ientry = LoadTree(jentry);
0048     if (ientry < 0) break;
0049     nb = fChain->GetEntry(jentry);   nbytes += nb;
0050     //
0051
0052     double pathSum     = 0.0; // x_i
0053     double normDensSum = 0.0; // rho_i x x_i
0054     //    cout << " Steps " << Nsteps << endl;
0055
0056     // loop on steps
0057     for(Int_t iStep=0; iStep<Nsteps; iStep++) {
0058       //      cout << iStep << endl;
0059       // x_i
0060       double x_i = sqrt(
0061             (FinalX[iStep]-InitialX[iStep]) * (FinalX[iStep]-InitialX[iStep])
0062             +
0063             (FinalY[iStep]-InitialY[iStep]) * (FinalY[iStep]-InitialY[iStep])
0064             +
0065             (FinalZ[iStep]-InitialZ[iStep]) * (FinalZ[iStep]-InitialZ[iStep])
0066             );
0067       // rho_i
0068       double rho_i = MaterialDensity[iStep];
0069       /*
0070     cout << "################################" << endl;
0071     cout << "\t x_i = "   << x_i   << " mm"    << endl;
0072     cout << "\t rho_i = " << rho_i << " g/cm3" << endl;
0073     cout << "################################" << endl;
0074       */
0075       pathSum     += x_i;
0076       normDensSum += (rho_i * x_i);
0077       //
0078     } // step loop
0079
0080     if(Nsteps!=0) {
0081       // average density: Sum(x_i x rho_i) / Sum(x_i)
0082       double averageDensity = normDensSum / pathSum;
0083       /*
0084     cout << "################################"         << endl;
0085     cout << "\t eta = "              << ParticleEta    << endl;
0086     cout << "\t Sum(x_i) = "         << pathSum        << " mm"         << endl;
0087     cout << "\t Sum(x_i x rho_i) = " << normDensSum    << " mm x g/cm3" << endl;
0088     cout << "\t <rho> = "            << averageDensity << " g/cm3"      << endl;
0089     cout << "################################"         << endl;
0090       */
0091       prof_density_vs_eta->Fill(fabs(ParticleEta),averageDensity);
0092       //
0093     }
0094
0095   } // event loop
0096
0097   // Draw plots
0098   MakePlots("Gigi");
0099   //
0100
0101 }
0102
0103
0104 void Density::MakePlots(TString suffix) {
0105   //
0106   TGaxis::SetMaxDigits(3);
0107   gStyle->SetOptStat(0);
0108   gStyle->SetOptFit(0);
0109   gStyle->SetOptLogy(0);
0110   //
0111   //
0112   // canvas
0113   TCanvas can_Gigi("can_Gigi","can_Gigi",1300,800);
0114   //  can_Gigi.Range(0,0,25,25);
0115   can_Gigi.SetFillColor(kWhite);
0116   //
0117
0118   // Draw
0119   can_Gigi.cd();
0120   prof_density_vs_eta->SetMarkerColor(kBlue);
0121   prof_density_vs_eta->SetMarkerStyle(20);
0122   prof_density_vs_eta->Draw("E1");
0123   //
0124
0125   // Store
0126   can_Gigi.Update();
0127   can_Gigi.SaveAs( Form("%s/AverageDensity_%s.eps",  theDirName.Data(), suffix.Data() ) );
0128   can_Gigi.SaveAs( Form("%s/AverageDensity_%s.gif",  theDirName.Data(), suffix.Data() ) );
0129   //
0130 }
0131
0132 void Density::helpfulCommands() {
0133   cout << "########################################" << endl;
0134   cout << "a.Loop()"                                 << endl;
0135   cout << "########################################" << endl;
0136 }