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