]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FASTSIM/AliFastGlauber.cxx
Now it is possible to select kPyCharmppMNRwmi
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
index ceb05e18da41d98b50715b99f0f72a0ded5f3ef4..020b456d89202faf8f860ce426744918d998132c 100644 (file)
@@ -16,7 +16,7 @@
 /* $Id$ */
 //
 // Utility class to make simple Glauber type calculations 
-//           for SYMMTRIC collision geometries (AA):
+//           for SYMMETRIC collision geometries (AA):
 // Impact parameter, production points, reaction plane dependence
 //
 // The SimulateTrigger method can be used for simple MB and hard-process
 // from AliRoot
 #include "AliFastGlauber.h"
 // from root
-#include <TStyle.h>
-#include <TH1F.h>
-#include <TH2F.h>
+#include <Riostream.h>
+#include <TCanvas.h>
 #include <TF1.h>
 #include <TF2.h>
-#include <TLegend.h>
-#include <TCanvas.h>
-#include <TRandom.h>
 #include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TLegend.h>
+#include <TMath.h>
 #include <TROOT.h>
-#include <Riostream.h>
+#include <TRandom.h>
+#include <TStyle.h>
 
 ClassImp(AliFastGlauber)
 
@@ -82,17 +83,27 @@ TF1*    AliFastGlauber::fgWIntRadius     = NULL;
 TF2*    AliFastGlauber::fgWKParticipants = NULL; 
 TF1*    AliFastGlauber::fgWParticipants  = NULL; 
 TF2*    AliFastGlauber::fgWAlmondCurrent = NULL;    
-TF2     AliFastGlauber::fgWAlmondFixedB[40]; 
+TF2*    AliFastGlauber::fgWAlmondFixedB[40]; 
 const Int_t AliFastGlauber::fgkMCInts = 100000;
 Int_t AliFastGlauber::fgCounter = 0;       
 
-AliFastGlauber::AliFastGlauber() : fName()
-{
-  //  Default Constructor
-  //
+AliFastGlauber::AliFastGlauber(): 
+    fWSr0(0.),
+    fWSd(0.), 
+    fWSw(0.), 
+    fWSn(0.), 
+    fSigmaHard(0.),
+    fSigmaNN(0.),  
+    fA(0),         
+    fBmin(0.),     
+    fBmax(0.),     
+    fEllDef(0),    
+    fName()     
+{
+  //  Default Constructor 
   fgCounter++;
   if(fgCounter>1)
-    Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter);
+    Error("AliFastGlauber","More than one instance (%d) is not supported, check your code!",fgCounter);
 
   //  Defaults for Pb
   SetMaxImpact();
@@ -100,10 +111,29 @@ AliFastGlauber::AliFastGlauber() : fName()
   SetPbPbLHC();
 }
 
+AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
+    :TObject(gl),
+     fWSr0(0.),
+     fWSd(0.), 
+     fWSw(0.), 
+     fWSn(0.), 
+     fSigmaHard(0.),
+     fSigmaNN(0.),  
+     fA(0),         
+     fBmin(0.),     
+     fBmax(0.),     
+     fEllDef(0),    
+     fName()     
+{
+// Copy constructor
+    gl.Copy(*this);
+}
+
 AliFastGlauber::~AliFastGlauber()
 {
+// Destructor
   fgCounter--;
-  //if(fgCounter==0) Reset();
+  for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
 }
 
 void AliFastGlauber::SetAuAuRhic()
@@ -148,7 +178,7 @@ void AliFastGlauber::Init(Int_t mode)
   fgWSb->SetParameter(2, fWSw);
   fgWSb->SetParameter(3, fWSn);
 
-  fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
+  fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
   fgWSbz->SetParameter(0, fWSr0);
   fgWSbz->SetParameter(1, fWSd);
   fgWSbz->SetParameter(2, fWSw);
@@ -247,7 +277,7 @@ void AliFastGlauber::Init(Int_t mode)
     for(Int_t k=0; k<40; k++) {
       sprintf(almondName,"WAlmondFixedB%d",k);
       fgWAlmondCurrent = (TF2*)ff->Get(almondName);
-      new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
+      fgWAlmondFixedB[k] = fgWAlmondCurrent;
     }
     delete ff;
   }
@@ -268,7 +298,7 @@ void AliFastGlauber::Init(Int_t mode)
   fgWPathLength->SetParameter(2, 0);     //Pathlength definition
 }
 
-void AliFastGlauber::Reset()
+void AliFastGlauber::Reset() const
 {
   //
   // Reset dynamic allocated formulas
@@ -331,7 +361,7 @@ void AliFastGlauber::DrawOverlap() const
   TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
   c2->cd();
   Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
-  TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0,max);
+  TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
   h2f->SetStats(0);
   h2f->GetXaxis()->SetTitle("b [fm]");
   h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
@@ -681,8 +711,11 @@ Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
   //  MC Integration
   //
   Double_t y = 0;
+  
+
   for (Int_t i = 0; i < fgkMCInts; i++)
     {
+       
       const Double_t kphi = TMath::Pi() * gRandom->Rndm();
       const Double_t kb1  = fgBMax      * gRandom->Rndm();     
       y += fgWStarfi->Eval(kb1, kphi);
@@ -884,11 +917,11 @@ Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
   } // radial steps
 
   Double_t y=0.;
-  if (!kiopt)  // My length definition (is exact for hard disk)
-    if(w) y= 2. * rw / w; 
-  else {
-    const Double_t knorm=fgWSta->Eval(1e-4);
-    if(knorm) y =  TMath::Sqrt(2. * rw * kDr / knorm / knorm);
+  if (!kiopt) { // My length definition (is exact for hard disk)
+      if(w) y= 2. * rw / w; 
+  else {
+      const Double_t knorm=fgWSta->Eval(1e-4);
+      if(knorm) y =  TMath::Sqrt(2. * rw * kDr / knorm / knorm);
   }
   return y; //fm
 }
@@ -999,6 +1032,44 @@ Double_t AliFastGlauber::Binaries(Double_t b) const
   return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
 }
 
+Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
+{
+//
+// Calculate the mean overlap for impact parameter range b1 .. b2
+//
+    Double_t sum  = 0.;
+    Double_t sumc = 0.;
+    Double_t b    = b1;
+    
+    while (b < b2-0.005) {
+       Double_t  nc = GetNumberOfCollisions(b);
+       sum  += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
+       sumc += 10. * fgWSgeo->Eval(b) * 0.01;
+       b += 0.01;
+    }
+    return (sum / CrossSection(b1, b2));
+}
+
+
+Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
+{
+//
+// Calculate the mean number of collisions per event for impact parameter range b1 .. b2
+//
+    Double_t sum  = 0.;
+    Double_t sumc = 0.;
+    Double_t b    = b1;
+    
+    while (b < b2-0.005) {
+       Double_t  nc = GetNumberOfCollisions(b);
+       sum  += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
+       sumc += 10. * fgWSgeo->Eval(b) * 0.01;
+       b += 0.01;
+    }
+    return (sum / CrossSection(b1, b2));
+}
+
+
 Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
 {
   //
@@ -1217,7 +1288,7 @@ void AliFastGlauber::GetRandomBHard(Double_t& b)
   b = fgWSbinary->GetRandom();
   Int_t bin = 2*(Int_t)b;
   if( (b-(Int_t)b) > 0.5) bin++;
-  fgWAlmondCurrent = &fgWAlmondFixedB[bin]; 
+  fgWAlmondCurrent = fgWAlmondFixedB[bin]; 
   return;
 }
 
@@ -1470,13 +1541,13 @@ void AliFastGlauber::PlotAlmonds() const
   gStyle->SetPalette(1,0);
   c->Divide(2,2);
   c->cd(1);
-  fgWAlmondFixedB[0].Draw("cont1");
+  fgWAlmondFixedB[0]->Draw("cont1");
   c->cd(2);
-  fgWAlmondFixedB[10].Draw("cont1");
+  fgWAlmondFixedB[10]->Draw("cont1");
   c->cd(3);
-  fgWAlmondFixedB[20].Draw("cont1");
+  fgWAlmondFixedB[20]->Draw("cont1");
   c->cd(4);
-  fgWAlmondFixedB[30].Draw("cont1");
+  fgWAlmondFixedB[30]->Draw("cont1");
   return;
 }
 
@@ -1763,3 +1834,18 @@ void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
   return;
 }
 
+AliFastGlauber& AliFastGlauber::operator=(const  AliFastGlauber& rhs)
+{
+// Assignment operator
+    rhs.Copy(*this);
+    return *this;
+}
+
+void AliFastGlauber::Copy(TObject&) const
+{
+    //
+    // Copy 
+    //
+    Fatal("Copy","Not implemented!\n");
+}
+