src/framework/PartonShowerGenerator.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #include "PartonShowerGenerator.h"
0017 #include "PartonShower.h"
0018 #include "JetEnergyLoss.h"
0019 #include "JetScapeLogger.h"
0020
0021 #include <iostream>
0022
0023 using namespace std;
0024
0025 namespace Jetscape {
0026
0027 void PartonShowerGenerator::DoShower(JetEnergyLoss &j) {
0028   double tStart = 0;
0029   double currentTime = 0;
0030
0031   VERBOSESHOWER(8) << "Hard Parton from Initial Hard Process ...";
0032   VERBOSEPARTON(6, *j.GetShowerInitiatingParton());
0033
0034   // consider pointers for speed up ...
0035   vector<Parton> pIn;
0036   vector<Parton> pOut;
0037   vector<Parton> pInTemp;
0038   vector<Parton> pOutTemp;
0039   vector<Parton> pInTempModule;
0040
0041   vector<node> vStartVec;
0042   vector<node> vStartVecOut;
0043   vector<node> vStartVecTemp;
0044
0045   node vStart;
0046   node vEnd;
0047
0048   pIn.push_back(*j.GetShowerInitiatingParton());
0049
0050   // Add here the Hard Shower emitting parton ...
0051   vStart = j.GetShower()->new_vertex(make_shared<Vertex>());
0052   vEnd = j.GetShower()->new_vertex(make_shared<Vertex>());
0053   j.GetShower()->new_parton(
0054       vStart, vEnd, make_shared<Parton>(*j.GetShowerInitiatingParton()));
0055
0056   // start then the recursive shower ...
0057   vStartVec.push_back(vEnd);
0058   //vStartVecTemp.push_back(vEnd);
0059
0060   // ISSUE: Probably not yet 100% wrt to time step evolution ...
0061   // Logic mistake to remove the original ones when no split occured !!??? Follow up!!!!
0062   REMARK << "DoShower() Splitting including time evolution (allowing "
0063             "non-splits at later times) implemeted correctly (should be made "
0064             "nicer/pointers). To be checked!!!";
0065
0066   // --------------------------------------------
0067
0068   do {
0069     VERBOSESHOWER(7) << "Current time = " << currentTime << " with #Input "
0070                      << pIn.size();
0071     currentTime += j.GetDeltaT();
0072
0073     // --------------------------------------------
0074
0075     for (int i = 0; i < pIn.size(); i++) {
0076       //DEBUG:
0077       //cout<<currentTime<<" pIn size = "<<pIn.size()<<" "<<i<<" "<<pIn[i].pt()<<endl;
0078
0079       pInTemp.push_back(pIn[i]);
0080       pInTempModule.push_back(pIn[i]);
0081
0082       j.SentInPartons(j.GetDeltaT(), currentTime, pIn[i].pt(), pInTempModule,
0083                       pOutTemp);
0084
0085       vStart = vStartVec[i];
0086       vStartVecTemp.push_back(vStart);
0087
0088       //DEBUG:
0089       //cout<<vStart<<endl;
0090       // --------------------------------------------
0091       for (int k = 0; k < pOutTemp.size(); k++) {
0092         vEnd = j.GetShower()->new_vertex(
0093             make_shared<Vertex>(0, 0, 0, currentTime));
0094         j.GetShower()->new_parton(vStart, vEnd,
0095                                   make_shared<Parton>(pOutTemp[k]));
0096
0097         //DEBUG:
0098         //cout<<vStart<<"-->"<<vEnd<<endl;
0099         //cout<<pOutTemp[k];
0100         //cout<<vStartVec.size()<<endl;
0101         //cout<<pInTempModule.size()<<endl;
0102
0103         vStartVecOut.push_back(vEnd);
0104         pOut.push_back(pOutTemp[k]);
0105
0106         // --------------------------------------------
0107         // Add new roots from ElossModules ...
0108         // (maybe add for clarity a new vector in the signal!???)
0109         // Otherwise keep track of input size (so far always 1
0110         // and check if size > 1 and create additional root nodes to that vertex ...
0111         // Simple Test here below:
0112         // DEBUG:
0113         //cout<<"In JetEnergyloss : "<<pInTempModule.size()<<endl;
0114
0115         if (pInTempModule.size() > 1) {
0116           VERBOSESHOWER(7) << pInTempModule.size() - 1
0117                            << " new root node(s) to be added ...";
0118           //cout<<pInTempModule.size()-1<<" new root node(s) to be added ..."<<endl;
0119
0120           for (int l = 1; l < pInTempModule.size(); l++) {
0121             node vNewRootNode = j.GetShower()->new_vertex(
0122                 make_shared<Vertex>(0, 0, 0, currentTime - j.GetDeltaT()));
0123             j.GetShower()->new_parton(vNewRootNode, vEnd,
0124                                       make_shared<Parton>(pInTempModule[l]));
0125           }
0126         }
0127         // --------------------------------------------
0128
0129         if (k == 0) {
0130           pInTemp.pop_back();
0131           vStartVecTemp.pop_back();
0132         }
0133       }
0134       // --------------------------------------------
0135
0136       pOutTemp.clear();
0137       pInTempModule.clear();
0138     }
0139
0140     // --------------------------------------------
0141
0142     pIn.clear();
0143
0144     pIn.insert(pIn.end(), pInTemp.begin(), pInTemp.end());
0145     pIn.insert(pIn.end(), pOut.begin(), pOut.end());
0146
0147     pOut.clear();
0148     pInTemp.clear();
0149
0150     vStartVec.clear();
0151
0152     vStartVec.insert(vStartVec.end(), vStartVecTemp.begin(),
0153                      vStartVecTemp.end());
0154     vStartVec.insert(vStartVec.end(), vStartVecOut.begin(), vStartVecOut.end());
0155
0156     vStartVecOut.clear();
0157     vStartVecTemp.clear();
0158   } while (currentTime < j.GetMaxT()); //other criteria (how to include; TBD)
0159
0160   // --------------------------------------------
0161
0162   // real ceck using mom. consveration at vertex ...!!!!
0163   //pShower->save(&(cout<<BOLDCYAN));
0164
0165   pIn.clear();
0166   pOut.clear();
0167   pInTemp.clear();
0168   pOutTemp.clear();
0169   vStartVec.clear();
0170 }
0171
0172 } // end namespace Jetscape