1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // Class to read events from external (TNtupla) file
18 // Events -> neutron removal by EM dissociation of Pb nuclei
19 // Data from RELDIS code (by I. Pshenichov)
22 #include <TParticle.h>
24 #include <TVirtualMC.h>
25 #include <TDatabasePDG.h>
27 #include "AliGenReaderEMD.h"
31 ClassImp(AliGenReaderEMD)
33 AliGenReaderEMD::AliGenReaderEMD():
77 for(int i=0; i<70; i++){
78 fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
79 fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
81 fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
82 fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
84 fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
85 fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
86 fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
87 fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
88 fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
89 fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
91 fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
92 fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
93 fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
94 fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
99 if(fPcToTrack==kAll) printf("\n\t *** AliGenReaderEMD will track all produced particles \n\n");
100 else if(fPcToTrack==kNotNucleons) printf("\n\t *** AliGenReaderEMD will track all produced particles except nucleons\n\n");
101 else if(fPcToTrack==kOnlyNucleons) printf("\n\t *** AliGenReaderEMD will track only nucleons\n\n");
105 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
106 AliGenReader(reader),
152 // -----------------------------------------------------------------------------------
153 AliGenReaderEMD::~AliGenReaderEMD()
158 // -----------------------------------------------------------------------------------
159 AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
161 // Assignment operator
166 // -----------------------------------------------------------------------------------
167 void AliGenReaderEMD::Copy(TObject&) const
172 Fatal("Copy","Not implemented!\n");
175 // -----------------------------------------------------------------------------------
176 void AliGenReaderEMD::Init()
179 // Reset the existing file environment and open a new root file
183 pFile = new TFile(fFileName);
185 printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
187 fTreeNtuple = (TTree*)gDirectory->Get("h2032");
188 fNcurrent = fStartEvent;
190 TTree *Ntu=fTreeNtuple;
192 // Set branch addresses
194 Ntu->SetBranchAddress("Nleft",&fNnAside);
195 Ntu->SetBranchAddress("Eleft",&fEnAside);
196 Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
197 Ntu->SetBranchAddress("Pxl", fPxnAside);
198 Ntu->SetBranchAddress("Pyl", fPynAside);
199 Ntu->SetBranchAddress("Pzl", fPznAside);
200 Ntu->SetBranchAddress("Nright",&fNnCside);
201 Ntu->SetBranchAddress("Eright",&fEnCside);
202 Ntu->SetBranchAddress("Pxr", fPxnCside);
203 Ntu->SetBranchAddress("Pyr", fPynCside);
204 Ntu->SetBranchAddress("Pzr", fPznCside);
206 Ntu->SetBranchAddress("Nleft_p",&fNpAside);
207 Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
208 Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
209 Ntu->SetBranchAddress("Pxl_p", fPxpAside);
210 Ntu->SetBranchAddress("Pyl_p", fPypAside);
211 Ntu->SetBranchAddress("Pzl_p", fPzpAside);
212 Ntu->SetBranchAddress("Nright_p",&fNpCside);
213 Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
214 Ntu->SetBranchAddress("Pxr_p", fPxpCside);
215 Ntu->SetBranchAddress("Pyr_p", fPypCside);
216 Ntu->SetBranchAddress("Pzr_p", fPzpCside);
218 Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
219 Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
220 Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
221 Ntu->SetBranchAddress("Pxl_pp", fPxppAside);
222 Ntu->SetBranchAddress("Pyl_pp", fPyppAside);
223 Ntu->SetBranchAddress("Pzl_pp", fPzppAside);
224 Ntu->SetBranchAddress("Nright_pp",&fNppCside);
225 Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
226 Ntu->SetBranchAddress("Pxr_pp", fPxppCside);
227 Ntu->SetBranchAddress("Pyr_pp", fPyppCside);
228 Ntu->SetBranchAddress("Pzr_pp", fPzppCside);
230 Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
231 Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
232 Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
233 Ntu->SetBranchAddress("Pxl_pm", fPxpmAside);
234 Ntu->SetBranchAddress("Pyl_pm", fPypmAside);
235 Ntu->SetBranchAddress("Pzl_pm", fPzpmAside);
236 Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
237 Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
238 Ntu->SetBranchAddress("Pxr_pm", fPxpmCside);
239 Ntu->SetBranchAddress("Pyr_pm", fPypmCside);
240 Ntu->SetBranchAddress("Pzr_pm", fPzpmCside);
242 Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
243 Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
244 Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
245 Ntu->SetBranchAddress("Pxl_p0", fPxp0Aside);
246 Ntu->SetBranchAddress("Pyl_p0", fPyp0Aside);
247 Ntu->SetBranchAddress("Pzl_p0", fPzp0Aside);
248 Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
249 Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
250 Ntu->SetBranchAddress("Pxr_p0", fPxp0Cside);
251 Ntu->SetBranchAddress("Pyr_p0", fPyp0Cside);
252 Ntu->SetBranchAddress("Pzr_p0", fPzp0Cside);
254 Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
255 Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
256 Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
257 Ntu->SetBranchAddress("Pxl_et", fPxetaAside);
258 Ntu->SetBranchAddress("Pyl_et", fPyetaAside);
259 Ntu->SetBranchAddress("Pzl_et", fPzetaAside);
260 Ntu->SetBranchAddress("Nright_et",&fNetaCside);
261 Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
262 Ntu->SetBranchAddress("Pxr_et", fPxetaCside);
263 Ntu->SetBranchAddress("Pyr_et", fPyetaCside);
264 Ntu->SetBranchAddress("Pzr_et", fPzetaCside);
266 Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
267 Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
268 Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
269 Ntu->SetBranchAddress("Pxl_om", fPxomegaAside);
270 Ntu->SetBranchAddress("Pyl_om", fPyomegaAside);
271 Ntu->SetBranchAddress("Pzl_om", fPzomegaAside);
272 Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
273 Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
274 Ntu->SetBranchAddress("Pxr_om", fPxomegaCside);
275 Ntu->SetBranchAddress("Pyr_om", fPyomegaCside);
276 Ntu->SetBranchAddress("Pzr_om", fPzomegaCside);
279 // -----------------------------------------------------------------------------------
280 Int_t AliGenReaderEMD::NextEvent()
282 // Read the next event
284 fNparticle = 0; fOffset=0;
286 TFile* pFile = fTreeNtuple->GetCurrentFile();
290 Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
291 if(fNcurrent < nentries) {
292 fTreeNtuple->GetEvent(fNcurrent);
293 if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
295 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
296 nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
298 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
299 nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
300 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
303 printf("\t #### Putting %d particles in the stack\n", nTracks);
304 /*if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t %d+%d neutrons, %d+%d protons\n",
305 fNnAside,fNnCside, fNpAside,fNpCside);
306 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
307 fNppAside,fNppCside,fNpmAside,fNpmCside,
308 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);*/
315 // -----------------------------------------------------------------------------------
316 TParticle* AliGenReaderEMD::NextParticle()
318 // Read the next particle
319 Float_t p[4]={0.,0.,0.,0.};
322 if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
323 if(fNparticle<fNnAside){
324 p[0] = fPxnAside[fNparticle];
325 p[1] = fPynAside[fNparticle];
326 p[2] = fPznAside[fNparticle];
328 // printf(" pc%d n sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
330 else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
331 p[0] = fPxnCside[fNparticle];
332 p[1] = fPynCside[fNparticle];
333 p[2] = fPznCside[fNparticle];
335 // printf(" pc%d n sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
337 else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
338 p[0] = fPxpAside[fNparticle];
339 p[1] = fPypAside[fNparticle];
340 p[2] = fPzpAside[fNparticle];
342 // printf(" pc%d p sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
344 else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
345 p[0] = fPxpCside[fNparticle];
346 p[1] = fPypCside[fNparticle];
347 p[2] = fPzpCside[fNparticle];
349 // printf(" pc%d p sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
351 fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
352 } //**********************************************************************************************
353 if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
354 if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
355 p[0] = fPxppAside[fNparticle];
356 p[1] = fPyppAside[fNparticle];
357 p[2] = fPzppAside[fNparticle];
358 pdgCode = fppPDGCode;
359 // printf(" pc%d pi+ sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
361 if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
362 p[0] = fPxppCside[fNparticle];
363 p[1] = fPyppCside[fNparticle];
364 p[2] = fPzppCside[fNparticle];
365 pdgCode = fppPDGCode;
366 // printf(" pc%d pi+ sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
368 if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
369 p[0] = fPxpmAside[fNparticle];
370 p[1] = fPypmAside[fNparticle];
371 p[2] = fPzpmAside[fNparticle];
372 pdgCode = fpmPDGCode;
373 // printf(" pc%d pi- sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
375 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
376 p[0] = fPxpmCside[fNparticle];
377 p[1] = fPypmCside[fNparticle];
378 p[2] = fPzpmCside[fNparticle];
379 pdgCode = fpmPDGCode;
380 // printf(" pc%d pi- sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
382 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
383 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
384 p[0] = fPxp0Aside[fNparticle];
385 p[1] = fPyp0Aside[fNparticle];
386 p[2] = fPzp0Aside[fNparticle];
387 pdgCode = fp0PDGCode;
388 // printf(" pc%d pi0 sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
390 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside &&
391 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
392 p[0] = fPxp0Cside[fNparticle];
393 p[1] = fPyp0Cside[fNparticle];
394 p[2] = fPzp0Cside[fNparticle];
395 pdgCode = fp0PDGCode;
396 // printf(" pc%d pi0 sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
398 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
399 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
400 p[0] = fPxetaAside[fNparticle];
401 p[1] = fPyetaAside[fNparticle];
402 p[2] = fPzetaAside[fNparticle];
403 pdgCode = fetaPDGCode;
404 // printf(" pc%d eta sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
406 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside &&
407 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
408 p[0] = fPxetaCside[fNparticle];
409 p[1] = fPyetaCside[fNparticle];
410 p[2] = fPzetaCside[fNparticle];
411 pdgCode = fetaPDGCode;
412 // printf(" pc%d eta sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
414 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
415 fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
416 p[0] = fPxomegaAside[fNparticle];
417 p[1] = fPyomegaAside[fNparticle];
418 p[2] = fPzomegaAside[fNparticle];
419 pdgCode = fomegaPDGCode;
420 // printf(" pc%d omega sideA: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
422 if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
423 && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
424 p[0] = fPxomegaCside[fNparticle];
425 p[1] = fPyomegaCside[fNparticle];
426 p[2] = fPzomegaCside[fNparticle];
427 pdgCode = fomegaPDGCode;
428 // printf(" pc%d omega sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
433 Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
434 Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
435 p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
438 Warning("Generate","Particle %d E = %f GeV mass = %f GeV ",pdgCode,p[3],amass);
441 //printf(" Pc %d: PDGcode %d p(%1.2f, %1.2f, %1.2f, %1.3f)\n",
442 // fNparticle,pdgCode,p[0], p[1], p[2], p[3]);
444 TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
445 p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
446 if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);