EVGEN/AliGenCocktailAfterBurner.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //
  19 // Container class for AliGenerator and AfterBurners
  20 // (which are AliGenerators as well) through recursion.
  21 // The container is itself an AliGenerator a
  22 // what is stored are not the pointers to the generators directly
  23 // but to objects of type
  24 // AliGenCocktailAfterBurner entry.
  25 // The class provides also iterator functionality.
  26 // Author: andreas.morsch@cern.ch and piotr.skowronski@cern.ch
  27 //
  28 // 24.09.2001  Piotr Skowronski
  29 //             debug -> gDebug,
  30 //             fNEvents replaced with gAlice->GetEventsPerRun()
  31 //
  32
  33
  34 #include <Riostream.h>
  35
  36 #include <TList.h>
  37 #include <TObjArray.h>
  38 #include <TParticle.h>
  39
  40 #include "AliGenCocktailAfterBurner.h"
  41 #include "AliGenCocktailEntry.h"
  42 #include "AliStack.h"
  43
  44
  45 ClassImp(AliGenCocktailAfterBurner)
  46 /*********************************************************************/
  47 /*********************************************************************/
  48
  49 AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()
  50 {
  51 // Constructor
  52     if (gDebug > 0)
  53        cout<<"AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()"<<endl;
  54     SetName("AliGenCocktailAfterBurner");
  55     SetTitle("AliGenCocktailAfterBurner");
  56     fInternalStacks =0;
  57     fActiveStack = 0;
  58     fCurrentGenerator = 0;
  59     fCurrentEvent =0;
  60     fAfterBurnerEntries = new TList();
  61     fNAfterBurners = 0;
  62     fGenerationDone = kFALSE;
  63
  64     fActiveEvent = -1;
  65     fNBgEvents = 0;
  66 }
  67 /*********************************************************************/
  68 AliGenCocktailAfterBurner::AliGenCocktailAfterBurner(const AliGenCocktailAfterBurner& in)
  69 {
  70  //cpy ctor
  71 }
  72
  73 /*********************************************************************/
  74
  75 AliGenCocktailAfterBurner::~AliGenCocktailAfterBurner()
  76   {
  77 //destructor
  78
  79     if (fInternalStacks) //delete stacks
  80      {
  81        fInternalStacks->SetOwner();
  82        delete fInternalStacks;
  83     }
  84     if (fAfterBurnerEntries) delete fAfterBurnerEntries; //delete entries
  85   }
  86 /*********************************************************************/
  87 /*********************************************************************/
  88
  89 void AliGenCocktailAfterBurner::
  90 AddAfterBurner(AliGenerator *AfterBurner, char* Name, Float_t RateExp)
  91 {
  92 //
  93 //  Forward parameters to the new AfterBurner
  94
  95     if (gDebug>0)cout<<"AliGenCocktailAfterBurner::AddAfterBurner  Named "<<Name<<endl;
  96
  97     if(TestBit(kPtRange))
  98         AfterBurner->SetPtRange(fPtMin,fPtMax);
  99     if(TestBit(kMomentumRange))
 100         AfterBurner->SetMomentumRange(fPMin,fPMax);
 101
 102     AfterBurner->SetYRange(fYMin,fYMax);
 103     AfterBurner->SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
 104     AfterBurner->SetThetaRange(fThetaMin*180/TMath::Pi(),fThetaMax*180/TMath::Pi());
 105     AfterBurner->SetOrigin(fOrigin[0], fOrigin[1], fOrigin[2]);
 106     AfterBurner->SetSigma(fOsigma[0], fOsigma[1], fOsigma[2]);
 107     AfterBurner->SetVertexSmear(fVertexSmear);
 108     AfterBurner->SetTrackingFlag(fTrackIt);
 109 //
 110 //  Add AfterBurner to list
 111
 112     AliGenCocktailEntry *entry =
 113         new AliGenCocktailEntry(AfterBurner, Name, RateExp);
 114     fAfterBurnerEntries->Add(entry);
 115     fNAfterBurners++;
 116 //
 117
 118 }
 119 /*********************************************************************/
 120 /*********************************************************************/
 121
 122 void AliGenCocktailAfterBurner::Init()
 123 {
 124 // Initialisation
 125     fGenerationDone = kFALSE;
 126     if (fInternalStacks) //delete stacks
 127      {
 128        fInternalStacks->SetOwner();
 129        fInternalStacks->Delete(); //clean after previous generation cycle
 130      }
 131
 132     this->AliGenCocktail::Init();
 133
 134     if (gDebug>0) cout<<"AliGenCocktailAfterBurner::Init"<<endl;
 135     TIter next(fAfterBurnerEntries);
 136     AliGenCocktailEntry *entry;
 137     //
 138     // Loop over generators and initialize
 139     while((entry = (AliGenCocktailEntry*)next())) {
 140         entry->Generator()->Init();
 141     }
 142 }
 143 /*********************************************************************/
 144 /*********************************************************************/
 145
 146 void AliGenCocktailAfterBurner::Generate()
 147 {
 148 //
 149 // Generate event
 150 //  Firsts runs each generator for all events
 151 //  than after burners ones for each event
 152 //
 153 //  It generates and processes all events during
 154 //  first call only.
 155 //  In next calls it just returns already generated
 156 //  and processed events to the gAlice
 157
 158     if (gDebug>0)
 159       cout<<"#####################################"<<endl
 160           <<"#AliGenCocktailAfterBurner::Generate#"<<endl
 161           <<"#####################################"<<endl;
 162
 163     Int_t i; //iterator
 164     AliStack * stack;
 165
 166     if (fGenerationDone)
 167     {//if generation is done (in first call)
 168      //just copy particles from the stack to the gAlice
 169       SetTracks(++fCurrentEvent);
 170       cout<<"Returning event "<<fCurrentEvent<<endl;
 171       return;
 172     }
 173     else
 174     { //Here we are in the first call of the method
 175       fCurrentEvent=0;
 176       Int_t numberOfEvents = gAlice->GetEventsPerRun();
 177       //Create stacks
 178       fInternalStacks = new TObjArray(numberOfEvents + fNBgEvents); //Create array of internal stacks
 179       for(i=0;i<numberOfEvents + fNBgEvents;i++)
 180        {
 181         stack = new AliStack(10000);
 182         stack->Reset();
 183         fInternalStacks->Add(stack);
 184        }
 185 /*********************************************************************/
 186       TIter next(fEntries);
 187       AliGenCocktailEntry *entry;
 188       AliGenCocktailEntry *e1;
 189       AliGenCocktailEntry *e2;
 190       TObjArray *partArray;
 191   //
 192   // Loop over generators and generate events
 193       Int_t igen=0;
 194       while((entry = (AliGenCocktailEntry*)next()))
 195       {
 196         igen++;
 197         cout<<"Generator "<<igen<<"  : "<<entry->GetName()<<endl;
 198 /***********************************************/
 199 //First generator for all evenets, than second for all events, etc...
 200         for(i=0;i<numberOfEvents + fNBgEvents;i++)
 201           {
 202             cout<<"                  EVENT "<<i<<endl;
 203             stack = GetStack(i);
 204             partArray = stack->Particles();
 205             fCurrentGenerator = entry->Generator();
 206             fCurrentGenerator->SetStack(stack);
 207             if (igen ==1)
 208               {
 209                 entry->SetFirst(0);
 210               }
 211             else
 212               {
 213                 entry->SetFirst((partArray->GetEntriesFast())+1);
 214               }
 215                 fCurrentGenerator->Generate();
 216                 entry->SetLast(partArray->GetEntriesFast());
 217            }
 218 /***********************************************/
 219       }
 220       next.Reset();
 221       while((entry = (AliGenCocktailEntry*)next()))
 222         {
 223           entry->PrintInfo();
 224         }
 225       for ( entry=FirstGenerator();entry;entry=NextGenerator() )
 226         {
 227           entry->PrintInfo();
 228         }
 229       for (FirstGeneratorPair(e1,e2); (e1&&e2); NextGeneratorPair(e1,e2) )
 230         {
 231           printf("\n -----------------------------");
 232           e1->PrintInfo();
 233           e2->PrintInfo();
 234         }
 235
 236
 237       /***********************************************/
 238       /*******After Burners Processing****************/
 239       /***********************************************/
 240       TIter nextAfterBurner(fAfterBurnerEntries);
 241       AliGenCocktailEntry *afterBurnerEntry;
 242       Int_t iab =0; //number of current after burner / counter
 243
 244       cout<<"\n\nRunning After Burners"<<endl;
 245       while((afterBurnerEntry = (AliGenCocktailEntry*)nextAfterBurner()))
 246         {
 247           cout<<"After Burner "<<iab++<<"  :"<<afterBurnerEntry->GetName()<<endl;
 248           fCurrentGenerator = afterBurnerEntry->Generator();
 249           fCurrentGenerator->Generate();
 250         }
 251       cout<<endl<<"Finished. Processed "<<iab<<" After Burners"<<endl;
 252
 253       /***********************************************/
 254       /***********************************************/
 255       /***********************************************/
 256
 257       fGenerationDone=kTRUE;
 258       SetTracks(0); //copy event 0 to gAlice stack
 259
 260 /*********************************************************************/
 261
 262     }//else generated
 263 }
 264 /*********************************************************************/
 265 /*********************************************************************/
 266
 267 AliGenCocktailAfterBurner& AliGenCocktailAfterBurner::operator=(const  AliGenCocktailAfterBurner& rhs)
 268 {
 269 // Assignment operator
 270     return *this;
 271 }
 272 /*********************************************************************/
 273 /*********************************************************************/
 274
 275 AliStack* AliGenCocktailAfterBurner::GetStack(Int_t n) const
 276 {
 277 //Returns the pointer to the N'th stack (event)
 278   if( ( n<0 ) || ( n>=GetNumberOfEvents() ) )
 279     {
 280       Fatal("AliGenCocktailAfterBurner::GetStack","Asked for non existing stack (%d)",n);
 281       return 0;
 282     }
 283     return ((AliStack*) fInternalStacks->At(n) );
 284 }
 285 /*********************************************************************/
 286 /*********************************************************************/
 287
 288 void AliGenCocktailAfterBurner::SetActiveEventNumber(Int_t actev)
 289 {
 290 //Set Active Events Number and Active Stack
 291 //There is only one active event number
 292 //Made fo convinience of work with AfterBurners (HBT processor)
 293
 294     fActiveEvent = actev;
 295     fActiveStack = GetStack(actev);
 296 }
 297 /*********************************************************************/
 298 /*********************************************************************/
 299
 300 void AliGenCocktailAfterBurner::SetTracks(Int_t stackno)
 301 {
 302 //Method which copies tracks from given stack to the
 303 //gAlice's stack
 304     AliStack* instack = GetStack(stackno);
 305     Int_t done;
 306     Int_t parent;
 307     Int_t pdg;
 308     Double_t px, py, pz, e, vx, vy, vz, tof, polx, poly, polz;
 309     TMCProcess mech;
 310     Int_t ntr;
 311     Float_t weight;
 312     TVector3 pol;
 313
 314     TParticle * p;
 315     Int_t n = instack->GetNtrack();
 316     if (gDebug)
 317     {
 318       cout<<"AliGenCocktailAfterBurner::SetTracks("<<stackno<<"). Number of particles is: "<<n<<"\n";
 319     }
 320
 321     for(Int_t i = 0; i < n; i++)
 322     {
 323
 324       p = instack->Particle(i);
 325       done = !p->TestBit(kDoneBit);
 326       parent = p->GetMother(0);
 327       pdg = p->GetPdgCode();
 328       px = p->Px();
 329       py = p->Py();
 330       pz = p->Pz();
 331       e  = p->Energy();
 332       vx = p->Vx();
 333       vy = p->Vy();
 334       vz = p->Vz();
 335       tof = p->T();
 336       p->GetPolarisation(pol);
 337       polx = pol.X();
 338       poly = pol.Y();
 339       polz = pol.Z();
 340       mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
 341       weight = p->GetWeight();
 342
 343       gAlice->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
 344                        polx, poly, polz, mech, ntr, weight);
 345     }
 346 }
 347 /*********************************************************************/
 348 /*********************************************************************/
 349
 350 TMCProcess AliGenCocktailAfterBurner::IntToMCProcess(Int_t no)
 351 {
 352  //Mothod used to convert uniqueID (integer) to TMCProcess type
 353     const TMCProcess kMCprocesses[kMaxMCProcess] =
 354     {
 355      kPNoProcess, kPMultipleScattering, kPEnergyLoss, kPMagneticFieldL,
 356      kPDecay, kPPair, kPCompton, kPPhotoelectric, kPBrem, kPDeltaRay,
 357      kPAnnihilation, kPHadronic, kPNoProcess, kPEvaporation, kPNuclearFission,
 358      kPNuclearAbsorption, kPPbarAnnihilation, kPNCapture, kPHElastic,
 359      kPHInhelastic, kPMuonNuclear, kPTOFlimit,kPPhotoFission, kPNoProcess,
 360      kPRayleigh, kPNoProcess, kPNoProcess, kPNoProcess, kPNull, kPStop
 361     };
 362
 363     for (Int_t i = 0;i<kMaxMCProcess;i++)
 364     {
 365       if (kMCprocesses[i] == no)
 366         {
 367           return kMCprocesses[i];
 368         }
 369     }
 370     return kPNoProcess;
 371 }