STARLIGHT/starlight/TStarLight/testsl.C

   1 // -*- C++ -*-
   2 // $Id$
   3
   4 // creates histogram of generated rho0 kinematical data
   5 // writes a text file "rho0.txt" which should be identical to ../../test/rho0.txt
   6
   7 void testsl() {
   8   gSystem->Load("libStarLight");
   9   gSystem->Load("libAliStarLight.so");
  10
  11   TStarLight* sl = new TStarLight("starlight generator", "title", "");
  12
  13 #if 1
  14   sl->SetParameter("BEAM_1_Z     =    1    #Z of projectile");
  15   sl->SetParameter("BEAM_1_A     =    1    #A of projectile");
  16   sl->SetParameter("BEAM_2_Z     =   82    #Z of target");
  17   sl->SetParameter("BEAM_2_A     =  208    #A of target");
  18   sl->SetParameter("BEAM_1_GAMMA = 4264.4  #Gamma of the colliding ions");
  19   sl->SetParameter("BEAM_2_GAMMA = 1682.4  #Gamma of the colliding ions");
  20 #else
  21   sl->SetParameter("BEAM_1_Z     =   82    #Z of projectile");
  22   sl->SetParameter("BEAM_1_A     =  208    #A of projectile");
  23   sl->SetParameter("BEAM_2_Z     =   82    #Z of target");
  24   sl->SetParameter("BEAM_2_A     =  208    #A of target");
  25   sl->SetParameter("BEAM_2_GAMMA = 1470  #Gamma of the colliding ions");
  26   sl->SetParameter("BEAM_1_GAMMA = 1470  #Gamma of the colliding ions");
  27 #endif
  28   sl->SetParameter("W_MAX        =   15    #Max value of w");
  29   sl->SetParameter("W_MIN        =    -1   #Min value of w");
  30   sl->SetParameter("W_N_BINS     =   40    #Bins i w");
  31   sl->SetParameter("RAP_MAX      =    8.   #max y");
  32   sl->SetParameter("RAP_N_BINS   =   80    #Bins i y");
  33   sl->SetParameter("CUT_PT       =    0    #Cut in pT? 0 = (no, 1 = yes)");
  34   sl->SetParameter("PT_MIN       =    1.0  #Minimum pT in GeV");
  35   sl->SetParameter("PT_MAX       =    3.0  #Maximum pT in GeV");
  36   sl->SetParameter("CUT_ETA      =    0    #Cut in pseudorapidity? (0 = no, 1 = yes)");
  37   sl->SetParameter("ETA_MIN      =  -10    #Minimum pseudorapidity");
  38   sl->SetParameter("ETA_MAX      =   10    #Maximum pseudorapidity");
  39   sl->SetParameter("PROD_MODE    =    4    #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), # 4 = incoherent vector meson, 5 = A+A DPMJet single, 6 = A+A DPMJet double, 7 = p+A DPMJet single, 8 = p+A Pythia single )");
  40   sl->SetParameter("PROD_PID     =  113    #Channel of interest (not relevant for photonuclear processes)");
  41   sl->SetParameter("RND_SEED     = 34533   #Random number seed");
  42   sl->SetParameter("BREAKUP_MODE  =   5    #Controls the nuclear breakup");
  43   sl->SetParameter("INTERFERENCE  =   0    #Interference (0 = off, 1 = on)");
  44   sl->SetParameter("IF_STRENGTH   =   1.   #% of intefernce (0.0 - 0.1)");
  45   sl->SetParameter("COHERENT      =   1    #Coherent=1,Incoherent=0");
  46   sl->SetParameter("INCO_FACTOR   =   1.   #percentage of incoherence");
  47   sl->SetParameter("INT_PT_MAX    =   0.24 #Maximum pt considered, when interference is turned on");
  48   sl->SetParameter("INT_PT_N_BINS = 120    #Number of pt bins when interference is turned on");
  49
  50   sl->InitStarLight();
  51   sl->PrintInputs(std::cout);
  52   TClonesArray tca("TParticle", 100);
  53
  54   TLorentzVector v[2], vSum;
  55   TH1* hM  = new TH1D("hM",  "STARLIGHT;M#(){#pi^{+}#pi^{-}}",     100, 0., 2.);
  56   TH1* hPt = new TH1D("hPt", "STARLIGHT;P_{T}#(){#pi^{+}#pi^{-}}", 100, 0., 1.);
  57   TH1* hY  = new TH1D("hY",  "STARLIGHT;Y#(){#pi^{+}#pi^{-}}",     160,-8., 8.);
  58
  59   std::ofstream ofs("rho0.txt");
  60   TParticle *p;
  61   for (Int_t counter(0); counter<1000;) {
  62     sl->GenerateEvent();
  63     sl->BoostEvent();
  64     sl->ImportParticles(&tca, "ALL");
  65     Bool_t isOK = kTRUE;
  66     for (Int_t i=0; i<tca.GetEntries() && isOK; ++i) {
  67       p = (TParticle*)tca.At(i);
  68       p->Momentum(v[i]);
  69 //       isOK = TMath::Abs(v[i].Rapidity()) <= 0.9;
  70     }
  71     tca.Clear();
  72     if (!isOK) continue;
  73     Printf("counter, %d", counter, isOK);
  74     ++counter;
  75     vSum = v[0] + v[1];
  76     ofs << std::fixed << std::setprecision(4)
  77         << vSum.M() << " " << vSum.Perp() << " " << vSum.Rapidity() << " "
  78         << v[0].Eta() << " " << v[0].Px() << " " << v[0].Py() << " " << v[0].Pz() << " "
  79         << v[1].Eta() << " " << v[1].Px() << " " << v[1].Py() << " " << v[1].Pz()
  80         << std::endl;
  81     hM->Fill(vSum.M());
  82     hPt->Fill(vSum.Perp());
  83     hY->Fill(vSum.Rapidity());
  84   }
  85   TFile::Open("sl.root", "RECREATE");
  86   sl->Write();
  87   gFile->Write();
  88
  89   hM->Draw();
  90   c1->SaveAs("SL.pdf(");
  91   hPt->Draw();
  92   c1->SaveAs("SL.pdf");
  93   hY->Draw();
  94   c1->SaveAs("SL.pdf)");
  95 }