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