Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtPDL.cpp
1 //--------------------------------------------------------------------------
2 //
3 // Environment:
4 //      This software is part of the EvtGen package developed jointly
5 //      for the BaBar and CLEO collaborations.  If you use all or part
6 //      of it, please give an appropriate acknowledgement.
7 //
8 // Copyright Information: See EvtGen/COPYRIGHT
9 //      Copyright (C) 1998      Caltech, UCSB
10 //
11 // Module: EvtPDL.cc
12 //
13 // Description: routines to store particle properties in EvtPDL structure.
14 //
15 // Modification history:
16 //
17 //    DJL/RYD     September 25, 1996         Module created
18 //
19 //------------------------------------------------------------------------
20 // 
21 #include "EvtGenBase/EvtPatches.hh"
22 #include <iostream>
23 #include <fstream>
24 #include <stdlib.h>
25 #include <ctype.h>
26 #include <string.h>
27 #include "EvtGenBase/EvtPDL.hh"
28 #include "EvtGenBase/EvtPartProp.hh"
29 #include "EvtGenBase/EvtId.hh"
30 #include "EvtGenBase/EvtParticle.hh"
31 #include "EvtGenBase/EvtReport.hh"
32 using std::endl;
33 using std::fstream;
34 using std::ifstream;
35
36 static int first=1;
37
38 unsigned int EvtPDL::_firstAlias;
39 int EvtPDL::_nentries;
40
41 std::map<std::string,int> EvtPDL::_particleNameLookup;
42
43 EvtPDL::EvtPDL() {
44
45   if (first!=0) { 
46     first=0;
47     _nentries=0;
48     _firstAlias=999999;
49   }
50
51 }
52
53
54 EvtPDL::~EvtPDL(){
55   
56 }
57
58 void EvtPDL::read(const char* fname)
59 {
60   readPDT(fname);
61 }
62
63 void EvtPDL::readPDT(const std::string fname){
64
65
66   ifstream indec;
67   
68   indec.open(fname.c_str());
69
70   char cmnd[100];
71   char xxxx[100];
72
73   char pname[100];
74   int  stdhepid;
75   double mass;
76   double pwidth;
77   double pmaxwidth;
78   int    chg3;  
79   int    spin2;
80   double ctau;
81   int    lundkc;
82   EvtId i;
83
84   if (!indec) {
85     report(ERROR,"EvtGen") << "Could not open:"<<fname.c_str()<<"EvtPDL"<<endl;
86     return;
87   }
88
89   do{
90
91     char ch,ch1;
92
93     do{
94
95       indec.get(ch);  
96       if (ch=='\n') indec.get(ch);
97       if (ch!='*') {
98         indec.putback(ch);
99       }
100       else{
101         while (indec.get(ch1),ch1!='\n');
102       }
103     } while(ch=='*');
104
105     indec >> cmnd;
106
107     if (strcmp(cmnd,"end")){
108
109       if (!strcmp(cmnd,"add")){
110
111         indec >> xxxx;
112         indec >> xxxx;
113         indec >> pname;
114         indec >> stdhepid;
115         indec >> mass;
116         indec >> pwidth;
117         indec >> pmaxwidth;
118         indec >> chg3;  
119         indec >> spin2;
120         indec >> ctau;
121         indec >> lundkc;
122
123
124         i=EvtId(_nentries,_nentries);
125
126         EvtPartProp tmp;
127
128         tmp.setSpinType(EvtSpinType::SCALAR);
129         
130         
131         if (spin2==0) tmp.setSpinType(EvtSpinType::SCALAR);
132         if (spin2==1) tmp.setSpinType(EvtSpinType::DIRAC);
133         if (spin2==2) tmp.setSpinType(EvtSpinType::VECTOR);
134         if (spin2==3) tmp.setSpinType(EvtSpinType::RARITASCHWINGER);
135         if (spin2==4) tmp.setSpinType(EvtSpinType::TENSOR);
136         if (spin2==5) tmp.setSpinType(EvtSpinType::SPIN5HALF);
137         if (spin2==6) tmp.setSpinType(EvtSpinType::SPIN3);
138         if (spin2==7) tmp.setSpinType(EvtSpinType::SPIN7HALF);
139         if (spin2==8) tmp.setSpinType(EvtSpinType::SPIN4);
140         if (spin2==2 && mass < 0.0001 ) tmp.setSpinType(EvtSpinType::PHOTON);
141         if (spin2==1 && mass < 0.0001 ) tmp.setSpinType(EvtSpinType::NEUTRINO);
142         
143         
144         if (!strcmp(pname,"string")){
145           tmp.setSpinType(EvtSpinType::STRING);
146         }
147         
148         if (!strcmp(pname,"vpho")){
149           tmp.setSpinType(EvtSpinType::VECTOR);
150         }
151         
152         
153         tmp.setId(i);
154         tmp.setIdChgConj(EvtId(-1,-1));
155         tmp.setStdHep(stdhepid);
156         tmp.setLundKC(lundkc);
157         tmp.setName(pname);
158         if (_particleNameLookup.find(std::string(pname))!=
159             _particleNameLookup.end()) {
160             report(ERROR,"EvtGen")<<"The particle name:"<<pname<<" is already defined."<<endl;
161             report(ERROR,"EvtGen") << "Will terminate execution.";
162             ::abort();
163         }
164         _particleNameLookup[std::string(pname)]=_nentries;
165         tmp.setctau(ctau);
166         tmp.setChg3(chg3);
167         
168         tmp.initLineShape(mass,pwidth,pmaxwidth);
169
170
171         partlist().push_back(tmp);
172         _nentries++;
173
174       }
175
176       // if find a set read information and discard it
177
178       if (!strcmp(cmnd,"set")){
179
180         indec >> xxxx;
181         indec >> xxxx;
182         indec >> xxxx;
183         indec >> xxxx;
184       }
185
186     }
187
188   }while(strcmp(cmnd,"end"));
189   
190   setUpConstsPdt();
191
192 }
193
194
195 void EvtPDL::aliasChgConj(EvtId a,EvtId abar){
196
197   if (EvtPDL::chargeConj(EvtId(a.getId(),a.getId()))!=
198                  EvtId(abar.getId(),abar.getId())) {
199
200     report(ERROR,"EvtGen")<<"Can't charge conjugate the two aliases:"
201                           <<EvtPDL::name(a).c_str()<<" and "<<EvtPDL::name(abar).c_str()<<endl;
202       
203     ::abort();
204
205   }
206
207   partlist()[a.getAlias()].setIdChgConj(abar);
208   partlist()[abar.getAlias()].setIdChgConj(a);
209
210 }
211
212 EvtId EvtPDL::chargeConj(EvtId id){
213   EvtId idchg=partlist()[id.getAlias()].getIdChgConj();
214
215   if (idchg!=EvtId(-1,-1)) return idchg;
216
217   if (id.getId()!=id.getAlias()){
218     if (chargeConj(EvtId(id.getId(),id.getId()))==EvtId(id.getId(),id.getId())){
219     
220       partlist()[id.getAlias()].setIdChgConj(id);
221       return id;
222     }
223   }
224
225   if (id.getAlias()!=id.getId()) {
226
227     report(ERROR,"EvtGen")<<"Trying to charge conjugate alias particle:"
228                           <<name(id).c_str()<<" without defining the alias!"<<endl;
229       
230     ::abort();
231
232   }
233
234   for (size_t i=0;i<partlist().size();i++){
235     if (partlist()[i].getStdHep()==-partlist()[id.getId()].getStdHep()){
236       partlist()[id.getId()].setIdChgConj(partlist()[i].getId());
237       return partlist()[i].getId();
238     }
239   }
240   
241   partlist()[id.getId()].setIdChgConj(id);
242   return id;
243   
244 }
245
246 EvtId EvtPDL::evtIdFromStdHep(int stdhep){
247
248   for (size_t i=0;i<partlist().size();i++){
249     if (partlist()[i].getStdHep()==stdhep)
250       return partlist()[i].getId();
251   }
252   
253   return EvtId(-1,-1);
254   
255 }
256
257
258
259 void EvtPDL::alias(EvtId num,const std::string& newname){
260   
261   if ( _firstAlias < partlist().size() ) {
262     for(size_t i=_firstAlias;i<partlist().size();i--){
263       if (newname==partlist()[i].getName()){
264         report(WARNING,"EvtGen")<<"Redefining alias:"<<newname.c_str()<<" will be ignored!"<<endl;
265         return;
266       }
267     }
268   }
269   else{
270     _firstAlias=partlist().size();
271   }
272
273   partlist().push_back(partlist()[num.getId()]);
274   int entry=partlist().size()-1;
275   partlist()[entry].setName(newname);
276   if (_particleNameLookup.find(std::string(newname))!=
277       _particleNameLookup.end()){
278             report(ERROR,"EvtGen")<<"The particle name:"<<newname<<" is already defined."<<endl;
279             report(ERROR,"EvtGen") << "Will terminate execution.";
280             ::abort();
281   }
282   _particleNameLookup[std::string(newname)]=entry;
283   partlist()[entry].setId(EvtId(num.getId(),entry));
284   //Lange - Dec7, 2003. Unset the charge conjugate.
285   partlist()[entry].setIdChgConj(EvtId(-1,-1));
286
287 }
288
289 EvtId EvtPDL::getId(const std::string& name ){
290
291   std::map<std::string,int>::iterator it=_particleNameLookup.find(std::string(name));
292   if (it==_particleNameLookup.end()) return EvtId(-1,-1);
293
294   return partlist()[it->second].getId();
295   
296 }
297
298 void EvtPDL::setUpConstsPdt(){
299
300 }
301
302
303 // Function to get EvtId from LundKC ( == Pythia Hep Code , KF ) 
304 EvtId EvtPDL::evtIdFromLundKC(int pythiaId){
305
306   unsigned int i;
307
308   for (i=0;i<partlist().size();i++){
309     if (partlist()[i].getLundKC()==pythiaId)
310       return partlist()[i].getId();
311   }
312   
313   return EvtId(-1,-1);
314   
315 }
316  
317 double EvtPDL::getMeanMass(EvtId i ) { 
318   return partlist()[i.getId()].getMass(); 
319 }
320
321 double EvtPDL::getMass(EvtId i ) {
322   return partlist()[i.getId()].rollMass();
323 }
324
325 double EvtPDL::getRandMass(EvtId i, EvtId *parId, int nDaug, EvtId *dauId, 
326                            EvtId *othDaugId,double maxMass, 
327                            double *dauMasses ) {
328   return partlist()[i.getId()].getRandMass(parId,nDaug,dauId,
329                                           othDaugId,maxMass,dauMasses);
330 }
331
332 double EvtPDL::getMassProb(EvtId i, double mass, double massPar, int nDaug, 
333                            double *massDau) { 
334   return partlist()[i.getId()].getMassProb(mass,massPar,nDaug,massDau);
335 }
336
337 double EvtPDL::getMaxMass(EvtId i ) {
338   return partlist()[i.getId()].getMassMax();
339 }
340
341 double EvtPDL::getMinMass(EvtId i ) { 
342   return partlist()[i.getId()].getMassMin();
343 }
344
345 double EvtPDL::getMaxRange(EvtId i ) {
346   return partlist()[i.getId()].getMaxRange();
347 }
348
349 double EvtPDL::getWidth(EvtId i ) { 
350   return partlist()[i.getId()].getWidth();
351 }
352
353 double EvtPDL::getctau(EvtId i ) {
354   return partlist()[i.getId()].getctau();
355 }
356
357 int EvtPDL::getStdHep(EvtId id ) {
358   return partlist()[id.getId()].getStdHep();
359 }
360
361 int EvtPDL::getLundKC(EvtId id ) {
362   return partlist()[id.getId()].getLundKC();
363 }
364
365 int EvtPDL::chg3(EvtId i ) {
366   return partlist()[i.getId()].getChg3();
367 }
368
369 EvtSpinType::spintype EvtPDL::getSpinType(EvtId i ) {
370   return partlist()[i.getId()].getSpinType();
371 }
372
373 std::string EvtPDL::name(EvtId i) { 
374   return partlist()[i.getAlias()].getName();
375 }
376
377 size_t EvtPDL::entries() { 
378   return partlist().size();
379 }
380
381 EvtId EvtPDL::getEntry(int i) {
382   return partlist()[i].getId();
383 }
384
385 void EvtPDL::reSetMass(EvtId i, double mass) {
386   partlist()[i.getId()].reSetMass(mass);
387 }
388
389 void EvtPDL::reSetWidth(EvtId i, double width) { 
390   partlist()[i.getId()].reSetWidth(width);
391 }
392
393 void EvtPDL::reSetMassMin(EvtId i, double mass) { 
394   partlist()[i.getId()].reSetMassMin(mass);
395 }
396
397 void EvtPDL::reSetMassMax(EvtId i,double mass) { 
398   partlist()[i.getId()].reSetMassMax(mass);
399 }
400
401 void EvtPDL::reSetBlatt(EvtId i,double blatt) {
402   partlist()[i.getId()].reSetBlatt(blatt);
403 }
404
405 void EvtPDL::reSetBlattBirth(EvtId i,double blatt) {
406   partlist()[i.getId()].reSetBlattBirth(blatt);
407 }
408
409 void EvtPDL::includeBirthFactor(EvtId i,bool yesno) {
410   partlist()[i.getId()].includeBirthFactor(yesno);
411 }
412
413 void EvtPDL::includeDecayFactor(EvtId i,bool yesno) {
414   partlist()[i.getId()].includeDecayFactor(yesno);
415 }
416
417 void EvtPDL::changeLS(EvtId i, std::string &newLS ) { 
418   partlist()[i.getId()].newLineShape(newLS);
419 }
420
421 void EvtPDL::setPWForDecay(EvtId i, int spin, EvtId d1, EvtId d2) {  
422   partlist()[i.getId()].setPWForDecay(spin,d1,d2);
423 }
424
425 void EvtPDL::setPWForBirthL(EvtId i, int spin, EvtId par, EvtId othD) {  
426   partlist()[i.getId()].setPWForBirthL(spin,par,othD);
427 }