TEvtGen/Photos/PhotosBranch.cxx

   1 #include <vector>
   2 #include <list>
   3 #include "PH_HEPEVT_Interface.h"
   4 #include "PhotosParticle.h"
   5 #include "PhotosBranch.h"
   6 #include "Photos.h"
   7 #include "Log.h"
   8 using std::vector;
   9 using std::list;
  10 using std::endl;
  11
  12 namespace Photospp
  13 {
  14
  15 PhotosBranch::PhotosBranch(PhotosParticle* p)
  16 {
  17         daughters = p->getDaughters();
  18
  19         //Suppress if somehow got stable particle
  20         if(daughters.size()==0)
  21         {
  22                 Log::Debug(1)<<"Stable particle."<<endl;
  23                 suppression = 1;
  24                 forcing     = 0;
  25                 particle    = NULL;
  26                 return;
  27         }
  28         else if(daughters.at(0)->getMothers().size()==1)
  29         {
  30                 // Regular case - one mother
  31                 Log::Debug(1)<<"Regular case."<<endl;
  32                 particle  = p;
  33                 mothers   = p->findProductionMothers();
  34         }
  35         else
  36         {
  37                 // Advanced case - branch with multiple mothers - no mid-particle
  38                 Log::Debug(1)<<"Advanced case."<<endl;
  39                 particle  = NULL;
  40                 mothers   = daughters.at(0)->getMothers();
  41         }
  42
  43   //--------------------------------------------------
  44   // Finalize suppression/forcing checks
  45   // NOTE: if user forces decay of specific particle,
  46   //       this overrides any suppresion
  47   //--------------------------------------------------
  48
  49         forcing = checkForcingLevel();
  50         if(!forcing) suppression = checkSuppressionLevel();
  51         else         suppression = 0;
  52
  53         // Even if forced or passed suppression check, we still have to check few things
  54         if(!suppression)
  55         {
  56                 // Check momentum conservation
  57                 suppression=!checkMomentumConservation();
  58                 if(suppression) Log::Warning()<<"Branching ignored due to 4-momentum non conservation"<<endl;
  59
  60                 // Check if advanced case has only one daughter
  61                 if(!particle && daughters.size()==1) suppression=-1;
  62
  63                 // If any of special cases is true, we're not forcing this branch
  64                 if(suppression) forcing=0;
  65         }
  66 }
  67
  68 void PhotosBranch::process()
  69 {
  70         Log::Debug(703)<<"   Processing barcode: "<<( (particle) ? particle->getBarcode() : ( (mothers.size()) ? mothers.at(0)->getBarcode() : -1) )<<endl;
  71         /*
  72         cout<<"Particles send to photos (with barcodes in brackets):"<<endl;
  73         vector<PhotosParticle *> get = getParticles();
  74         for(int i=0;i<(int)get.size();i++) cout<<"ID: "<<get.at(i)->getPdgID()<<" ("<<get.at(i)->getBarcode()<<"), "; cout<<endl;
  75         */
  76         int index = PH_HEPEVT_Interface::set(this);
  77         PH_HEPEVT_Interface::prepare();
  78         PHOTOS_MAKE_C(index);
  79         PH_HEPEVT_Interface::complete();
  80         PH_HEPEVT_Interface::get();
  81         checkMomentumConservation();
  82 }
  83
  84 vector<PhotosParticle *> PhotosBranch::getParticles()
  85 {
  86         vector<PhotosParticle *> ret = mothers;
  87         if(particle) ret.push_back(particle);
  88         ret.insert(ret.end(),daughters.begin(),daughters.end());
  89         return ret;
  90 }
  91
  92 bool PhotosBranch::checkMomentumConservation()
  93 {
  94         if(particle)           return particle->checkMomentumConservation();
  95         if(mothers.size()>0)   return mothers.at(0)->checkMomentumConservation();
  96         return true;
  97 }
  98
  99 vector<PhotosBranch *> PhotosBranch::createBranches(vector<PhotosParticle *> particles)
 100 {
 101         Log::Debug(700)<<"PhotosBranch::createBranches - filtering started"<<endl;
 102         list<PhotosParticle *> list(particles.begin(),particles.end());
 103         vector<PhotosBranch *> branches;
 104
 105         // First - add all forced decays
 106         if(Photos::forceBremList)
 107         {
 108                 std::list<PhotosParticle *>::iterator it;
 109                 for(it=list.begin();it!=list.end();it++)
 110                 {
 111                         PhotosBranch *branch = new PhotosBranch(*it);
 112                         int forcing = branch->getForcingStatus();
 113                         if(forcing)
 114                         {
 115                                 Log::Debug(701)<<" Forced: "<<(*it)->getPdgID()<<" (barcode: "<<(*it)->getBarcode()<<") with forcing status= "<<forcing<<endl;
 116                                 branches.push_back(branch);
 117                                 it = list.erase(it);
 118                                 --it;
 119                                 // If forcing consecutive decays
 120                                 if(forcing==2)
 121                                 {
 122                                         PhotosParticle *p = branch->getDecayingParticle();
 123                                         if(!p)
 124                                         {
 125                                                 if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
 126                                                 else continue;
 127                                         }
 128                                         vector<PhotosParticle *> tree = p->getDecayTree();
 129                                         //Add branches for all particles from the list - max O(n*m)
 130                                         std::list<PhotosParticle *>::iterator it2;
 131                                         for(it2=list.begin();it2!=list.end();it2++)
 132                                         {
 133                                                 for(int i=0;i<(int)tree.size();i++)
 134                                                 {
 135                                                         if(tree.at(i)->getBarcode()==(*it2)->getBarcode())
 136                                                         {
 137                                                                 PhotosBranch *b = new PhotosBranch(*it2);
 138                                                                 branches.push_back(b);
 139                                                                 // If we were to delete our next particle in line
 140                                                                 if(it==it2) --it;
 141                                                                 it2 = list.erase(it2);
 142                                                                 --it2;
 143                                                                 break;
 144                                                         }
 145                                                 }
 146                                         }
 147                                 }
 148                         }
 149                         else delete branch;
 150                 }
 151         }
 152         // Quit if we're suppressing everything
 153         if(Photos::isSuppressed) return branches;
 154         // Now - check if remaining decays are suppressed
 155         while(!list.empty())
 156         {
 157                 PhotosParticle *particle = list.front();
 158                 list.pop_front();
 159                 if(!particle) continue;
 160
 161                 PhotosBranch *branch = new PhotosBranch(particle);
 162                 int suppression = branch->getSuppressionStatus();
 163                 if(!suppression) branches.push_back(branch);
 164                 else
 165                 {
 166                         Log::Debug(702)<<"  Suppressed: "<<particle->getPdgID()<<" (barcode: "<<particle->getBarcode()<<") with suppression status= "<<suppression<<endl;
 167                         //If suppressing consecutive decays
 168                         if(suppression==2)
 169                         {
 170                                 PhotosParticle *p = branch->getDecayingParticle();
 171                                 if(!p)
 172                                 {
 173                                         if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
 174                                         else continue;
 175                                 }
 176                                 vector<PhotosParticle *> tree = p->getDecayTree();
 177                                 //Remove all particles from the list - max O(n*m)
 178                                 std::list<PhotosParticle *>::iterator it;
 179                                 for(it=list.begin();it!=list.end();it++)
 180                                 {
 181                                         for(int i=0;i<(int)tree.size();i++)
 182                                         {
 183                                                 if(tree.at(i)->getBarcode()==(*it)->getBarcode())
 184                                                 {
 185                                                         it = list.erase(it);
 186                                                         --it;
 187                                                         break;
 188                                                 }
 189                                         }
 190                                 }
 191                         }
 192                         delete branch;
 193                         continue;
 194                 }
 195
 196                 //In case we don't have mid-particle erase rest of the mothers from list
 197                 if(!branch->getDecayingParticle())
 198                 {
 199                         vector<PhotosParticle *> mothers = branch->getMothers();
 200                         for(int i=0;i<(int)mothers.size();i++)
 201                         {
 202                                 PhotosParticle *m = mothers.at(i);
 203                                 if(m->getBarcode()==particle->getBarcode()) continue;
 204                                 std::list<PhotosParticle *>::iterator it;
 205                                 for(it=list.begin();it!=list.end();it++)
 206                                         if(m->getBarcode()==(*it)->getBarcode())
 207                                         {
 208                                                 it = list.erase(it);
 209                                                 break;
 210                                         }
 211                         }
 212                 }
 213         }
 214         return branches;
 215 }
 216
 217 int PhotosBranch::checkList(bool forceOrSuppress)
 218 {
 219         vector< vector<int>* > *list = (forceOrSuppress) ? Photos::forceBremList : Photos::supBremList;
 220         if(!list) return 0;
 221
 222         // Can't check without pdgid
 223         int motherID;
 224         if(particle) motherID = particle->getPdgID();
 225         else
 226         {
 227                 if(mothers.size()==0) return 0;
 228                 motherID = mothers.at(0)->getPdgID();
 229         }
 230
 231         // Create list of daughters
 232         vector<int> dID;
 233         for(int j=0;j<(int)daughters.size();j++) dID.push_back(daughters[j]->getPdgID());
 234
 235         vector< vector<int> *> &patternList = *list;
 236
 237         // Check if the mother and list of daughters matches any of the declared patterns
 238         for(int j=0; j<(int)patternList.size();j++)
 239         {
 240                 // Skip patterns that don't have our mother
 241                 if(motherID!=(*patternList[j])[0]) continue;
 242
 243                 // Compare decay daughters with pattern - max O(n*m)
 244                 vector<int> &pattern = *patternList[j];
 245                 bool fullMatch=true;
 246                 for(int k = 1; k<(int)pattern.size()-1; k++)
 247                 {
 248                         bool oneMatch=false;
 249                         for(int l=0;l<(int)dID.size(); l++)
 250                                 if(pattern[k]==dID[l]) { oneMatch=true; break; }
 251                         if(!oneMatch) { fullMatch=false; break; }
 252                 }
 253                 // Check if the matching pattern is set for consecutive suppression
 254                 /*
 255                   Currently minimal matching is used.
 256                   If e.g. 25 -> -15 is suppressed, then 25 -> 15,-15 and 25 -> 15,-15,22 etc. is suppressed too
 257                   For exact matching (then suppress 25 -> 15,-15 ; 25 -> 15,-15,22 etc. must be done separately) uncoment line ...:
 258                 */
 259                 // if(pattern.size()<=2 || (fullMatch && dID.size()==pattern.size()-2) )
 260                 // ...and comment out line:
 261                 if(pattern.size()<=2 || fullMatch)
 262                         return (pattern.back()==1) ? 2 : 1;
 263         }
 264         return 0;
 265 }
 266
 267 } // namespace Photospp