main.cpp 14.9 KB
Newer Older
Bryan Berthou's avatar
refs#16    
Bryan Berthou committed
1
2
#include <ElementaryUtils/logger/CustomException.h>
#include <ElementaryUtils/logger/LoggerManager.h>
Pawel Sznajder's avatar
Pawel Sznajder committed
3
4
5
6
7
8
9
10
11
12
#include <ElementaryUtils/parameters/Parameters.h>
#include <partons/beans/convol_coeff_function/ConvolCoeffFunctionResult.h>
#include <partons/beans/convol_coeff_function/GAM2/GAM2ConvolCoeffFunctionKinematic.h>
#include <partons/beans/gpd/GPDType.h>
#include <partons/beans/List.h>
#include <partons/beans/PerturbativeQCDOrderType.h>
#include <partons/beans/PolarizationType.h>
#include <partons/modules/convol_coeff_function/ConvolCoeffFunctionModule.h>
#include <partons/modules/convol_coeff_function/GAM2/GAM2CFFStandard.h>
#include <partons/modules/gpd/GPDGK16.h>
13
14
15
#include <partons/modules/process/GAM2/GAM2ProcessGPSSW21.h>
#include <partons/modules/xi_converter/GAM2/GAM2XiConverterExact.h>
#include <partons/modules/scales/GAM2/GAM2ScalesMgg2Multiplier.h>
Pawel Sznajder's avatar
Pawel Sznajder committed
16
#include <partons/ModuleObjectFactory.h>
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
17
#include <partons/Partons.h>
18
#include <QtCore/qcoreapplication.h>
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
19
20
#include <string>
#include <vector>
21
#include <tuple>
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
22

Pawel Sznajder's avatar
Pawel Sznajder committed
23
using namespace PARTONS;
24

Pawel Sznajder's avatar
Pawel Sznajder committed
25
26
27
/*
 * Parse XML scenarios.
 */
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
28
std::vector<std::string> parseArguments(int argc, char** argv) {
29
    std::vector<std::string> xmlScenarios(argc - 1);
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
30

31
    for (unsigned int i = 1; i < argc; i++) {
32
        xmlScenarios[i - 1] = argv[i];
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
33
34
35
36
37
    }

    return xmlScenarios;
}

Pawel Sznajder's avatar
Pawel Sznajder committed
38
39
40
/*
 * Main function.
 */
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
41
42
int main(int argc, char** argv) {

43
    // Init Qt4
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
44
    QCoreApplication a(argc, argv);
Pawel Sznajder's avatar
Pawel Sznajder committed
45
    Partons* pPartons = 0;
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
46
47

    try {
Pawel Sznajder's avatar
Pawel Sznajder committed
48

Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
49
        // Init PARTONS application
Pawel Sznajder's avatar
Pawel Sznajder committed
50
        pPartons = Partons::getInstance();
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
51
52
        pPartons->init(argc, argv);

Pawel Sznajder's avatar
Pawel Sznajder committed
53
54
55
56
57
58
59
60
61
62
        // Create GPD module with the BaseModuleFactory
        GPDModule* pGPDModule =
                Partons::getInstance()->getModuleObjectFactory()->newGPDModule(
                        GPDGK16::classId);

        // Create CFF module with the BaseModuleFactory
        GAM2ConvolCoeffFunctionModule* pGAM2CFFModule =
                static_cast<GAM2ConvolCoeffFunctionModule*>(Partons::getInstance()->getModuleObjectFactory()->newModuleObject(
                        GAM2CFFStandard::classId));

63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
        // Create XiConverterModule
        PARTONS::GAM2XiConverterModule* pXiConverterModule =
                PARTONS::Partons::getInstance()->getModuleObjectFactory()->newGAM2XiConverterModule(
                        PARTONS::GAM2XiConverterExact::classId);

        // Create ScalesModule
        PARTONS::GAM2ScalesModule* pScalesModule =
                PARTONS::Partons::getInstance()->getModuleObjectFactory()->newGAM2ScalesModule(
                        PARTONS::GAM2ScalesMgg2Multiplier::classId);

        // Set its lambda parameter, so MuF2 = MuR2 = lambda * Mgg2
        pScalesModule->configure(
                ElemUtils::Parameter(
                        PARTONS::GAM2ScalesMgg2Multiplier::PARAMETER_NAME_LAMBDA,
                        1.));

        // Create ProcessModule
        PARTONS::GAM2ProcessModule* pProcessModule =
                PARTONS::Partons::getInstance()->getModuleObjectFactory()->newGAM2ProcessModule(
                        PARTONS::GAM2ProcessGPSSW21::classId);

Pawel Sznajder's avatar
Pawel Sznajder committed
84
85
86
        // Create parameters to configure later DVCSCFFModel with PerturbativeQCD = LO
        ElemUtils::Parameters parameters(
                PerturbativeQCDOrderType::PARAMETER_NAME_PERTURBATIVE_QCD_ORDER_TYPE,
87
                PerturbativeQCDOrderType::NLO);
Pawel Sznajder's avatar
Pawel Sznajder committed
88
89
90
91
92
93

        // Configure DVCSCFFModule with previous parameters.
        pGAM2CFFModule->configure(parameters);

        // Link modules (set physics assumptions of your computation)
        pGAM2CFFModule->setGPDModule(pGPDModule);
94
95
96
        pProcessModule->setScaleModule(pScalesModule);
        pProcessModule->setXiConverterModule(pXiConverterModule);
        pProcessModule->setConvolCoeffFunctionModule(pGAM2CFFModule);
Pawel Sznajder's avatar
Pawel Sznajder committed
97

Pawel Sznajder's avatar
Pawel Sznajder committed
98
        //GPD
Pawel Sznajder's avatar
Pawel Sznajder committed
99
100
        List<GPDType> gpdList;
        gpdList.add(GPDType::H);
Oskar Grocholski's avatar
Oskar Grocholski committed
101
        gpdList.add(GPDType::E);
Pawel Sznajder's avatar
Pawel Sznajder committed
102
        gpdList.add(GPDType::Ht);
Pawel Sznajder's avatar
Pawel Sznajder committed
103

Pawel Sznajder's avatar
Pawel Sznajder committed
104
        //polarization states
Oskar Grocholski's avatar
Oskar Grocholski committed
105
        // X and Y PLUS polarizations !
Pawel Sznajder's avatar
Pawel Sznajder committed
106
107
108
109
110
111
112
113
        std::vector<
                std::tuple<PolarizationType::Type, PolarizationType::Type,
                        PolarizationType::Type> > polarizations;

        polarizations.push_back(
                std::make_tuple(PolarizationType::LIN_TRANS_X_PLUS,
                        PolarizationType::LIN_TRANS_X_PLUS,
                        PolarizationType::LIN_TRANS_X_PLUS));
114
115
116
117
118
119
120
121
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_X_PLUS));
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS));
Pawel Sznajder's avatar
Pawel Sznajder committed
122
123
        polarizations.push_back(
                std::make_tuple(PolarizationType::LIN_TRANS_X_PLUS,
Oskar Grocholski's avatar
Oskar Grocholski committed
124
125
                        PolarizationType::LIN_TRANS_Y_MINUS,
                        PolarizationType::LIN_TRANS_Y_MINUS));
126
127
128
129
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_X_PLUS));
Pawel Sznajder's avatar
Pawel Sznajder committed
130
        polarizations.push_back(
Oskar Grocholski's avatar
Oskar Grocholski committed
131
132
                std::make_tuple(PolarizationType::LIN_TRANS_Y_MINUS,
                        PolarizationType::LIN_TRANS_Y_MINUS,
Pawel Sznajder's avatar
Pawel Sznajder committed
133
134
                        PolarizationType::LIN_TRANS_X_PLUS));
        polarizations.push_back(
Oskar Grocholski's avatar
Oskar Grocholski committed
135
                std::make_tuple(PolarizationType::LIN_TRANS_Y_MINUS,
136
                        PolarizationType::LIN_TRANS_X_PLUS,
Oskar Grocholski's avatar
Oskar Grocholski committed
137
                        PolarizationType::LIN_TRANS_Y_MINUS));
138
139
140
141
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS));
Pawel Sznajder's avatar
Pawel Sznajder committed
142

143
144
//  OLDER CODE
//  u - dependence
Pawel Sznajder's avatar
Pawel Sznajder committed
145

Oskar Grocholski's avatar
Now??    
Oskar Grocholski committed
146
        size_t n = 20;
Pawel Sznajder's avatar
Pawel Sznajder committed
147
148
149
150
        double min = -3.;
        double max = -1.;

        double capS = 20.;
151
        double capMgg2 = 4.;
Pawel Sznajder's avatar
Pawel Sznajder committed
152
153
        double capM = 0.938272013;

Pawel Sznajder's avatar
Pawel Sznajder committed
154
155
156
157
158
159
160
        double xi = (-pow(capM,2) + capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS);
        xi = -((pow(capM,2) - capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS));
        double s = ((pow(capM,2) - capS)*(-pow(capM,2) + capMgg2 - 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);
        s = -((pow(capM,2) - capS)*(pow(capM,2) - capMgg2 + 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);

        double t = - pow(2*xi*capM, 2)/(1. - pow(xi, 2));
        double E = (s - pow(capM, 2)) / (2*capM);
Oskar Grocholski's avatar
Oskar Grocholski committed
161

162
//         Simplified kinematics
Pawel Sznajder's avatar
Pawel Sznajder committed
163
164
165
//        double xi = capMgg2 / (2. * ( capS - std::pow(capM, 2) ) - capMgg2);
//        double t = - 4. * std::pow( xi * capM , 2) / ( 1. - xi*xi );
//        double E = ( capS - std::pow(capM, 2) ) / 2. / capM;
Pawel Sznajder's avatar
Pawel Sznajder committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180

        std::cout << "xi: " << xi << std::endl;
        std::cout << "t: " << t << std::endl;
        std::cout << "E: " << E << std::endl;

        for (size_t i = 0; i <= n; i++) {

            double uPrim = min + i * (max - min) / double(n);

            double Mgg2 = 4.;

            double result = 0.;

            for (size_t j = 0; j < polarizations.size(); j++) {

181
                result +=
Pawel Sznajder's avatar
Pawel Sznajder committed
182
183
184
185
                pProcessModule->compute(
                        std::get<0>(polarizations.at(j)),
                        std::get<1>(polarizations.at(j)),
                        std::get<2>(polarizations.at(j)),
186
                        NumA::Vector3D(0., 0., 0.),
Pawel Sznajder's avatar
Pawel Sznajder committed
187
188
189
                        GAM2ObservableKinematic(t, uPrim, capMgg2, E, 0., 0.), gpdList).getValue().makeSameUnitAs(PhysicalUnit::PB).getValue();
            }

Oskar Grocholski's avatar
Oskar Grocholski committed
190
            std::cout << uPrim << "\t" << result/2. << std::endl; // over 2, not 8
Pawel Sznajder's avatar
Pawel Sznajder committed
191
192
        }

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
// t-DEPENDENCE PART
/*
 * Figure 4 from the 2 photons paper
 * t-dependence, unpolarized differential cross section
 * SgN = 20 GeV^2, u' = -2 GeV^2, Mgg2 = 3 GeV^2
 * t from -1 to -0.1 GeV^2
 */
//        n = 50;
//        min = -1.;
//        max = -0.02618 ;
//
//        capS = 20.;
//        capMgg2 = 3.;
//        capM = 0.938272013;
//
//        xi = (-pow(capM,2) + capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS);
//        xi = -((pow(capM,2) - capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS));
//        s = ((pow(capM,2) - capS)*(-pow(capM,2) + capMgg2 - 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);
//        s = -((pow(capM,2) - capS)*(pow(capM,2) - capMgg2 + 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);
//
////        double t = - pow(2*xi*capM, 2)/(1. - pow(xi, 2));
//        E = (s - pow(capM, 2)) / (2*capM);
//
//        double uPrim = -2.;


        // Simplified kinematics
//        double xi = capMgg2 / (2. * ( capS - std::pow(capM, 2) ) - capMgg2);
//        double t = - 4. * std::pow( xi * capM , 2) / ( 1. - xi*xi );
//        double E = ( capS - std::pow(capM, 2) ) / 2. / capM;

//        std::cout << "t-dependence" << std::endl << "==========================" << std::endl;
//        std::cout << "xi: " << xi << std::endl;
//        std::cout << "u: " << uPrim << std::endl;
//        std::cout << "E: " << E << std::endl;
//
//        for (size_t i = 0; i <= n; i++) {
//
//            t = min + i * (max - min) / double(n);
//
//            double Mgg2 = 3.;
//
//            double result = 0.;
//
//            for (size_t j = 0; j < polarizations.size(); j++) {
//
//                result +=
//                pProcessModule->compute(
//                        std::get<0>(polarizations.at(j)),
//                        std::get<1>(polarizations.at(j)),
//                        std::get<2>(polarizations.at(j)),
//                        NumA::Vector3D(0., 0., 0.),
//                        GAM2ObservableKinematic(t, uPrim, capMgg2, E, 0., 0.), gpdList).getValue().makeSameUnitAs(PhysicalUnit::PB).getValue();
//            }
//
//            std::cout << t << "\t" << result/2. << std::endl; // over 2, not 8
//        }

        // END OF THE t-DEPENDENCE PART


/*
 * Figure 7 from the paper
 * Angular dependence of d sigma / dt du dMgg2 dphi on phi
 * phi is an angle between polarization vector of the initial photon
 * and the transverse component of the outgoing photon 1
 * Sgn = 20, Mgg2 = 3, u = -2 [GeV^2], t = t_min
 */

        //polarization states
        // X and Y PLUS polarizations for outgoing photons
        // Incoming photon -
//        std::vector<
//                std::tuple<PolarizationType::Type, PolarizationType::Type,
//                        PolarizationType::Type> > polarizations;
//
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::UNDEFINED,
//                        PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_X_PLUS));
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::UNDEFINED,
//                        PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_X_PLUS));
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::UNDEFINED,
//                        PolarizationType::LIN_TRANS_X_PLUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS));
//        polarizations.push_back(
//                std::make_tuple(PolarizationType::UNDEFINED,
//                        PolarizationType::LIN_TRANS_Y_MINUS,
//                        PolarizationType::LIN_TRANS_Y_MINUS));
//
//        size_t n = 20;
//        double min = 0.;
//        double max = 2. * 3.14159265;
//
//        double capS = 20.;
//        double capMgg2 = 3.;
//        double capM = 0.938272013;
//
//        double xi = (-pow(capM,2) + capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS);
//        xi = -((pow(capM,2) - capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS)))/(2*pow(capM,2) - capMgg2 + 2*capS));
//        double s = ((pow(capM,2) - capS)*(-pow(capM,2) + capMgg2 - 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);
//        s = -((pow(capM,2) - capS)*(pow(capM,2) - capMgg2 + 3*capS + sqrt(pow(capM,4) + pow(capMgg2 - capS,2) - 2*pow(capM,2)*(capMgg2 + capS))))/(4.*capS);
//
//        double t = - pow(2*xi*capM, 2)/(1. - pow(xi, 2));
//        double E = (s - pow(capM, 2)) / (2*capM);
//
//        double uPrim = -2.;
//
//                for (size_t i = 0; i <= n; i++) {
//
//                    double phi = min + i * (max - min) / double(n);
//
//                    double Mgg2 = capMgg2;
//
//                    double result = 0.;
//
//                    for (size_t j = 0; j < polarizations.size(); j++) {
//
//                        result +=
//                        pProcessModule->compute(
//                                std::get<0>(polarizations.at(j)),
//                                std::get<1>(polarizations.at(j)),
//                                std::get<2>(polarizations.at(j)),
//                                NumA::Vector3D(0., 0., 0.),
//                                GAM2ObservableKinematic(t, uPrim, capMgg2, E, 0., 0.), gpdList).getValue().makeSameUnitAs(PhysicalUnit::PB).getValue();
//                    }
//
//                    std::cout << phi << "\t" << result / 2. / 3.14159265 << std::endl;
//                }
325
326
327

// Remove pointer references
// Module pointers are managed by PARTONS
Pawel Sznajder's avatar
Pawel Sznajder committed
328
329
330
331
332
333
334
        Partons::getInstance()->getModuleObjectFactory()->updateModulePointerReference(
                pGAM2CFFModule, 0);
        pGAM2CFFModule = 0;

        Partons::getInstance()->getModuleObjectFactory()->updateModulePointerReference(
                pGPDModule, 0);
        pGPDModule = 0;
335
336
337
338
339
340
341
342

    }
    // Appropriate catching of exceptions is crucial for working of PARTONS.
    // PARTONS defines its own type of exception, which allows to display class name and function name
    // where the exception has occurred, but also a human readable explanation.
    catch (const ElemUtils::CustomException &e) {

        // Display what happened
Bryan Berthou's avatar
refs#16    
Bryan Berthou committed
343
        pPartons->getLoggerManager()->error(e);
344

Nabil Chouika's avatar
typo    
Nabil Chouika committed
345
        // Close PARTONS application properly
Nabil Chouika's avatar
refs#16    
Nabil Chouika committed
346
347
348
        if (pPartons) {
            pPartons->close();
        }
349
350
351
352
353
    }
    // In a case of standard exception.
    catch (const std::exception &e) {

        // Display what happened
Bryan Berthou's avatar
refs#16    
Bryan Berthou committed
354
        pPartons->getLoggerManager()->error("main", __func__, e.what());
355

Nabil Chouika's avatar
typo    
Nabil Chouika committed
356
        // Close PARTONS application properly
Nabil Chouika's avatar
refs#16    
Nabil Chouika committed
357
358
359
        if (pPartons) {
            pPartons->close();
        }
Bryan Berthou's avatar
refs#16    
Bryan Berthou committed
360
    }
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
361

Nabil Chouika's avatar
typo    
Nabil Chouika committed
362
    // Close PARTONS application properly
Bryan Berthou's avatar
refs#16    
Bryan Berthou committed
363
364
    if (pPartons) {
        pPartons->close();
Bryan Berthou's avatar
refs#16  
Bryan Berthou committed
365
366
367
368
    }

    return 0;
}