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