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