test/genkine/gen/fastGen.C

   1 // Example: generation of kinematics tree with selected properties.
   2 // Below we select events containing the decays D* -> D0 pi, D0 -> K- pi+
   3 // inside the barrel part of the ALICE detector (45 < theta < 135)
   4
   5 // To be able to compile, please add -I<ALIROOT_INSTALL_PATH>/include
   6 // to your .rootrc
   7 #if !defined(__CINT__) || defined(__MAKECINT__)
   8 #include <Riostream.h>
   9 #include <TH1F.h>
  10 #include <TStopwatch.h>
  11 #include <TDatime.h>
  12 #include <TRandom.h>
  13 #include <TDatabasePDG.h>
  14 #include <TParticle.h>
  15 #include <TArrayI.h>
  16 #include <TTree.h>
  17
  18 #include "AliGenerator.h"
  19 #include "AliPDG.h"
  20 #include "AliRunLoader.h"
  21 #include "AliRun.h"
  22 #include "AliStack.h"
  23 #include "AliHeader.h"
  24 #include "AliGenPythia.h"
  25 #include "AliPythia.h"
  26 #endif
  27
  28 Float_t EtaToTheta(Float_t arg);
  29 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar);
  30
  31 void fastGen(Int_t nev = 1, const char* filename = "galice.root")
  32 {
  33   AliPDG::AddParticlesToPdgDataBase();
  34   TDatabasePDG::Instance();
  35
  36
  37
  38   // Run loader
  39   AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
  40
  41   rl->SetCompressionLevel(2);
  42   rl->SetNumberOfEventsPerFile(nev);
  43   rl->LoadKinematics("RECREATE");
  44   rl->MakeTree("E");
  45   gAlice->SetRunLoader(rl);
  46
  47   //  Create stack
  48   rl->MakeStack();
  49   AliStack* stack = rl->Stack();
  50
  51   //  Header
  52   AliHeader* header = rl->GetHeader();
  53
  54   //  Create and Initialize Generator
  55
  56   // Example of charm generation taken from Config_PythiaHeavyFlavours.C
  57   AliGenPythia *gener = new AliGenPythia(-1);
  58   gener->SetEnergyCMS(14000.);
  59   gener->SetMomentumRange(0,999999);
  60   gener->SetPhiRange(0., 360.);
  61   gener->SetThetaRange(0.,180.);
  62   //  gener->SetProcess(kPyCharmppMNR); // Correct Pt distribution, wrong mult
  63   gener->SetProcess(kPyMb); // Correct multiplicity, wrong Pt
  64   gener->SetStrucFunc(kCTEQ4L);
  65   gener->SetPtHard(2.1,-1.0);
  66   gener->SetFeedDownHigherFamily(kFALSE);
  67   gener->SetStack(stack);
  68   gener->Init();
  69
  70   // Go to galice.root
  71   rl->CdGAFile();
  72
  73   // Forbid some decays. Do it after gener->Init(0, because
  74   // the initialization of the generator includes reading of the decay table.
  75
  76   AliPythia * py= AliPythia::Instance();
  77   py->SetMDME(737,1,0); //forbid D*+->D+ + pi0
  78   py->SetMDME(738,1,0);//forbid D*+->D+ + gamma
  79
  80   // Forbid all D0 decays except D0->K- pi+
  81   for(Int_t d=747; d<=762; d++){
  82     py->SetMDME(d,1,0);
  83   }
  84   // decay 763 is D0->K- pi+
  85   for(Int_t d=764; d<=807; d++){
  86     py->SetMDME(d,1,0);
  87   }
  88
  89   //
  90   //                        Event Loop
  91   //
  92
  93   TStopwatch timer;
  94   timer.Start();
  95   for (Int_t iev = 0; iev < nev; iev++) {
  96
  97     cout <<"Event number "<< iev << endl;
  98
  99     //  Initialize event
 100     header->Reset(0,iev);
 101     rl->SetEventNumber(iev);
 102     stack->Reset();
 103     rl->MakeTree("K");
 104
 105     //  Generate event
 106     Int_t nprim = 0;
 107     Int_t ntrial = 0;
 108     Int_t ndstar = 0;
 109
 110
 111
 112     //-------------------------------------------------------------------------------------
 113
 114     while(!ndstar) {
 115       // Selection of events with D*
 116       stack->Reset();
 117       stack->ConnectTree(rl->TreeK());
 118       gener->Generate();
 119       ntrial++;
 120       nprim = stack->GetNprimary();
 121
 122       for(Int_t ipart =0; ipart < nprim; ipart++){
 123         TParticle * part = stack->Particle(ipart);
 124         if(part)    {
 125
 126           if (TMath::Abs(part->GetPdgCode())== 413) {
 127
 128             TArrayI daughtersId;
 129
 130             GetFinalDecayProducts(ipart,*stack,daughtersId);
 131
 132             Bool_t kineOK = kTRUE;
 133
 134             Double_t thetaMin = TMath::Pi()/4;
 135             Double_t thetaMax = 3*TMath::Pi()/4;
 136
 137             for (Int_t id=1; id<=daughtersId[0]; id++) {
 138               TParticle * daughter = stack->Particle(daughtersId[id]);
 139               if (!daughter) {
 140                 kineOK = kFALSE;
 141                 break;
 142               }
 143
 144               Double_t theta = daughter->Theta();
 145               if (theta<thetaMin || theta>thetaMax) {
 146                 kineOK = kFALSE;
 147                 break;
 148               }
 149             }
 150
 151             if (!kineOK) continue;
 152
 153             part->Print();
 154             ndstar++;
 155
 156           }
 157         }
 158       }
 159     }
 160
 161     cout << "Number of particles " << nprim << endl;
 162     cout << "Number of trials " << ntrial << endl;
 163
 164     //  Finish event
 165     header->SetNprimary(stack->GetNprimary());
 166     header->SetNtrack(stack->GetNtrack());
 167
 168     //      I/O
 169     stack->FinishEvent();
 170     header->SetStack(stack);
 171     rl->TreeE()->Fill();
 172     rl->WriteKinematics("OVERWRITE");
 173
 174   } // event loop
 175   timer.Stop();
 176   timer.Print();
 177
 178   //                         Termination
 179   //  Generator
 180   gener->FinishRun();
 181   //  Write file
 182   rl->WriteHeader("OVERWRITE");
 183   gener->Write();
 184   rl->Write();
 185 }
 186
 187
 188
 189 Float_t EtaToTheta(Float_t arg){
 190   return (180./TMath::Pi())*2.*atan(exp(-arg));
 191 }
 192
 193
 194 void GetFinalDecayProducts(Int_t ind, AliStack & stack , TArrayI & ar){
 195
 196   // Recursive algorithm to get the final decay products of a particle
 197   //
 198   // ind is the index of the particle in the AliStack
 199   // stack is the particle stack from the generator
 200   // ar contains the indexes of the final decay products
 201   // ar[0] is the number of final decay products
 202
 203   if (ind<0 || ind>stack.GetNtrack()) {
 204     cerr << "Invalid index of the particle " << ind << endl;
 205     return;
 206   }
 207   if (ar.GetSize()==0) {
 208     ar.Set(10);
 209     ar[0] = 0;
 210   }
 211
 212   TParticle * part = stack.Particle(ind);
 213
 214   Int_t iFirstDaughter = part->GetFirstDaughter();
 215   if( iFirstDaughter<0) {
 216     // This particle is a final decay product, add its index to the array
 217     ar[0]++;
 218     if (ar.GetSize() <= ar[0]) ar.Set(ar.GetSize()+10); // resize if needed
 219     ar[ar[0]] = ind;
 220     return;
 221   }
 222
 223   Int_t iLastDaughter = part->GetLastDaughter();
 224
 225   for (Int_t id=iFirstDaughter; id<=iLastDaughter;id++) {
 226     // Now search for final decay products of the daughters
 227     GetFinalDecayProducts(id,stack,ar);
 228   }
 229 }