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