STARLIGHT/starlight/utils/analyse.cxx

   1 #include "analyse.h"
   2 #include <TMath.h>
   3 #include <iostream>
   4
   5 using namespace std;
   6 Analyse::Analyse() :
   7   fInfile("slight.out"),
   8   fNEvents(1)
   9 {
  10   //Constructor
  11
  12   //Creating histograms
  13   fPtEl = new TH1F("PtEl", "Transverse momentum e+/e-", 100, 0, 2.);
  14   fPtMu = new TH1F("PtMu", "Transverse momentum mu+/mu-", 100, 0, 2.);
  15   fPtPi = new TH1F("PtPi", "Transverse momentum pi+/pi-", 100, 0, 2.);
  16
  17   fRapEl = new TH1F("RapEl", "Rapidity e+/e-", 200, -10, 10);
  18   fRapMu = new TH1F("RapMu", "Rapidity mu+/mu-", 200, -10, 10);
  19   fRapPi = new TH1F("RapPi", "Rapidity pi+/pi-", 200, -10, 10);
  20
  21   fInvMassEl = new TH1F("InvMassEl", "Invariant mass", 100, 0, 5);
  22   fInvMassMu = new TH1F("InvMassMu", "Invariant mass", 100, 0, 5);
  23   fInvMassPi = new TH1F("InvMassPi", "Invariant mass", 100, 0, 5);
  24
  25   fPt1 = new TH1F("fPt1", "Transverse momentum track 1", 100, 0, 2.);
  26   fPt2 = new TH1F("fPt2", "Transverse momentum track 2", 100, 0, 2.);
  27 }
  28
  29 Analyse::Analyse(char* infile, Int_t nEvents) :
  30   fInfile(infile),
  31   fNEvents(nEvents)
  32 {
  33   //Special constructor
  34
  35   //Creating histograms
  36   fPtEl = new TH1F("PtEl", "Transverse momentum e+/e-", 100, 0, 2.);
  37   fPtMu = new TH1F("PtMu", "Transverse momentum mu+/mu-", 100, 0, 2.);
  38   fPtPi = new TH1F("PtPi", "Transverse momentum pi+/pi-", 100, 0, 2.);
  39
  40   fRapEl = new TH1F("RapEl", "Rapidity e+/e-", 200, -10, 10);
  41   fRapMu = new TH1F("RapMu", "Rapidity mu+/mu-", 200, -10, 10);
  42   fRapPi = new TH1F("RapPi", "Rapidity pi+/pi-", 200, -10., 10.);
  43
  44   fInvMassEl = new TH1F("InvMassEl", "Invariant mass", 100, 0, 5);
  45   fInvMassMu = new TH1F("InvMassMu", "Invariant mass", 100, 0, 5);
  46   fInvMassPi = new TH1F("InvMassPi", "Invariant mass", 100, 0., 5.);
  47
  48
  49   fPt1 = new TH1F("fPt1", "Transverse momentum track 1", 100, 0, 2.);
  50   fPt2 = new TH1F("fPt2", "Transverse momentum track 2", 100, 0, 2.);
  51 }
  52
  53 Analyse::~Analyse()
  54 {
  55   //Destructor
  56   delete fPtEl;
  57   delete fPtMu;
  58   delete fPtPi;
  59
  60   delete fRapEl;
  61   delete fRapMu;
  62   delete fRapPi;
  63
  64   delete fInvMassEl;
  65   delete fInvMassMu;
  66   delete fInvMassPi;
  67
  68   delete fPt1;
  69   delete fPt2;
  70
  71
  72 }
  73
  74 Int_t Analyse::Init()
  75 {
  76
  77
  78   cout << "Opening textfile " << fInfile << endl;
  79   if( !(filelist=fopen(fInfile,"r")) ){
  80     cout<<"Couldn't open input file: "<<fInfile<<endl;
  81     return -1;
  82   }
  83   cout << "Done opening textfile" << endl;
  84   return 0;
  85 }
  86
  87 Int_t Analyse::NextEvent()
  88 {
  89   char linelabel[20];
  90   int i1=0;
  91   int i2=0;
  92   int i3=0;
  93   double x1=0.0;
  94   double x2=0.0;
  95   double x3=0.0;
  96   double x4=0.0;
  97   int ntrk=0;
  98   int nvtx=0;
  99   // Event line
 100   fscanf(filelist,"%s %d %d %d ",linelabel,&i1,&ntrk,&i2);
 101   fNParticles = ntrk;
 102   cout<<linelabel<<"  "<<i1<<"  "<<ntrk<<"  "<<i2<<"   "<<fNParticles<<endl;
 103   // Vertex line
 104   fscanf(filelist,"%s %lf %lf %lf %lf %d %d %d %d",
 105          linelabel,&x1,&x2,&x3,&x4,&i1,&i2,&i3,&nvtx);
 106   cout<<linelabel<<"  "<<x1<<"  "<<x2<<"  "<<x3<<"  "<<x4<<"  "<<i1<<"  "<<i2<<"  "<<i3<<"  "<<nvtx<<endl;
 107   if(ntrk != nvtx)cout<<"ERROR: ntrk = "<<ntrk<<"  nvtx = "<<nvtx<<endl;
 108   //
 109   return fNParticles;
 110 }
 111
 112 TParticle* Analyse::NextParticle()
 113 {
 114   char tracklabel[20];
 115   int i1=0;
 116   int i2=0;
 117   int i3=0;
 118   int i4=0;
 119   Int_t idpart = 0;
 120   Double_t px = 0.0;
 121   Double_t py = 0.0;
 122   Double_t pz = 0.0;
 123   Double_t ep = 0.0;
 124
 125   cout<<"In NextParticle: fNparticles = "<<fNParticles<<endl;
 126
 127   // for ( int itk=0; itk < fNParticles; itk++){
 128
 129
 130     fscanf(filelist,"%s %d %le %le %le %d %d %d %d",
 131          tracklabel,&i1,&px,&py,&pz,&i2,&i3,&i4,&idpart);
 132     cout<<"   "<<tracklabel<<"  "<<i1<<"  "<<px<<"  "<<py<<"  "<<pz<<"  "<<i2<<"  "<<i3<<"  "<<i4<<"  "<<idpart<<endl;
 133
 134     TParticle *particle =
 135       new TParticle(idpart, 0, -1, -1, -1, -1, px, py, pz, ep, 0., 0., 0., 0.);
 136
 137     //  }
 138
 139   return particle;
 140 }
 141
 142 void Analyse::doAnalysis()
 143 {
 144
 145   Int_t check = Init();
 146   if(check < 0) return;
 147   //Doing the analysis
 148   for(Int_t ev = 0; ev < fNEvents; ev++){
 149
 150     const Int_t ntracks = NextEvent();
 151     //Array of TLorentzVectors. One vector for each tracks
 152     TLorentzVector* vecArr[ntracks];
 153     TParticle* partArr[ntracks];
 154     //Looping over the tracks of the event
 155     for(Int_t tr = 0; tr < ntracks; tr++){
 156       //Getting a TParticle from the TClonesArray
 157       TParticle *part = NextParticle();
 158       Double_t mpi = 0.13957018;
 159       Double_t energy = TMath::Sqrt(mpi*mpi+part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
 160       //Creating a new TLorentzVector and setting px, py, pz and E.
 161       vecArr[tr] = new TLorentzVector;
 162       vecArr[tr]->SetPxPyPzE(part->Px(), part->Py(), part->Pz(), energy);
 163       cout << "particle " << tr << ": px: " << part->Px() << " py: " << part->Py() << " pz: " << part->Pz() << " Energy: " << energy << endl;
 164       partArr[tr] = part;
 165     }
 166
 167     fPt1->Fill(vecArr[0]->Pt());
 168     fPt2->Fill(vecArr[1]->Pt());
 169
 170     //Creating a new TLorentzVector, which is the sum of the elements in vecArr
 171     TLorentzVector sum;
 172     for(Int_t i = 0; i < ntracks; i++){
 173       sum += *vecArr[i];
 174     }
 175     //Filling the histograms depending on particle type
 176
 177     if(partArr[0]->GetPdgCode() == 11 || partArr[0]->GetPdgCode() == -11){
 178       fPtEl->Fill(sum.Pt());
 179       fRapEl->Fill(sum.Rapidity());
 180       fInvMassEl->Fill(sum.M());
 181     }
 182     else if(partArr[0]->GetPdgCode() == 13 || partArr[0]->GetPdgCode() == -13){
 183       fPtMu->Fill(sum.Pt());
 184       fRapMu->Fill(sum.Rapidity());
 185       fInvMassMu->Fill(sum.M());
 186     }
 187     else if(partArr[0]->GetPdgCode() == 211 || partArr[0]->GetPdgCode() == -211){
 188       fPtPi->Fill(sum.Pt());
 189       cout << "sum.Rapidity: " << sum.Rapidity() << endl;
 190       cout << "sum.M(): " << sum.M() << endl;
 191       fRapPi->Fill(sum.Rapidity());
 192       fInvMassPi->Fill(sum.M());
 193     }
 194
 195   }
 196   //Writing the histograms to file
 197   TFile file("histograms.root", "RECREATE");
 198   fPtEl->Write();
 199   fRapEl->Write();
 200   fInvMassEl->Write();
 201   fPtMu->Write();
 202   fRapMu->Write();
 203   fInvMassMu->Write();
 204   fPtPi->Write();
 205   fRapPi->Write();
 206   fInvMassPi->Write();
 207   fPt1->Write();
 208   fPt2->Write();
 209
 210
 211 }