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