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