PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx

   1 #include "TChain.h"\r
   2 #include "TList.h"\r
   3 #include "TCanvas.h"\r
   4 #include "TParticle.h"\r
   5 #include "TLorentzVector.h"\r
   6 #include "TGraphErrors.h"\r
   7 #include "TH1F.h"\r
   8 #include "TH2F.h"\r
   9 #include "TArrayF.h"\r
  10 #include "TF1.h"\r
  11 #include "TRandom.h"\r
  12 #include "TFile.h"\r
  13 \r
  14 #include "AliAnalysisManager.h"\r
  15 #include "AliLog.h"\r
  16 \r
  17 #include "AliAnalysisTaskToyModel.h"\r
  18 #include "AliBalance.h"\r
  19 \r
  20 \r
  21 // Analysis task for the toy model analysis\r
  22 // Authors: Panos.Christakoglou@nikhef.nl\r
  23 \r
  24 using std::cout;\r
  25 using std::endl;\r
  26 \r
  27 ClassImp(AliAnalysisTaskToyModel)\r
  28 \r
  29 //________________________________________________________________________\r
  30 AliAnalysisTaskToyModel::AliAnalysisTaskToyModel() \r
  31 : TObject(),\r
  32   fUseDebug(kFALSE),\r
  33   fBalance(0),\r
  34   fRunShuffling(kFALSE), fShuffledBalance(0),\r
  35   fList(0), fListBF(0), fListBFS(0),\r
  36   fHistEventStats(0),\r
  37   fHistNumberOfAcceptedParticles(0),\r
  38   fHistReactionPlane(0),\r
  39   fHistEtaTotal(0), fHistEta(0),\r
  40   fHistRapidity(0),\r
  41   fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),\r
  42   fHistPhi(0),\r
  43   fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),\r
  44   fHistPt(0),\r
  45   fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),\r
  46   fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),\r
  47   fNetChargeMean(0.0), fNetChargeSigma(0.0),\r
  48   fPtMin(0.0), fPtMax(100.0),\r
  49   fEtaMin(-1.0), fEtaMax(1.0),\r
  50   fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),\r
  51   fUseAllCharges(kFALSE), fParticleMass(0.0),\r
  52   fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),\r
  53   fReactionPlane(0.0),\r
  54   fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0), \r
  55   fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),\r
  56   fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),\r
  57   fPionPercentage(0.8), fPionMass(0.0),\r
  58   fPtSpectraPions(0), fTemperaturePions(100.),\r
  59   fAzimuthalAnglePions(0), fDirectedFlowPions(0.0), \r
  60   fEllipticFlowPions(0.0), fTriangularFlowPions(0.0), \r
  61   fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),\r
  62   fKaonPercentage(0.8), fKaonMass(0.0),\r
  63   fPtSpectraKaons(0), fTemperatureKaons(100.),\r
  64   fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0), \r
  65   fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),\r
  66   fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),\r
  67   fProtonPercentage(0.8), fProtonMass(0.0),\r
  68   fPtSpectraProtons(0), fTemperatureProtons(100.),\r
  69   fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0), \r
  70   fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),\r
  71   fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),\r
  72   fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1) {\r
  73   // Constructor\r
  74 }\r
  75 \r
  76 //________________________________________________________________________\r
  77 AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {\r
  78   //Destructor\r
  79   delete fPtSpectraAllCharges;\r
  80   delete fAzimuthalAngleAllCharges;\r
  81   delete fPtSpectraPions;\r
  82   delete fAzimuthalAnglePions;\r
  83   delete fPtSpectraKaons;\r
  84   delete fAzimuthalAngleKaons;\r
  85   delete fPtSpectraProtons;\r
  86   delete fAzimuthalAngleProtons;\r
  87 }\r
  88 \r
  89 //________________________________________________________________________\r
  90 void AliAnalysisTaskToyModel::Init() {\r
  91   //Initialize objects\r
  92   //==============gRandom Seed=======================//\r
  93   gRandom->SetSeed(0);  //seed is set to a random value (depending on machine clock)\r
  94   //==============gRandom Seed=======================//\r
  95 \r
  96   //==============Particles and spectra==============//\r
  97   TParticle *pion = new TParticle();\r
  98   pion->SetPdgCode(211);\r
  99   fPionMass = pion->GetMass();\r
 100 \r
 101   TParticle *kaon = new TParticle();\r
 102   kaon->SetPdgCode(321);\r
 103   fKaonMass = kaon->GetMass();\r
 104   \r
 105   TParticle *proton = new TParticle();\r
 106   proton->SetPdgCode(2212);\r
 107   fProtonMass = proton->GetMass();\r
 108 \r
 109   if(fUseAllCharges) {\r
 110     fParticleMass = fPionMass;\r
 111     fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 112     fPtSpectraAllCharges->SetParName(0,"Mass");\r
 113     fPtSpectraAllCharges->SetParName(1,"Temperature");\r
 114   }\r
 115   else {\r
 116     fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 117     fPtSpectraPions->SetParName(0,"Mass");\r
 118     fPtSpectraPions->SetParName(1,"Temperature");\r
 119     \r
 120     fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 121     fPtSpectraKaons->SetParName(0,"Mass");\r
 122     fPtSpectraKaons->SetParName(1,"Temperature");\r
 123     \r
 124     fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);\r
 125     fPtSpectraProtons->SetParName(0,"Mass");\r
 126     fPtSpectraProtons->SetParName(1,"Temperature");\r
 127   }\r
 128   //==============Particles and spectra==============//\r
 129 \r
 130   //==============Flow values==============//\r
 131   if(fUseAllCharges) {\r
 132     if(fUseDebug) {\r
 133       Printf("Directed flow: %lf",fDirectedFlowAllCharges);\r
 134       Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);\r
 135       Printf("Triangular flow: %lf",fTriangularFlowAllCharges);\r
 136       Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);\r
 137       Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);\r
 138     }\r
 139 \r
 140     fAzimuthalAngleAllCharges = new TF1("fAzimuthalAngleAllCharges","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
 141     fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");\r
 142     fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");\r
 143     fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow"); \r
 144     fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");\r
 145     fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");\r
 146     fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");\r
 147   }\r
 148   else {\r
 149     fAzimuthalAnglePions = new TF1("fAzimuthalAnglePions","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
 150     fAzimuthalAnglePions->SetParName(0,"Reaction Plane");\r
 151     fAzimuthalAnglePions->SetParName(1,"Directed flow");\r
 152     fAzimuthalAnglePions->SetParName(2,"Elliptic flow"); \r
 153     fAzimuthalAnglePions->SetParName(3,"Triangular flow");\r
 154     fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");\r
 155     fAzimuthalAnglePions->SetParName(5,"Pentangular flow");\r
 156     \r
 157     fAzimuthalAngleKaons = new TF1("fAzimuthalAngleKaons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
 158     fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");\r
 159     fAzimuthalAngleKaons->SetParName(1,"Directed flow");\r
 160     fAzimuthalAngleKaons->SetParName(2,"Elliptic flow"); \r
 161     fAzimuthalAngleKaons->SetParName(3,"Triangular flow");\r
 162     fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");\r
 163     fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");\r
 164     \r
 165     fAzimuthalAngleProtons = new TF1("fAzimuthalAngleProtons","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2*(x-[0]))+2.*[3]*TMath::Cos(3*(x-[0]))+2.*[4]*TMath::Cos(4*(x-[0]))+2.*[5]*TMath::Cos(5*(x-[0]))",0.,2.*TMath::Pi());\r
 166     fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");\r
 167     fAzimuthalAngleProtons->SetParName(1,"Directed flow");\r
 168     fAzimuthalAngleProtons->SetParName(2,"Elliptic flow"); \r
 169     fAzimuthalAngleProtons->SetParName(3,"Triangular flow");\r
 170     fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");\r
 171     fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");\r
 172   }\r
 173   //==============Flow values==============//\r
 174 }\r
 175 \r
 176 //________________________________________________________________________\r
 177 void AliAnalysisTaskToyModel::CreateOutputObjects() {\r
 178   // Create histograms\r
 179   // Called once\r
 180   if(!fBalance) {\r
 181     fBalance = new AliBalance();\r
 182     fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
 183   }\r
 184   if(fRunShuffling) {\r
 185     if(!fShuffledBalance) {\r
 186       fShuffledBalance = new AliBalance();\r
 187       fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
 188     }\r
 189   }\r
 190 \r
 191   //QA list\r
 192   fList = new TList();\r
 193   fList->SetName("listQA");\r
 194   fList->SetOwner();\r
 195 \r
 196   //Balance Function list\r
 197   fListBF = new TList();\r
 198   fListBF->SetName("listBF");\r
 199   fListBF->SetOwner();\r
 200 \r
 201   if(fRunShuffling) {\r
 202     fListBFS = new TList();\r
 203     fListBFS->SetName("listBFShuffled");\r
 204     fListBFS->SetOwner();\r
 205   }\r
 206 \r
 207   //==============QA================//\r
 208   //Event stats.\r
 209   TString gCutName[4] = {"Total","Offline trigger",\r
 210                          "Vertex","Analyzed"};\r
 211   fHistEventStats = new TH1F("fHistEventStats",\r
 212                              "Event statistics;;N_{events}",\r
 213                              4,0.5,4.5);\r
 214   for(Int_t i = 1; i <= 4; i++)\r
 215     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 216   fList->Add(fHistEventStats);\r
 217 \r
 218   fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);\r
 219   fList->Add(fHistNumberOfAcceptedParticles);\r
 220   \r
 221   fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());\r
 222   fList->Add(fHistReactionPlane);\r
 223 \r
 224   //Pseudo-rapidity\r
 225   fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.); \r
 226   fList->Add(fHistEtaTotal);\r
 227 \r
 228   fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5); \r
 229   fList->Add(fHistEta);\r
 230 \r
 231   //Rapidity\r
 232   fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5); \r
 233   fList->Add(fHistRapidity);\r
 234   fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5); \r
 235   fList->Add(fHistRapidityPions);\r
 236   fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5); \r
 237   fList->Add(fHistRapidityKaons);\r
 238   fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5); \r
 239   fList->Add(fHistRapidityProtons);\r
 240 \r
 241   //Phi\r
 242   fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 243   fList->Add(fHistPhi);\r
 244 \r
 245   fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 246   fList->Add(fHistPhiPions);\r
 247   fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 248   fList->Add(fHistPhiKaons);\r
 249   fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());\r
 250   fList->Add(fHistPhiProtons);\r
 251 \r
 252 \r
 253   //Pt\r
 254   fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 255   fList->Add(fHistPt);\r
 256 \r
 257   fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 258   fList->Add(fHistPtPions);\r
 259   fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 260   fList->Add(fHistPtKaons);\r
 261   fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);\r
 262   fList->Add(fHistPtProtons);\r
 263 \r
 264   //==============Balance function histograms================//\r
 265   // Initialize histograms if not done yet\r
 266   if(!fBalance->GetHistNp(0)){\r
 267     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 268     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 269     fBalance->InitHistograms();\r
 270   }\r
 271 \r
 272   if(fRunShuffling) {\r
 273     if(!fShuffledBalance->GetHistNp(0)) {\r
 274       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 275       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 276       fShuffledBalance->InitHistograms();\r
 277     }\r
 278   }\r
 279 \r
 280   for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
 281     fListBF->Add(fBalance->GetHistNp(a));\r
 282     fListBF->Add(fBalance->GetHistNn(a));\r
 283     fListBF->Add(fBalance->GetHistNpn(a));\r
 284     fListBF->Add(fBalance->GetHistNnn(a));\r
 285     fListBF->Add(fBalance->GetHistNpp(a));\r
 286     fListBF->Add(fBalance->GetHistNnp(a));\r
 287 \r
 288     if(fRunShuffling) {\r
 289       fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
 290       fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
 291       fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
 292       fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
 293       fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
 294       fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
 295     }  \r
 296   }\r
 297 \r
 298   // Post output data.\r
 299   //PostData(1, fList);\r
 300   //PostData(2, fListBF);\r
 301   //if(fRunShuffling) PostData(3, fListBFS);\r
 302 }\r
 303 \r
 304 //________________________________________________________________________\r
 305 void AliAnalysisTaskToyModel::Run(Int_t nEvents) {\r
 306   // Main loop\r
 307   // Called for each event\r
 308   Double_t vCharge = 0;\r
 309   Double_t vY = 0.0;\r
 310   Double_t vEta = 0.0;\r
 311   Double_t vPhi = 0.0;\r
 312   Double_t vP[3] = {0.,0.,0.};\r
 313   Double_t vPt = 0.0;\r
 314   Double_t vE = 0.0;\r
 315   Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;\r
 316 \r
 317   if(fUseAllCharges) {\r
 318     fPtSpectraAllCharges->SetParameter(0,fParticleMass);\r
 319     fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);\r
 320 \r
 321     fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);\r
 322     fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);\r
 323     fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);\r
 324     fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);\r
 325     fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);\r
 326   }\r
 327   else {\r
 328     fPtSpectraPions->SetParameter(0,fPionMass);\r
 329     fPtSpectraPions->SetParameter(1,fTemperaturePions);\r
 330     fPtSpectraKaons->SetParameter(0,fKaonMass);\r
 331     fPtSpectraKaons->SetParameter(1,fTemperatureKaons);\r
 332     fPtSpectraProtons->SetParameter(0,fProtonMass);\r
 333     fPtSpectraProtons->SetParameter(1,fTemperatureProtons);\r
 334 \r
 335     fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);\r
 336     fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);\r
 337     fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);\r
 338     fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);\r
 339     fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);\r
 340 \r
 341     fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);\r
 342     fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);\r
 343     fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);\r
 344     fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);\r
 345     fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);\r
 346 \r
 347     fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);\r
 348     fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);\r
 349     fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);\r
 350     fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);\r
 351     fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);\r
 352   }\r
 353 \r
 354   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {\r
 355     // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)\r
 356     vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled\r
 357     vector<Double_t> *chargeVector[9];          // original charge\r
 358     for(Int_t i = 0; i < 9; i++){\r
 359       chargeVectorShuffle[i] = new vector<Double_t>;\r
 360       chargeVector[i]        = new vector<Double_t>;\r
 361     }\r
 362 \r
 363     fHistEventStats->Fill(1);\r
 364     fHistEventStats->Fill(2);\r
 365     fHistEventStats->Fill(3);\r
 366 \r
 367     if((iEvent%1000) == 0) \r
 368       cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;\r
 369     \r
 370     //Multiplicities\r
 371     Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));\r
 372     Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));\r
 373     \r
 374     Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);\r
 375     Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;\r
 376     if(fUseDebug) \r
 377       Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
 378 \r
 379     //Int_t nGeneratedPositive = 0, nGeneratedNegative = 0;\r
 380     Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;\r
 381     \r
 382     //Randomization of the reaction plane\r
 383     //fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();\r
 384     fReactionPlane = 0.0;\r
 385     if(fUseAllCharges) \r
 386       fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);\r
 387     else {\r
 388       fAzimuthalAnglePions->SetParameter(0,fReactionPlane);\r
 389       fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);\r
 390       fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);\r
 391     }\r
 392     \r
 393     Int_t gNumberOfAcceptedParticles = 0;\r
 394     Int_t gNumberOfAcceptedPositiveParticles = 0;\r
 395     Int_t gNumberOfAcceptedNegativeParticles = 0;\r
 396  \r
 397     //Generate positive particles\r
 398     for(Int_t iParticleCount = 0; iParticleCount < nGeneratedPositive; iParticleCount++) {\r
 399       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
 400       if(fUseDebug) \r
 401         Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);\r
 402 \r
 403       //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
 404       vEta = gRandom->Gaus(0.0,4.0);\r
 405 \r
 406       //Fill QA histograms (full phase space)\r
 407       fHistEtaTotal->Fill(vEta);\r
 408 \r
 409       vCharge = 1.0;\r
 410       //nGeneratedPositive += 1;\r
 411       \r
 412       //Acceptance\r
 413       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
 414 \r
 415       if(!fUseAllCharges) {\r
 416         //Decide the specie\r
 417         Double_t randomNumberSpecies = gRandom->Rndm();\r
 418         if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
 419           nGeneratedPions += 1;\r
 420           vPt = fPtSpectraPions->GetRandom();\r
 421           vPhi = fAzimuthalAnglePions->GetRandom();\r
 422           fParticleMass = fPionMass;\r
 423           isPion = kTRUE;\r
 424         }\r
 425         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
 426           nGeneratedKaons += 1;\r
 427           vPt = fPtSpectraKaons->GetRandom();\r
 428           vPhi = fAzimuthalAngleKaons->GetRandom();\r
 429           fParticleMass = fKaonMass;\r
 430           isKaon = kTRUE;\r
 431         }\r
 432         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
 433           nGeneratedProtons += 1;\r
 434           vPt = fPtSpectraProtons->GetRandom();\r
 435           vPhi = fAzimuthalAngleProtons->GetRandom();\r
 436           fParticleMass = fProtonMass;\r
 437           isProton = kTRUE;\r
 438         }\r
 439       }\r
 440       else {\r
 441         vPt = fPtSpectraAllCharges->GetRandom();\r
 442         vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
 443       }\r
 444       \r
 445       vP[0] = vPt*TMath::Cos(vPhi);\r
 446       vP[1] = vPt*TMath::Sin(vPhi);\r
 447       vP[2] = vPt*TMath::SinH(vEta);\r
 448       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 449                         TMath::Power(vP[0],2) +\r
 450                         TMath::Power(vP[1],2) +\r
 451                         TMath::Power(vP[2],2));\r
 452       \r
 453       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
 454       \r
 455       //pt coverage\r
 456       if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
 457       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 458 \r
 459       //acceptance filter\r
 460       if(fUseAcceptanceParameterization) {\r
 461         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 462         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
 463           continue;\r
 464       }\r
 465 \r
 466       gNumberOfAcceptedPositiveParticles += 1;\r
 467 \r
 468       //Fill QA histograms (acceptance);\r
 469       fHistEta->Fill(vEta);\r
 470       fHistRapidity->Fill(vY);\r
 471       fHistPhi->Fill(vPhi);\r
 472       fHistPt->Fill(vPt);\r
 473       if(isPion) {\r
 474         fHistRapidityPions->Fill(vY);\r
 475         fHistPhiPions->Fill(vPhi);\r
 476         fHistPtPions->Fill(vPt);\r
 477       }\r
 478       else if(isKaon) {\r
 479         fHistRapidityKaons->Fill(vY);\r
 480         fHistPhiKaons->Fill(vPhi);\r
 481         fHistPtKaons->Fill(vPt);\r
 482       }\r
 483       else if(isProton) {\r
 484         fHistRapidityProtons->Fill(vY);\r
 485         fHistPhiProtons->Fill(vPhi);\r
 486         fHistPtProtons->Fill(vPt);\r
 487       }\r
 488 \r
 489       // fill charge vector\r
 490       chargeVector[0]->push_back(vCharge);\r
 491       chargeVector[1]->push_back(vY);\r
 492       chargeVector[2]->push_back(vEta);\r
 493       chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
 494       chargeVector[4]->push_back(vP[0]);\r
 495       chargeVector[5]->push_back(vP[1]);\r
 496       chargeVector[6]->push_back(vP[2]);\r
 497       chargeVector[7]->push_back(vPt);\r
 498       chargeVector[8]->push_back(vE);\r
 499       \r
 500       if(fRunShuffling) {\r
 501         chargeVectorShuffle[0]->push_back(vCharge);\r
 502         chargeVectorShuffle[1]->push_back(vY);\r
 503         chargeVectorShuffle[2]->push_back(vEta);\r
 504         chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
 505         chargeVectorShuffle[4]->push_back(vP[0]);\r
 506         chargeVectorShuffle[5]->push_back(vP[1]);\r
 507         chargeVectorShuffle[6]->push_back(vP[2]);\r
 508         chargeVectorShuffle[7]->push_back(vPt);\r
 509         chargeVectorShuffle[8]->push_back(vE);\r
 510       }\r
 511       gNumberOfAcceptedParticles += 1;\r
 512     }//generated positive particle loop\r
 513  \r
 514     //Generate negative particles\r
 515     for(Int_t iParticleCount = 0; iParticleCount < nGeneratedNegative; iParticleCount++) {\r
 516       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
 517       if(fUseDebug) \r
 518         Printf("Generating negative: %d(%d)",iParticleCount+1,nGeneratedNegative);\r
 519 \r
 520       //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
 521       vEta = gRandom->Gaus(0.0,4.0);\r
 522 \r
 523       //Fill QA histograms (full phase space)\r
 524       fHistEtaTotal->Fill(vEta);\r
 525 \r
 526       vCharge = -1.0;\r
 527       //nGeneratedNegative += 1;\r
 528       \r
 529       //Acceptance\r
 530       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
 531 \r
 532       if(!fUseAllCharges) {\r
 533         //Decide the specie\r
 534         Double_t randomNumberSpecies = gRandom->Rndm();\r
 535         if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
 536           nGeneratedPions += 1;\r
 537           vPt = fPtSpectraPions->GetRandom();\r
 538           vPhi = fAzimuthalAnglePions->GetRandom();\r
 539           fParticleMass = fPionMass;\r
 540           isPion = kTRUE;\r
 541         }\r
 542         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
 543           nGeneratedKaons += 1;\r
 544           vPt = fPtSpectraKaons->GetRandom();\r
 545           vPhi = fAzimuthalAngleKaons->GetRandom();\r
 546           fParticleMass = fKaonMass;\r
 547           isKaon = kTRUE;\r
 548         }\r
 549         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
 550           nGeneratedProtons += 1;\r
 551           vPt = fPtSpectraProtons->GetRandom();\r
 552           vPhi = fAzimuthalAngleProtons->GetRandom();\r
 553           fParticleMass = fProtonMass;\r
 554           isProton = kTRUE;\r
 555         }\r
 556       }\r
 557       else {\r
 558         vPt = fPtSpectraAllCharges->GetRandom();\r
 559         vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
 560       }\r
 561       \r
 562       vP[0] = vPt*TMath::Cos(vPhi);\r
 563       vP[1] = vPt*TMath::Sin(vPhi);\r
 564       vP[2] = vPt*TMath::SinH(vEta);\r
 565       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 566                         TMath::Power(vP[0],2) +\r
 567                         TMath::Power(vP[1],2) +\r
 568                         TMath::Power(vP[2],2));\r
 569       \r
 570       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
 571       \r
 572       //pt coverage\r
 573       if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
 574       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);\r
 575 \r
 576      //acceptance filter\r
 577       if(fUseAcceptanceParameterization) {\r
 578         Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 579         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
 580           continue;\r
 581       }\r
 582 \r
 583       gNumberOfAcceptedNegativeParticles += 1;\r
 584 \r
 585       //Fill QA histograms (acceptance);\r
 586       fHistEta->Fill(vEta);\r
 587       fHistRapidity->Fill(vY);\r
 588       fHistPhi->Fill(vPhi);\r
 589       fHistPt->Fill(vPt);\r
 590       if(isPion) {\r
 591         fHistRapidityPions->Fill(vY);\r
 592         fHistPhiPions->Fill(vPhi);\r
 593         fHistPtPions->Fill(vPt);\r
 594       }\r
 595       else if(isKaon) {\r
 596         fHistRapidityKaons->Fill(vY);\r
 597         fHistPhiKaons->Fill(vPhi);\r
 598         fHistPtKaons->Fill(vPt);\r
 599       }\r
 600       else if(isProton) {\r
 601         fHistRapidityProtons->Fill(vY);\r
 602         fHistPhiProtons->Fill(vPhi);\r
 603         fHistPtProtons->Fill(vPt);\r
 604       }\r
 605 \r
 606       // fill charge vector\r
 607       chargeVector[0]->push_back(vCharge);\r
 608       chargeVector[1]->push_back(vY);\r
 609       chargeVector[2]->push_back(vEta);\r
 610       chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
 611       chargeVector[4]->push_back(vP[0]);\r
 612       chargeVector[5]->push_back(vP[1]);\r
 613       chargeVector[6]->push_back(vP[2]);\r
 614       chargeVector[7]->push_back(vPt);\r
 615       chargeVector[8]->push_back(vE);\r
 616       \r
 617       if(fRunShuffling) {\r
 618         chargeVectorShuffle[0]->push_back(vCharge);\r
 619         chargeVectorShuffle[1]->push_back(vY);\r
 620         chargeVectorShuffle[2]->push_back(vEta);\r
 621         chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
 622         chargeVectorShuffle[4]->push_back(vP[0]);\r
 623         chargeVectorShuffle[5]->push_back(vP[1]);\r
 624         chargeVectorShuffle[6]->push_back(vP[2]);\r
 625         chargeVectorShuffle[7]->push_back(vPt);\r
 626         chargeVectorShuffle[8]->push_back(vE);\r
 627       }\r
 628       gNumberOfAcceptedParticles += 1;\r
 629     }//generated negative particle loop\r
 630     \r
 631     //Dynamical correlations\r
 632     Double_t vChargePrime = 0;\r
 633     Double_t vYPrime = 0.0;\r
 634     Double_t vEtaPrime = 0.0;\r
 635     Double_t vPhiPrime = 0.0;\r
 636     Double_t vPPrime[3] = {0.,0.,0.};\r
 637     Double_t vPtPrime = 0.0;\r
 638     Double_t vEPrime = 0.0;\r
 639     Int_t nGeneratedPositiveDynamicalCorrelations = 0;\r
 640     Int_t nGeneratedNegativeDynamicalCorrelations = 0;\r
 641     //Generate "correlated" particles \r
 642     if(fUseDynamicalCorrelations) {\r
 643       Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);\r
 644       for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {\r
 645         isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;\r
 646         \r
 647         //Pseudo-rapidity sampled from a Gaussian centered @ 0\r
 648         vEta = gRandom->Gaus(0.0,0.1);\r
 649         vCharge = 1.0;\r
 650         nGeneratedPositiveDynamicalCorrelations += 1;\r
 651         \r
 652         vEtaPrime = -vEta;\r
 653         vChargePrime = -1.0;\r
 654         nGeneratedNegativeDynamicalCorrelations += 1;\r
 655           \r
 656         //Acceptance\r
 657         if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;\r
 658         if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;\r
 659       \r
 660         if(!fUseAllCharges) {\r
 661           //Decide the specie\r
 662           Double_t randomNumberSpecies = gRandom->Rndm();\r
 663           if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {\r
 664             nGeneratedPions += 1;\r
 665             vPt = fPtSpectraPions->GetRandom();\r
 666             vPhi = fAzimuthalAnglePions->GetRandom();\r
 667             fParticleMass = fPionMass;\r
 668             isPion = kTRUE;\r
 669           }\r
 670           else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {\r
 671             nGeneratedKaons += 1;\r
 672             vPt = fPtSpectraKaons->GetRandom();\r
 673             vPhi = fAzimuthalAngleKaons->GetRandom();\r
 674             fParticleMass = fKaonMass;\r
 675             isKaon = kTRUE;\r
 676           }\r
 677           else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {\r
 678             nGeneratedProtons += 1;\r
 679             vPt = fPtSpectraProtons->GetRandom();\r
 680             vPtPrime = vPt;\r
 681             vPhi = fAzimuthalAngleProtons->GetRandom();\r
 682             fParticleMass = fProtonMass;\r
 683             isProton = kTRUE;\r
 684           }\r
 685         }\r
 686         else {\r
 687           vPt = fPtSpectraAllCharges->GetRandom();\r
 688           vPhi = fAzimuthalAngleAllCharges->GetRandom();\r
 689         }\r
 690         vPtPrime = vPt;\r
 691         vPhiPrime = vPhi;\r
 692 \r
 693         vP[0] = vPt*TMath::Cos(vPhi);\r
 694         vP[1] = vPt*TMath::Sin(vPhi);\r
 695         vP[2] = vPt*TMath::SinH(vEta);\r
 696         vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 697                           TMath::Power(vP[0],2) +\r
 698                           TMath::Power(vP[1],2) +\r
 699                           TMath::Power(vP[2],2));\r
 700         \r
 701         vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));\r
 702 \r
 703         vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);\r
 704         vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);\r
 705         vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);\r
 706         vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +\r
 707                           TMath::Power(vPPrime[0],2) +\r
 708                           TMath::Power(vPPrime[1],2) +\r
 709                           TMath::Power(vPPrime[2],2));\r
 710         \r
 711         vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));\r
 712       \r
 713         //pt coverage\r
 714         if((vPt < fPtMin) || (vPt > fPtMax)) continue;\r
 715         if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;\r
 716 \r
 717         //acceptance filter\r
 718         if(fUseAcceptanceParameterization) {\r
 719           Double_t gRandomNumberForAcceptance = gRandom->Rndm();\r
 720           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt)) \r
 721             continue;\r
 722           \r
 723           Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();\r
 724           if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime)) \r
 725             continue;\r
 726         }\r
 727       \r
 728         // fill charge vector (positive)\r
 729         chargeVector[0]->push_back(vCharge);\r
 730         chargeVector[1]->push_back(vY);\r
 731         chargeVector[2]->push_back(vEta);\r
 732         chargeVector[3]->push_back(TMath::RadToDeg()*vPhi);\r
 733         chargeVector[4]->push_back(vP[0]);\r
 734         chargeVector[5]->push_back(vP[1]);\r
 735         chargeVector[6]->push_back(vP[2]);\r
 736         chargeVector[7]->push_back(vPt);\r
 737         chargeVector[8]->push_back(vE);\r
 738         \r
 739         if(fRunShuffling) {\r
 740           chargeVectorShuffle[0]->push_back(vCharge);\r
 741           chargeVectorShuffle[1]->push_back(vY);\r
 742           chargeVectorShuffle[2]->push_back(vEta);\r
 743           chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhi);\r
 744           chargeVectorShuffle[4]->push_back(vP[0]);\r
 745           chargeVectorShuffle[5]->push_back(vP[1]);\r
 746           chargeVectorShuffle[6]->push_back(vP[2]);\r
 747           chargeVectorShuffle[7]->push_back(vPt);\r
 748           chargeVectorShuffle[8]->push_back(vE);\r
 749         }\r
 750 \r
 751         // fill charge vector (negative)\r
 752         chargeVector[0]->push_back(vChargePrime);\r
 753         chargeVector[1]->push_back(vYPrime);\r
 754         chargeVector[2]->push_back(vEtaPrime);\r
 755         chargeVector[3]->push_back(TMath::RadToDeg()*vPhiPrime);\r
 756         chargeVector[4]->push_back(vPPrime[0]);\r
 757         chargeVector[5]->push_back(vPPrime[1]);\r
 758         chargeVector[6]->push_back(vPPrime[2]);\r
 759         chargeVector[7]->push_back(vPtPrime);\r
 760         chargeVector[8]->push_back(vEPrime);\r
 761         \r
 762         if(fRunShuffling) {\r
 763           chargeVectorShuffle[0]->push_back(vChargePrime);\r
 764           chargeVectorShuffle[1]->push_back(vYPrime);\r
 765           chargeVectorShuffle[2]->push_back(vEtaPrime);\r
 766           chargeVectorShuffle[3]->push_back(TMath::RadToDeg()*vPhiPrime);\r
 767           chargeVectorShuffle[4]->push_back(vPPrime[0]);\r
 768           chargeVectorShuffle[5]->push_back(vPPrime[1]);\r
 769           chargeVectorShuffle[6]->push_back(vPPrime[2]);\r
 770           chargeVectorShuffle[7]->push_back(vPtPrime);\r
 771           chargeVectorShuffle[8]->push_back(vEPrime);\r
 772         }\r
 773 \r
 774         gNumberOfAcceptedParticles += 2;\r
 775       }//loop over the dynamical correlations\r
 776     }// usage of dynamical correlations\r
 777 \r
 778     if(fUseDebug) {\r
 779       Printf("=======================================================");\r
 780       Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);\r
 781       Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);\r
 782       Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);\r
 783       Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);\r
 784       if(!fUseAllCharges)\r
 785         Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);\r
 786       //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());\r
 787     }\r
 788 \r
 789     fHistEventStats->Fill(4);\r
 790     fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);\r
 791     fHistReactionPlane->Fill(fReactionPlane);\r
 792 \r
 793     //Calculate the balance function\r
 794     fBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVector);\r
 795     \r
 796     if(fRunShuffling) {\r
 797       // shuffle charges\r
 798       random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );\r
 799       fShuffledBalance->CalculateBalance(gNumberOfAcceptedParticles,chargeVectorShuffle);\r
 800     }\r
 801   }\r
 802 }      \r
 803 \r
 804 //________________________________________________________________________\r
 805 void  AliAnalysisTaskToyModel::FinishOutput() {\r
 806   //Printf("END BF");\r
 807 \r
 808   TFile *gOutput = TFile::Open("outputToyModel.root","recreate");\r
 809   fList->Write();\r
 810 \r
 811   if (!fBalance) {\r
 812     Printf("ERROR: fBalance not available");\r
 813     return;\r
 814   }  \r
 815   fListBF->Write();\r
 816 \r
 817   if(fRunShuffling) {\r
 818     if (!fShuffledBalance) {\r
 819       Printf("ERROR: fShuffledBalance not available");\r
 820       return;\r
 821     }\r
 822     fListBFS->Write();\r
 823   }\r
 824   gOutput->Close();\r
 825 }\r
 826 \r