[81] | 1 | #include "Examples/interface/Analysis_Ex.h"
|
---|
| 2 |
|
---|
| 3 | using namespace std;
|
---|
| 4 |
|
---|
| 5 | //******************************Debut de l'analyse****************************************
|
---|
| 6 | //*****************************************************************************************
|
---|
| 7 |
|
---|
| 8 | Analysis_Ex::Analysis_Ex(string CardWithCuts,string LogName)
|
---|
| 9 | {
|
---|
| 10 | string temp_string;
|
---|
| 11 | istringstream curstring;
|
---|
| 12 |
|
---|
| 13 | ifstream fichier_a_lire(CardWithCuts.c_str());
|
---|
| 14 | if(!fichier_a_lire.good()) {
|
---|
| 15 | cout << "DataCardname " << CardWithCuts << " not found, use default values" << endl;
|
---|
| 16 | return;
|
---|
| 17 | }
|
---|
| 18 |
|
---|
| 19 | while (getline(fichier_a_lire,temp_string)) {
|
---|
| 20 | curstring.clear(); // needed when using several times istringstream::str(string)
|
---|
| 21 | curstring.str(temp_string);
|
---|
| 22 | string varname;
|
---|
| 23 | float value;
|
---|
| 24 |
|
---|
| 25 | if(strstr(temp_string.c_str(),"#")) { }//remove comments
|
---|
| 26 | else if(strstr(temp_string.c_str(),"PT_ELEC")){curstring >> varname >> value; PT_ELEC = value;}
|
---|
| 27 | else if(strstr(temp_string.c_str(),"PT_MUON")){curstring >> varname >> value; PT_MUON = value;}
|
---|
| 28 | else if(strstr(temp_string.c_str(),"INV_MASS_LL")){curstring >> varname >> value; INV_MASS_LL = value;}
|
---|
| 29 | }
|
---|
| 30 |
|
---|
| 31 | ofstream f_out(LogName.c_str(),ofstream::app);
|
---|
| 32 |
|
---|
| 33 | f_out<<"*******************************************************************"<<endl;
|
---|
| 34 | f_out << left << setw(30) <<"Cut values used in the analysis: "<<""
|
---|
| 35 | << right << setw(37) <<"-------------------------------------"<<"\n";
|
---|
| 36 | f_out << left <<setw(50) << "Invariant mass of the leptons: "<<""
|
---|
| 37 | << right <<setw(17) << INV_MASS_LL <<"\n";
|
---|
| 38 | f_out<<"*******************************************************************"<<endl;
|
---|
| 39 |
|
---|
| 40 | }
|
---|
| 41 |
|
---|
| 42 | Analysis_Ex::~Analysis_Ex()
|
---|
| 43 | {
|
---|
| 44 | }
|
---|
| 45 |
|
---|
| 46 | void Analysis_Ex::Run(ExRootTreeReader *treeReaderGen, ExRootTreeReader *treeReaderRec, ExRootTreeReader *treeReaderTrig)
|
---|
| 47 | {
|
---|
| 48 | total=0;//initialisation of total number of events
|
---|
| 49 | cut_trig=0;
|
---|
| 50 | cut_1=0;//initialisation of counter for cut 1
|
---|
| 51 | cut_2=0;
|
---|
| 52 | //access the branches************************
|
---|
| 53 | //to get the generator level information
|
---|
| 54 | const TClonesArray *GEN = treeReaderGen->UseBranch("Particle");
|
---|
| 55 |
|
---|
| 56 | //to get the reconstructed level information
|
---|
| 57 | const TClonesArray *JET = treeReaderRec->UseBranch("Jet");
|
---|
| 58 | const TClonesArray *TAUJET = treeReaderRec->UseBranch("TauJet");
|
---|
| 59 | const TClonesArray *PHOTO = treeReaderRec->UseBranch("Photon");
|
---|
| 60 | const TClonesArray *ELEC = treeReaderRec->UseBranch("Electron");
|
---|
| 61 | const TClonesArray *MUON = treeReaderRec->UseBranch("Muon");
|
---|
| 62 | const TClonesArray *TRACKS = treeReaderRec->UseBranch("Tracks");
|
---|
| 63 | const TClonesArray *CALO = treeReaderRec->UseBranch("CaloTower");
|
---|
| 64 |
|
---|
| 65 | //to get the VFD reconstructed level information
|
---|
| 66 | const TClonesArray *ZDC = treeReaderRec->UseBranch("ZDChits");
|
---|
| 67 | const TClonesArray *RP220 = treeReaderRec->UseBranch("RP220hits");
|
---|
| 68 | const TClonesArray *FP420 = treeReaderRec->UseBranch("FP420hits");
|
---|
| 69 |
|
---|
| 70 | //to get the trigger information
|
---|
| 71 | const TClonesArray *TRIGGER = treeReaderTrig->UseBranch("TrigResult");
|
---|
| 72 |
|
---|
| 73 | //*******************************************
|
---|
| 74 |
|
---|
| 75 | //run on the events
|
---|
| 76 | Long64_t entry, allEntries = treeReaderRec->GetEntries();
|
---|
| 77 | cout << "** Chain contains " << allEntries << " events" << endl;
|
---|
| 78 | total=allEntries;
|
---|
| 79 |
|
---|
| 80 | //general information
|
---|
| 81 | float E,Px,Py,Pz;
|
---|
| 82 | float PT,Eta,Phi;
|
---|
| 83 |
|
---|
| 84 | //lepton information
|
---|
| 85 | bool IsolFlag;
|
---|
| 86 |
|
---|
| 87 | //bjet information
|
---|
| 88 | bool Btag;
|
---|
| 89 |
|
---|
| 90 | //Particle level information
|
---|
| 91 | int PID, Status, M1,M2,D1,D2;
|
---|
| 92 | float Charge, T, X, Y, Z, M;
|
---|
| 93 |
|
---|
| 94 | //VFD information
|
---|
| 95 | float S,q2,Tx,Ty;
|
---|
| 96 | int side;
|
---|
| 97 |
|
---|
| 98 | for(entry = 0; entry < allEntries; ++entry)
|
---|
| 99 | {
|
---|
| 100 | treeReaderGen->ReadEntry(entry);//access information of generated information
|
---|
| 101 | treeReaderRec->ReadEntry(entry);//access information of reconstructed information
|
---|
| 102 | treeReaderTrig->ReadEntry(entry);//access information of Trigger information
|
---|
| 103 |
|
---|
| 104 | //*****************************************************
|
---|
| 105 | //Example how to run on the generator level information
|
---|
| 106 | //*****************************************************
|
---|
| 107 | TIter itGen((TCollection*)GEN);
|
---|
| 108 | TRootGenParticle *gen;
|
---|
| 109 | itGen.Reset();
|
---|
| 110 | while( (gen = (TRootGenParticle*) itGen.Next()) )
|
---|
| 111 | {
|
---|
| 112 | PID = gen->PID; // particle HEP ID number
|
---|
| 113 | Status = gen->Status; // particle status
|
---|
| 114 | M1 = gen->M1; // particle 1st mother
|
---|
| 115 | M2 = gen->M2; // particle 2nd mother
|
---|
| 116 | D1 = gen->D1; // particle 1st daughter
|
---|
| 117 | D2 = gen->D2; // particle 2nd daughter
|
---|
| 118 | Charge = gen->Charge; // electrical charge
|
---|
| 119 |
|
---|
| 120 | T = gen->T; // particle vertex position (t component)
|
---|
| 121 | X = gen->X; // particle vertex position (x component)
|
---|
| 122 | Y = gen->Y; // particle vertex position (y component)
|
---|
| 123 | Z = gen->Z; // particle vertex position (z component)
|
---|
| 124 | M = gen->M; // particle mass
|
---|
| 125 | }
|
---|
| 126 |
|
---|
| 127 |
|
---|
| 128 | //***********************************************
|
---|
| 129 | //Example how to run on the reconstructed objects
|
---|
| 130 | //***********************************************
|
---|
| 131 |
|
---|
| 132 | //access the Electron branch
|
---|
| 133 | TIter itElec((TCollection*)ELEC);
|
---|
| 134 | TRootElectron *elec;
|
---|
| 135 | itElec.Reset();
|
---|
| 136 | while( (elec = (TRootElectron*) itElec.Next()) )
|
---|
| 137 | {
|
---|
| 138 | E = elec->E; // particle energy in GeV
|
---|
| 139 | Px = elec->Px; // particle momentum vector (x component) in GeV
|
---|
| 140 | Py = elec->Py; // particle momentum vector (y component) in GeV
|
---|
| 141 | Pz = elec->Pz; // particle momentum vector (z component) in GeV
|
---|
| 142 |
|
---|
| 143 | PT = elec->PT; // particle transverse momentum in GeV
|
---|
| 144 | Eta = elec->Eta; // particle pseudorapidity
|
---|
| 145 | Phi = elec->Phi; // particle azimuthal angle in rad
|
---|
| 146 | IsolFlag = elec->IsolFlag; // is the particule isolated?
|
---|
| 147 | }
|
---|
| 148 | //Running on the muon branch is identical:
|
---|
| 149 | TIter itMuon((TCollection*)MUON);
|
---|
| 150 | TRootMuon *muon;
|
---|
| 151 | itMuon.Reset();
|
---|
| 152 | while( (muon = (TRootMuon*) itMuon.Next()) ){}
|
---|
| 153 |
|
---|
| 154 | //access the Photon branch
|
---|
| 155 | TIter itGam((TCollection*)PHOTO);
|
---|
| 156 | TRootPhoton *gam;
|
---|
| 157 | itGam.Reset();
|
---|
| 158 | while( (gam = (TRootPhoton*) itGam.Next()) )
|
---|
| 159 | {
|
---|
| 160 | E = gam->E; // particle energy in GeV
|
---|
| 161 | Px = gam->Px; // particle momentum vector (x component) in GeV
|
---|
| 162 | Py = gam->Py; // particle momentum vector (y component) in GeV
|
---|
| 163 | Pz = gam->Pz; // particle momentum vector (z component) in GeV
|
---|
| 164 |
|
---|
| 165 | PT = gam->PT; // particle transverse momentum in GeV
|
---|
| 166 | Eta = gam->Eta; // particle pseudorapidity
|
---|
| 167 | Phi = gam->Phi; // particle azimuthal angle in rad
|
---|
| 168 | }
|
---|
| 169 |
|
---|
| 170 | //access the jet branch
|
---|
| 171 | TIter itJet((TCollection*)JET);
|
---|
| 172 | TRootJet *jet;
|
---|
| 173 | itJet.Reset();
|
---|
| 174 | while( (jet = (TRootJet*) itJet.Next()) )
|
---|
| 175 | {
|
---|
| 176 | E = jet->E; // particle energy in GeV
|
---|
| 177 | Px = jet->Px; // particle momentum vector (x component) in GeV
|
---|
| 178 | Py = jet->Py; // particle momentum vector (y component) in GeV
|
---|
| 179 | Pz = jet->Pz; // particle momentum vector (z component) in GeV
|
---|
| 180 |
|
---|
| 181 | PT = jet->PT; // particle transverse momentum in GeV
|
---|
| 182 | Eta = jet->Eta; // particle pseudorapidity
|
---|
| 183 | Phi = jet->Phi; // particle azimuthal angle in rad
|
---|
| 184 | Btag = jet->Btag; // is the jet BTagged
|
---|
| 185 | }
|
---|
| 186 | //Running on the tau-jet branch is identical:
|
---|
| 187 | TIter itTaujet((TCollection*)TAUJET);
|
---|
| 188 | TRootTauJet *taujet;
|
---|
| 189 | itTaujet.Reset();
|
---|
| 190 | while( (taujet = (TRootTauJet*) itTaujet.Next()) ){}
|
---|
| 191 |
|
---|
| 192 | //access the track branch
|
---|
| 193 | TIter itTrack((TCollection*)TRACKS);
|
---|
| 194 | TRootTracks *tracks;
|
---|
| 195 | itTrack.Reset();
|
---|
| 196 | while( (tracks = (TRootTracks*) itTrack.Next()) )
|
---|
| 197 | {
|
---|
| 198 | E = tracks->E; // particle energy in GeV
|
---|
| 199 | Px = tracks->Px; // particle momentum vector (x component) in GeV
|
---|
| 200 | Py = tracks->Py; // particle momentum vector (y component) in GeV
|
---|
| 201 | Pz = tracks->Pz; // particle momentum vector (z component) in GeV
|
---|
| 202 |
|
---|
| 203 | PT = tracks->PT; // particle transverse momentum in GeV
|
---|
| 204 | Eta = tracks->Eta; // particle pseudorapidity
|
---|
| 205 | Phi = tracks->Phi; // particle azimuthal angle in rad
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | //Running on the calo branch is identical:
|
---|
| 209 | TIter itCalo((TCollection*)CALO);
|
---|
| 210 | TRootCalo *calo;
|
---|
| 211 | itCalo.Reset();
|
---|
| 212 | while( (calo = (TRootCalo*) itCalo.Next()) ){}
|
---|
| 213 |
|
---|
| 214 | //***************************************************
|
---|
| 215 | //Example how to run on the VFD reconstructed objects
|
---|
| 216 | //***************************************************
|
---|
| 217 |
|
---|
| 218 | //access the ZDC branch
|
---|
| 219 | TIter itZdc((TCollection*)ZDC);
|
---|
| 220 | TRootZdcHits *zdc;
|
---|
| 221 | itZdc.Reset();
|
---|
| 222 | while( (zdc = (TRootZdcHits*) itZdc.Next()) )
|
---|
| 223 | {
|
---|
| 224 | E = zdc->E; // particle energy in GeV
|
---|
| 225 | Px = zdc->Px; // particle momentum vector (x component) in GeV
|
---|
| 226 | Py = zdc->Py; // particle momentum vector (y component) in GeV
|
---|
| 227 | Pz = zdc->Pz; // particle momentum vector (z component) in GeV
|
---|
| 228 |
|
---|
| 229 | PT = zdc->PT; // particle transverse momentum in GeV
|
---|
| 230 | Eta = zdc->Eta; // particle pseudorapidity
|
---|
| 231 | Phi = zdc->Phi; // particle azimuthal angle in rad
|
---|
| 232 |
|
---|
| 233 | T = zdc->T; // time of flight [s]
|
---|
| 234 | side = zdc->side; // -1 or +1
|
---|
| 235 | }
|
---|
| 236 |
|
---|
| 237 | //access the RP220 branch
|
---|
| 238 | TIter itRp220((TCollection*)RP220);
|
---|
| 239 | TRootRomanPotHits *rp220;
|
---|
| 240 | itRp220.Reset();
|
---|
| 241 | //TRootRomanPotHits.MakeClass("test_xav");
|
---|
| 242 | while( (rp220 = (TRootRomanPotHits*) itRp220.Next()) )
|
---|
| 243 | {
|
---|
| 244 | T = rp220->T; // time of flight to the detector [s]
|
---|
| 245 | S = rp220->S; // distance to the IP [m]
|
---|
| 246 | E = rp220->E; // reconstructed energy [GeV]
|
---|
| 247 | q2 = rp220->q2; // reconstructed squared momentum transfer [GeV^2]
|
---|
| 248 |
|
---|
| 249 | X = rp220->X; // horizontal distance to the beam [um]
|
---|
| 250 | Y = rp220->Y; // vertical distance to the beam [um]
|
---|
| 251 |
|
---|
| 252 | Tx = rp220->Tx; // angle of the momentum in the horizontal (x,z) plane [urad]
|
---|
| 253 | Ty = rp220->Ty; // angle of the momentum in the verical (y,z) plane [urad]
|
---|
| 254 |
|
---|
| 255 | side = rp220->side; // -1 or 1
|
---|
| 256 | }
|
---|
| 257 | //running on FP420 branch is identical
|
---|
| 258 | TIter itFp420((TCollection*)FP420);
|
---|
| 259 | TRootRomanPotHits *fp420;
|
---|
| 260 | itFp420.Reset();
|
---|
| 261 | while( (fp420 = (TRootRomanPotHits*) itFp420.Next()) ){}
|
---|
| 262 |
|
---|
| 263 | //*********************************************
|
---|
| 264 | //Example how to run on the trigger information
|
---|
| 265 | //*********************************************
|
---|
| 266 |
|
---|
| 267 | TRootTrigger *trig;
|
---|
| 268 | int NumTrigBit = TRIGGER->GetEntries();
|
---|
| 269 | //get the global response of the trigger
|
---|
| 270 |
|
---|
| 271 | bool GlobalResponse=false;
|
---|
| 272 | if(NumTrigBit!=0)GlobalResponse=true;
|
---|
| 273 | cout<<"GlobalResponse "<<GlobalResponse<<endl;
|
---|
| 274 | for(int i=0; i < NumTrigBit-1; i++){
|
---|
| 275 | trig = (TRootTrigger*)TRIGGER->At(i);
|
---|
| 276 | cout<<"The event has been accepted by the trigger number: "<<trig->Accepted<<endl;
|
---|
| 277 | }
|
---|
| 278 |
|
---|
| 279 | //********************************
|
---|
| 280 | //Example of a very small analysis
|
---|
| 281 | //********************************
|
---|
| 282 |
|
---|
| 283 | TLorentzVector Lept[2];
|
---|
| 284 |
|
---|
| 285 | if(NumTrigBit==0)continue; //event not accepted by the trigger
|
---|
| 286 | cut_trig++;//event accepted
|
---|
| 287 |
|
---|
| 288 | TSimpleArray<TRootElectron> el=SubArrayEl(ELEC,PT_ELEC);//the central isolated electrons, pt > PT_ELEC GeV
|
---|
| 289 | TSimpleArray<TRootMuon> mu=SubArrayMu(MUON,PT_MUON);//the central isolated electrons, pt > PT_MUON GeV
|
---|
| 290 |
|
---|
| 291 | Int_t numElec=el.GetEntries();
|
---|
| 292 |
|
---|
| 293 | if(el.GetEntries()+mu.GetEntries()!=2)continue;//Exactly 2 isolated leptons are needed
|
---|
| 294 | cut_1++;//event accepted
|
---|
| 295 | for(Int_t i=0;i < numElec; i++)Lept[i].SetPxPyPzE(el[i]->Px,el[i]->Py,el[i]->Pz,el[i]->E);
|
---|
| 296 | for(Int_t k = numElec; k < (numElec+mu.GetEntries()); k++)Lept[k].SetPxPyPzE(mu[k-numElec]->Px,mu[k-numElec]->Py,mu[k-numElec]->Pz,mu[k-numElec]->E);
|
---|
| 297 |
|
---|
| 298 | if((Lept[0]+Lept[1]).M() > INV_MASS_LL )continue;// the invariant mass should be < INV_MASS_LL
|
---|
| 299 | cut_2++;//event accepted
|
---|
| 300 | }
|
---|
| 301 |
|
---|
| 302 | }
|
---|
| 303 |
|
---|
| 304 | void Analysis_Ex::WriteOutput(string LogName)
|
---|
| 305 | {
|
---|
| 306 | ofstream f_out(LogName.c_str(),ofstream::app);
|
---|
| 307 |
|
---|
| 308 | f_out<<"*******************************************************************"<<endl;
|
---|
| 309 | f_out << left << setw(20) << "Numer of Events "<<""
|
---|
| 310 | << right << setw(15) << total <<"\n";
|
---|
| 311 | f_out << left << setw(17) << " Accepted by the trigger "<<""
|
---|
| 312 | << right << setw(20) << cut_trig <<"\n";
|
---|
| 313 | f_out << left << setw(17) <<" 2 leptons "<< ""
|
---|
| 314 | << right << setw(20) << cut_1 << ""
|
---|
| 315 | << right << setw(15) << cut_1/total << "\n";
|
---|
| 316 | f_out << left << setw(17) <<" Invariant mass "<< ""
|
---|
| 317 | << right << setw(20) << cut_2 << ""
|
---|
| 318 | << right << setw(15) << cut_2/total << "\n";
|
---|
| 319 | f_out<<"*******************************************************************"<<endl;
|
---|
| 320 | f_out<<" "<<endl;
|
---|
| 321 |
|
---|
| 322 | }
|
---|
| 323 |
|
---|
| 324 | TSimpleArray<TRootElectron> Analysis_Ex::SubArrayEl(const TClonesArray *ELEC,float pt)
|
---|
| 325 | {
|
---|
| 326 | TIter itElec((TCollection*)ELEC);
|
---|
| 327 | TRootElectron *elec;
|
---|
| 328 | itElec.Reset();
|
---|
| 329 | TSimpleArray<TRootElectron> array;
|
---|
| 330 | while( (elec = (TRootElectron*) itElec.Next()) )
|
---|
| 331 | {
|
---|
| 332 | if(elec->PT<pt)continue;
|
---|
| 333 | array.Add(elec);
|
---|
| 334 | }
|
---|
| 335 | return array;
|
---|
| 336 | }
|
---|
| 337 |
|
---|
| 338 | TSimpleArray<TRootMuon> Analysis_Ex::SubArrayMu(const TClonesArray *MUON,float pt)
|
---|
| 339 | {
|
---|
| 340 | TIter itMuon((TCollection*)MUON);
|
---|
| 341 | TRootMuon *muon;
|
---|
| 342 | itMuon.Reset();
|
---|
| 343 | TSimpleArray<TRootMuon> array;
|
---|
| 344 | while( (muon = (TRootMuon*) itMuon.Next()) )
|
---|
| 345 | {
|
---|
| 346 | if(muon->PT<pt)continue;
|
---|
| 347 | array.Add(muon);
|
---|
| 348 | }
|
---|
| 349 | return array;
|
---|
| 350 | }
|
---|
| 351 |
|
---|
| 352 |
|
---|