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        }
 220 /*********************************************************************/
 221       TIter next(fEntries);
 222       AliGenCocktailEntry *entry;
 223       AliGenCocktailEntry *e1;
 224       AliGenCocktailEntry *e2;
 225       const TObjArray *partArray;
 226   //
 227   // Loop over generators and generate events
 228       Int_t igen=0;
 229       while((entry = (AliGenCocktailEntry*)next()))
 230       {
 231         igen++;
 232         cout<<"Generator "<<igen<<"  : "<<entry->GetName()<<endl;
 233 /***********************************************/
 234 //First generator for all evenets, than second for all events, etc...
 235         for(i = 0; i < numberOfEvents + fNBgEvents; i++)
 236         {
 237             cout<<"                  EVENT "<<i << endl;
 238             stack = GetStack(i);
 239             partArray = stack->Particles();
 240             fCurrentGenerator = entry->Generator();
 241             fCurrentGenerator->SetStack(stack);
 242             if (igen ==1)
 243             {
 244                 entry->SetFirst(0);
 245             }
 246             else
 247             {
 248                 entry->SetFirst((partArray->GetEntriesFast())+1);
 249             }
 250             // Set the vertex for the generator
 251             Int_t ioff = 3 * i;
 252             fCurrentGenerator->SetVertex(eventVertex.At(ioff), eventVertex.At(ioff + 1), eventVertex.At(ioff + 2));
 253             fHeader = fHeaders[i];
 254             // Set the vertex for the cocktail
 255             TArrayF v(3);
 256             for (Int_t j=0; j<3; j++) v[j] = eventVertex.At(ioff + j);
 257             fHeader->SetPrimaryVertex(v);
 258             // Generate event
 259             fCurrentGenerator->Generate();
 260             //
 261             entry->SetLast(partArray->GetEntriesFast());
 262
 263             if (fCurrentGenerator->ProvidesCollisionGeometry())
 264               fCollisionGeometries[i] =
 265                 new AliCollisionGeometry(*(fCurrentGenerator->CollisionGeometry()));
 266         } // event loop
 267 /***********************************************/
 268       } // generator loop
 269       next.Reset();
 270       while((entry = (AliGenCocktailEntry*)next()))
 271         {
 272           entry->PrintInfo();
 273         }
 274       for ( entry=FirstGenerator();entry;entry=NextGenerator() )
 275         {
 276           entry->PrintInfo();
 277         }
 278       for (FirstGeneratorPair(e1,e2); (e1&&e2); NextGeneratorPair(e1,e2) )
 279         {
 280           printf("\n -----------------------------");
 281           e1->PrintInfo();
 282           e2->PrintInfo();
 283         }
 284
 285
 286       /***********************************************/
 287       /*******After Burners Processing****************/
 288       /***********************************************/
 289       TIter nextAfterBurner(fAfterBurnerEntries);
 290       AliGenCocktailEntry *afterBurnerEntry;
 291       Int_t iab =0; //number of current after burner / counter
 292
 293       cout<<"\n\nRunning After Burners"<<endl;
 294       while((afterBurnerEntry = (AliGenCocktailEntry*)nextAfterBurner()))
 295         {
 296           cout<<"After Burner "<<iab++<<"  :"<<afterBurnerEntry->GetName()<<endl;
 297           fCurrentGenerator = afterBurnerEntry->Generator();
 298           fCurrentGenerator->Generate();
 299         }
 300       cout<<endl<<"Finished. Processed "<<iab<<" After Burners"<<endl;
 301
 302       /***********************************************/
 303       /***********************************************/
 304       /***********************************************/
 305
 306       fGenerationDone=kTRUE;
 307       SetTracks(0); //copy event 0 to gAlice stack
 308
 309 /*********************************************************************/
 310       // Pass the header to gAlice
 311       fHeader = fHeaders[0];
 312       gAlice->SetGenEventHeader(fHeader);
 313     } //else generated
 314 }
 315 /*********************************************************************/
 316 /*********************************************************************/
 317
 318 AliStack* AliGenCocktailAfterBurner::GetStack(Int_t n) const
 319 {
 320 //Returns the pointer to the N'th stack (event)
 321   if( ( n<0 ) || ( n >= (GetNumberOfEvents()) ) )
 322     {
 323       Fatal("AliGenCocktailAfterBurner::GetStack","Asked for non existing stack (%d)",n);
 324       return 0;
 325     }
 326     return ((AliStack*) fInternalStacks->At(n) );
 327 }
 328
 329 /*********************************************************************/
 330 /*********************************************************************/
 331
 332 AliCollisionGeometry* AliGenCocktailAfterBurner::GetCollisionGeometry(Int_t n) const
 333 {
 334 //Returns the pointer to the N'th stack (event)
 335   if( ( n<0 ) || ( n>=GetNumberOfEvents() ) )
 336     {
 337       Fatal("AliGenCocktailAfterBurner::GetCollisionGeometry","Asked for non existing stack (%d)",n);
 338       return 0;
 339     }
 340     return fCollisionGeometries[n];
 341 }
 342
 343 /*********************************************************************/
 344 /*********************************************************************/
 345
 346 void AliGenCocktailAfterBurner::SetActiveEventNumber(Int_t actev)
 347 {
 348 //Set Active Events Number and Active Stack
 349 //There is only one active event number
 350 //Made fo convinience of work with AfterBurners (HBT processor)
 351
 352     fActiveEvent = actev;
 353     fActiveStack = GetStack(actev);
 354 }
 355 /*********************************************************************/
 356 /*********************************************************************/
 357
 358 void AliGenCocktailAfterBurner::SetTracks(Int_t stackno)
 359 {
 360 //Method which copies tracks from given stack to the
 361 //gAlice's stack
 362     AliStack* instack = GetStack(stackno);
 363     Int_t done;
 364     Int_t parent;
 365     Int_t pdg;
 366     Double_t px, py, pz, e, vx, vy, vz, tof, polx, poly, polz;
 367     TMCProcess mech;
 368     Int_t ntr;
 369     Float_t weight;
 370     TVector3 pol;
 371
 372     TParticle * p;
 373     Int_t n = instack->GetNtrack();
 374     if (gDebug)
 375     {
 376       cout<<"AliGenCocktailAfterBurner::SetTracks("<<stackno<<"). Number of particles is: "<<n<<"\n";
 377     }
 378
 379     for(Int_t i = 0; i < n; i++)
 380     {
 381
 382       p = instack->Particle(i);
 383       done = !p->TestBit(kDoneBit);
 384       parent = p->GetMother(0);
 385       pdg = p->GetPdgCode();
 386       px = p->Px();
 387       py = p->Py();
 388       pz = p->Pz();
 389       e  = p->Energy();
 390       vx = p->Vx();
 391       vy = p->Vy();
 392       vz = p->Vz();
 393       tof = p->T();
 394       p->GetPolarisation(pol);
 395       polx = pol.X();
 396       poly = pol.Y();
 397       polz = pol.Z();
 398       mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
 399       weight = p->GetWeight();
 400
 401       gAlice->GetMCApp()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
 402
 403       SetHighWaterMark(ntr) ;
 404
 405     }
 406 }
 407 /*********************************************************************/
 408 /*********************************************************************/
 409
 410 TMCProcess AliGenCocktailAfterBurner::IntToMCProcess(Int_t no)
 411 {
 412  //Mothod used to convert uniqueID (integer) to TMCProcess type
 413     const TMCProcess kMCprocesses[kMaxMCProcess] =
 414     {
 415      kPNoProcess, kPMultipleScattering, kPEnergyLoss, kPMagneticFieldL,
 416      kPDecay, kPPair, kPCompton, kPPhotoelectric, kPBrem, kPDeltaRay,
 417      kPAnnihilation, kPHadronic, kPNoProcess, kPEvaporation, kPNuclearFission,
 418      kPNuclearAbsorption, kPPbarAnnihilation, kPNCapture, kPHElastic,
 419      kPHInhelastic, kPMuonNuclear, kPTOFlimit,kPPhotoFission, kPNoProcess,
 420      kPRayleigh, kPNoProcess, kPNoProcess, kPNoProcess, kPNull, kPStop
 421     };
 422
 423     for (Int_t i = 0;i<kMaxMCProcess;i++)
 424     {
 425       if (kMCprocesses[i] == no)
 426         {
 427           return kMCprocesses[i];
 428         }
 429     }
 430     return kPNoProcess;
 431 }
 432