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