Changes in modules/PileUpMergerPythia8.cc [2d494a6:dacd9c5] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/PileUpMergerPythia8.cc
r2d494a6 rdacd9c5 29 29 #include "classes/DelphesClasses.h" 30 30 #include "classes/DelphesFactory.h" 31 #include "classes/Delphes TF2.h"31 #include "classes/DelphesFormula.h" 32 32 #include "classes/DelphesPileUpReader.h" 33 33 … … 56 56 57 57 PileUpMergerPythia8::PileUpMergerPythia8() : 58 f Function(0), fPythia(0), fItInputArray(0)58 fPythia(0), fItInputArray(0) 59 59 { 60 fFunction = new DelphesTF2;61 60 } 62 61 … … 65 64 PileUpMergerPythia8::~PileUpMergerPythia8() 66 65 { 67 delete fFunction;68 66 } 69 67 … … 74 72 const char *fileName; 75 73 76 fPileUpDistribution = GetInt("PileUpDistribution", 0);77 78 74 fMeanPileUp = GetDouble("MeanPileUp", 10); 79 80 fZVertexSpread = GetDouble("ZVertexSpread", 0.15); 81 fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09); 82 83 fInputBeamSpotX = GetDouble("InputBeamSpotX", 0.0); 84 fInputBeamSpotY = GetDouble("InputBeamSpotY", 0.0); 85 fOutputBeamSpotX = GetDouble("OutputBeamSpotX", 0.0); 86 fOutputBeamSpotY = GetDouble("OutputBeamSpotY", 0.0); 75 fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3; 87 76 88 77 fPTMin = GetDouble("PTMin", 0.0); 89 90 fFunction->Compile(GetString("VertexDistributionFormula", "0.0"));91 fFunction->SetRange(-fZVertexSpread, -fTVertexSpread, fZVertexSpread, fTVertexSpread);92 78 93 79 fileName = GetString("ConfigFile", "MinBias.cmnd"); … … 100 86 101 87 // create output arrays 102 fParticleOutputArray = ExportArray(GetString("ParticleOutputArray", "stableParticles")); 103 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices")); 88 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles")); 104 89 } 105 90 … … 118 103 TParticlePDG *pdgParticle; 119 104 Int_t pid, status; 120 Float_t x, y, z, t , vx, vy;105 Float_t x, y, z, t; 121 106 Float_t px, py, pz, e; 122 Double_t dz, dphi , dt;123 Int_t numberOfEvents, event, numberOfParticles, i;124 Candidate *candidate , *vertex;107 Double_t dz, dphi; 108 Int_t poisson, event, i; 109 Candidate *candidate; 125 110 DelphesFactory *factory; 126 111 127 const Double_t c_light = 2.99792458E8;128 129 112 fItInputArray->Reset(); 130 131 // --- Deal with primary vertex first ------132 133 fFunction->GetRandom2(dz, dt);134 135 dt *= c_light*1.0E3; // necessary in order to make t in mm/c136 dz *= 1.0E3; // necessary in order to make z in mm137 vx = 0.0;138 vy = 0.0;139 numberOfParticles = fInputArray->GetEntriesFast();140 113 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 141 114 { 142 vx += candidate->Position.X(); 143 vy += candidate->Position.Y(); 144 z = candidate->Position.Z(); 145 t = candidate->Position.T(); 146 candidate->Position.SetZ(z + dz); 147 candidate->Position.SetT(t + dt); 148 fParticleOutputArray->Add(candidate); 149 } 150 151 if(numberOfParticles > 0) 152 { 153 vx /= numberOfParticles; 154 vy /= numberOfParticles; 115 fOutputArray->Add(candidate); 155 116 } 156 117 157 118 factory = GetFactory(); 158 119 159 vertex = factory->NewCandidate(); 160 vertex->Position.SetXYZT(vx, vy, dz, dt); 161 fVertexOutputArray->Add(vertex); 120 poisson = gRandom->Poisson(fMeanPileUp); 162 121 163 // --- Then with pile-up vertices ------ 164 165 switch(fPileUpDistribution) 166 { 167 case 0: 168 numberOfEvents = gRandom->Poisson(fMeanPileUp); 169 break; 170 case 1: 171 numberOfEvents = gRandom->Integer(2*fMeanPileUp + 1); 172 break; 173 default: 174 numberOfEvents = gRandom->Poisson(fMeanPileUp); 175 break; 176 } 177 178 for(event = 0; event < numberOfEvents; ++event) 122 for(event = 0; event < poisson; ++event) 179 123 { 180 124 while(!fPythia->next()); 181 125 182 // --- Pile-up vertex smearing 183 184 fFunction->GetRandom2(dz, dt); 185 186 dt *= c_light*1.0E3; // necessary in order to make t in mm/c 187 dz *= 1.0E3; // necessary in order to make z in mm 188 126 dz = gRandom->Gaus(0.0, fZVertexSpread); 189 127 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); 190 128 191 vx = 0.0; 192 vy = 0.0; 193 numberOfParticles = fPythia->event.size(); 194 for(i = 1; i < numberOfParticles; ++i) 129 for(i = 1; i < fPythia->event.size(); ++i) 195 130 { 196 131 Pythia8::Particle &particle = fPythia->event[i]; … … 208 143 candidate->PID = pid; 209 144 210 candidate->Status = 1;145 candidate->Status = status; 211 146 212 147 pdgParticle = pdg->GetParticle(pid); … … 219 154 candidate->Momentum.RotateZ(dphi); 220 155 221 x -= fInputBeamSpotX; 222 y -= fInputBeamSpotY; 223 candidate->Position.SetXYZT(x, y, z + dz, t + dt); 156 candidate->Position.SetXYZT(x, y, z + dz, t); 224 157 candidate->Position.RotateZ(dphi); 225 candidate->Position += TLorentzVector(fOutputBeamSpotX, fOutputBeamSpotY, 0.0, 0.0);226 158 227 vx += candidate->Position.X(); 228 vy += candidate->Position.Y(); 229 230 fParticleOutputArray->Add(candidate); 159 fOutputArray->Add(candidate); 231 160 } 232 233 if(numberOfParticles > 0)234 {235 vx /= numberOfParticles;236 vy /= numberOfParticles;237 }238 239 vertex = factory->NewCandidate();240 vertex->Position.SetXYZT(vx, vy, dz, dt);241 vertex->IsPU = 1;242 243 fVertexOutputArray->Add(vertex);244 161 } 245 162 } 246 163 247 164 //------------------------------------------------------------------------------ 165
Note:
See TracChangeset
for help on using the changeset viewer.