PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskToyModel.cxx

   1 #include "TChain.h"
   2 #include "TList.h"
   3 #include "TCanvas.h"
   4 #include "TParticle.h"
   5 #include "TLorentzVector.h"
   6 #include "TGraphErrors.h"
   7 #include "TH1.h"
   8 #include "TH2.h"
   9 #include "TH3.h"
  10 #include "TArrayF.h"
  11 #include "TF1.h"
  12 #include "TRandom.h"
  13 #include "TFile.h"
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliLog.h"
  17
  18 #include "AliEventPoolManager.h"
  19
  20 #include "AliAnalysisTaskToyModel.h"
  21 #include "AliBalance.h"
  22 #include "AliBalancePsi.h"
  23 #include "AliAnalysisTaskTriggeredBF.h"
  24
  25
  26 // Analysis task for the toy model analysis
  27 // Authors: Panos.Christakoglou@nikhef.nl
  28
  29 using std::cout;
  30 using std::endl;
  31
  32 ClassImp(AliAnalysisTaskToyModel)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskToyModel::AliAnalysisTaskToyModel()
  36 : TObject(),
  37   fUseDebug(kFALSE),
  38   fBalance(0),
  39   fRunShuffling(kFALSE), fShuffledBalance(0),
  40   fRunMixing(kFALSE), fMixedBalance(0), fPoolMgr(0),
  41   fList(0), fListBF(0), fListBFS(0), fListBFM(0),
  42   fHistEventStats(0),
  43   fHistNumberOfAcceptedParticles(0),
  44   fHistReactionPlane(0),
  45   fHistEtaTotal(0), fHistEta(0),
  46   fHistEtaPhiPos(0), fHistEtaPhiNeg(0),
  47   fHistRapidity(0),
  48   fHistRapidityPions(0), fHistRapidityKaons(0), fHistRapidityProtons(0),
  49   fHistPhi(0),
  50   fHistPhiPions(0), fHistPhiKaons(0), fHistPhiProtons(0),
  51   fHistPt(0),
  52   fHistPtPions(0), fHistPtKaons(0), fHistPtProtons(0),
  53   fTotalMultiplicityMean(100.), fTotalMultiplicitySigma(0.0),
  54   fNetChargeMean(0.0), fNetChargeSigma(0.0),
  55   fPtMin(0.0), fPtMax(100.0),
  56   fEtaMin(-1.0), fEtaMax(1.0),
  57   fUseAcceptanceParameterization(kFALSE), fAcceptanceParameterization(0),
  58   fSimulateDetectorEffects(kFALSE),
  59   fNumberOfInefficientSectors(0),
  60   fInefficiencyFactorInPhi(1.0),
  61   fNumberOfDeadSectors(0),
  62   fEfficiencyDropNearEtaEdges(kFALSE),
  63   fEfficiencyMatrix(0),
  64   fUseAllCharges(kFALSE), fParticleMass(0.0),
  65   fPtSpectraAllCharges(0), fTemperatureAllCharges(100.),
  66   fReactionPlane(0.0),
  67   fAzimuthalAngleAllCharges(0), fDirectedFlowAllCharges(0.0),
  68   fEllipticFlowAllCharges(0.0), fTriangularFlowAllCharges(0.0),
  69   fQuandrangularFlowAllCharges(0.0), fPentangularFlowAllCharges(0.0),
  70   fPionPercentage(0.8), fPionMass(0.0),
  71   fPtSpectraPions(0), fTemperaturePions(100.),
  72   fAzimuthalAnglePions(0), fDirectedFlowPions(0.0),
  73   fEllipticFlowPions(0.0), fTriangularFlowPions(0.0),
  74   fQuandrangularFlowPions(0.0), fPentangularFlowPions(0.0),
  75   fKaonPercentage(0.8), fKaonMass(0.0),
  76   fPtSpectraKaons(0), fTemperatureKaons(100.),
  77   fAzimuthalAngleKaons(0), fDirectedFlowKaons(0.0),
  78   fEllipticFlowKaons(0.0), fTriangularFlowKaons(0.0),
  79   fQuandrangularFlowKaons(0.0), fPentangularFlowKaons(0.0),
  80   fProtonPercentage(0.8), fProtonMass(0.0),
  81   fPtSpectraProtons(0), fTemperatureProtons(100.),
  82   fAzimuthalAngleProtons(0), fDirectedFlowProtons(0.0),
  83   fEllipticFlowProtons(0.0), fTriangularFlowProtons(0.0),
  84   fQuandrangularFlowProtons(0.0), fPentangularFlowProtons(0.0),
  85   fUseDynamicalCorrelations(kFALSE), fDynamicalCorrelationsPercentage(0.1),
  86   fUseJets(kFALSE), fPtAssoc(0) {
  87   // Constructor
  88 }
  89
  90 //________________________________________________________________________
  91 AliAnalysisTaskToyModel::~AliAnalysisTaskToyModel() {
  92   //Destructor
  93   if(fUseAllCharges) {
  94     delete fPtSpectraAllCharges;
  95     delete fAzimuthalAngleAllCharges;
  96   }
  97   else {
  98     delete fPtSpectraPions;
  99     delete fAzimuthalAnglePions;
 100     delete fPtSpectraKaons;
 101     delete fAzimuthalAngleKaons;
 102     delete fPtSpectraProtons;
 103     delete fAzimuthalAngleProtons;
 104   }
 105   if(fUseJets) delete fPtAssoc;
 106 }
 107
 108 //________________________________________________________________________
 109 void AliAnalysisTaskToyModel::Init() {
 110   //Initialize objects
 111   //==============gRandom Seed=======================//
 112   gRandom->SetSeed(0);  //seed is set to a random value (depending on machine clock)
 113   //==============gRandom Seed=======================//
 114
 115   //==============Particles and spectra==============//
 116   TParticle *pion = new TParticle();
 117   pion->SetPdgCode(211);
 118   fPionMass = pion->GetMass();
 119
 120   TParticle *kaon = new TParticle();
 121   kaon->SetPdgCode(321);
 122   fKaonMass = kaon->GetMass();
 123
 124   TParticle *proton = new TParticle();
 125   proton->SetPdgCode(2212);
 126   fProtonMass = proton->GetMass();
 127
 128   if(fUseAllCharges) {
 129     fParticleMass = fPionMass;
 130     fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
 131     fPtSpectraAllCharges->SetParName(0,"Mass");
 132     fPtSpectraAllCharges->SetParName(1,"Temperature");
 133     //fPtSpectraAllCharges = new TF1("fPtSpectraAllCharges","(x^2/TMath::Sqrt(TMath::Power(x,2) + TMath::Power(0.139,2)))*TMath::Power((1. + x/[0]),-[1])",0.,20.);
 134     //fPtSpectraAllCharges->SetParName(0,"pt0");
 135     //fPtSpectraAllCharges->SetParName(1,"b");
 136   }
 137   else {
 138     fPtSpectraPions = new TF1("fPtSpectraPions","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
 139     fPtSpectraPions->SetParName(0,"Mass");
 140     fPtSpectraPions->SetParName(1,"Temperature");
 141
 142     fPtSpectraKaons = new TF1("fPtSpectraKaons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
 143     fPtSpectraKaons->SetParName(0,"Mass");
 144     fPtSpectraKaons->SetParName(1,"Temperature");
 145
 146     fPtSpectraProtons = new TF1("fPtSpectraProtons","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,5.);
 147     fPtSpectraProtons->SetParName(0,"Mass");
 148     fPtSpectraProtons->SetParName(1,"Temperature");
 149   }
 150   //==============Particles and spectra==============//
 151
 152   //==============Flow values==============//
 153   if(fUseAllCharges) {
 154     if(fUseDebug) {
 155       Printf("Directed flow: %lf",fDirectedFlowAllCharges);
 156       Printf("Elliptic flow: %lf",fEllipticFlowAllCharges);
 157       Printf("Triangular flow: %lf",fTriangularFlowAllCharges);
 158       Printf("Quandrangular flow: %lf",fQuandrangularFlowAllCharges);
 159       Printf("Pentangular flow: %lf",fPentangularFlowAllCharges);
 160     }
 161
 162     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());
 163     fAzimuthalAngleAllCharges->SetParName(0,"Reaction Plane");
 164     fAzimuthalAngleAllCharges->SetParName(1,"Directed flow");
 165     fAzimuthalAngleAllCharges->SetParName(2,"Elliptic flow");
 166     fAzimuthalAngleAllCharges->SetParName(3,"Triangular flow");
 167     fAzimuthalAngleAllCharges->SetParName(4,"Quandrangular flow");
 168     fAzimuthalAngleAllCharges->SetParName(5,"Pentangular flow");
 169   }
 170   else {
 171     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());
 172     fAzimuthalAnglePions->SetParName(0,"Reaction Plane");
 173     fAzimuthalAnglePions->SetParName(1,"Directed flow");
 174     fAzimuthalAnglePions->SetParName(2,"Elliptic flow");
 175     fAzimuthalAnglePions->SetParName(3,"Triangular flow");
 176     fAzimuthalAnglePions->SetParName(4,"Quandrangular flow");
 177     fAzimuthalAnglePions->SetParName(5,"Pentangular flow");
 178
 179     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());
 180     fAzimuthalAngleKaons->SetParName(0,"Reaction Plane");
 181     fAzimuthalAngleKaons->SetParName(1,"Directed flow");
 182     fAzimuthalAngleKaons->SetParName(2,"Elliptic flow");
 183     fAzimuthalAngleKaons->SetParName(3,"Triangular flow");
 184     fAzimuthalAngleKaons->SetParName(4,"Quandrangular flow");
 185     fAzimuthalAngleKaons->SetParName(5,"Pentangular flow");
 186
 187     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());
 188     fAzimuthalAngleProtons->SetParName(0,"Reaction Plane");
 189     fAzimuthalAngleProtons->SetParName(1,"Directed flow");
 190     fAzimuthalAngleProtons->SetParName(2,"Elliptic flow");
 191     fAzimuthalAngleProtons->SetParName(3,"Triangular flow");
 192     fAzimuthalAngleProtons->SetParName(4,"Quandrangular flow");
 193     fAzimuthalAngleProtons->SetParName(5,"Pentangular flow");
 194   }
 195   //==============Flow values==============//
 196
 197   //===================Jets===================//
 198   if(fUseJets) {
 199     fPtAssoc = new TF1("fPtAssoc","x*TMath::Exp(-TMath::Power([0]*[0]+x*x,0.5)/[1])",0.,20.);
 200     fPtAssoc->SetParName(0,"pt0");
 201     fPtAssoc->SetParName(1,"b");
 202     fPtAssoc->SetParameter(0,0.139);
 203     fPtAssoc->SetParameter(1,0.5);
 204     fPtAssoc->SetLineColor(1);
 205   }
 206
 207   //===================Jets===================//
 208
 209   //==============Efficiency matrix==============//
 210   if(fSimulateDetectorEffects) SetupEfficiencyMatrix();
 211   //==============Efficiency matrix==============//
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAnalysisTaskToyModel::SetupEfficiencyMatrix() {
 216   //Setup the efficiency matrix
 217   TH1F *hPt = new TH1F("hPt","",200,0.1,20.1);
 218   TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95);
 219   TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi());
 220   fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","",
 221                                hEta->GetNbinsX(),
 222                                hEta->GetXaxis()->GetXmin(),
 223                                hEta->GetXaxis()->GetXmax(),
 224                                hPt->GetNbinsX(),
 225                                hPt->GetXaxis()->GetXmin(),
 226                                hPt->GetXaxis()->GetXmax(),
 227                                hPhi->GetNbinsX(),
 228                                hPhi->GetXaxis()->GetXmin(),
 229                                hPhi->GetXaxis()->GetXmax());
 230
 231   //Efficiency in pt
 232   Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80};
 233   for(Int_t i=1;i<=20;i++) {
 234     hPt->SetBinContent(i,epsilon[i-1]);
 235     hPt->SetBinError(i,0.01);
 236   }
 237   for(Int_t i=21;i<=200;i++) {
 238     hPt->SetBinContent(i,epsilon[19]);
 239     hPt->SetBinError(i,0.01);
 240   }
 241
 242   //Efficiency in eta
 243   for(Int_t i=1;i<=hEta->GetNbinsX();i++) {
 244     hEta->SetBinContent(i,1.0);
 245     hEta->SetBinError(i,0.01);
 246   }
 247   if(fEfficiencyDropNearEtaEdges) {
 248     hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8);
 249     hEta->SetBinContent(3,0.9);
 250     hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8);
 251     hEta->SetBinContent(20,0.7);
 252   }
 253
 254   //Efficiency in phi
 255   for(Int_t i=1;i<=hPhi->GetNbinsX();i++) {
 256     hPhi->SetBinContent(i,1.0);
 257     hPhi->SetBinError(i,0.01);
 258   }
 259   for(Int_t i=1;i<=fNumberOfInefficientSectors;i++)
 260     hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi);
 261   for(Int_t i=1;i<=fNumberOfDeadSectors;i++)
 262     hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0);
 263
 264   //Fill the 3D efficiency map
 265   for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) {
 266     //cout<<"==================================="<<endl;
 267     for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {
 268       //cout<<"==================================="<<endl;
 269       for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {
 270         fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));
 271         //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;
 272       }
 273     }
 274   }
 275 }
 276
 277 //________________________________________________________________________
 278 void AliAnalysisTaskToyModel::CreateOutputObjects() {
 279   // Create histograms
 280   // Called once
 281
 282   // global switch disabling the reference
 283   // (to avoid "Replacing existing TH1" if several wagons are created in train)
 284   Bool_t oldStatus = TH1::AddDirectoryStatus();
 285   TH1::AddDirectory(kFALSE);
 286
 287   if(!fBalance) {
 288     fBalance = new AliBalancePsi();
 289     fBalance->SetDeltaEtaMax(2.0);
 290   }
 291   if(fRunShuffling) {
 292     if(!fShuffledBalance) {
 293       fShuffledBalance = new AliBalancePsi();
 294       fShuffledBalance->SetDeltaEtaMax(2.0);
 295     }
 296   }
 297   if(fRunMixing) {
 298     if(!fMixedBalance) {
 299       fMixedBalance = new AliBalancePsi();
 300       fMixedBalance->SetDeltaEtaMax(2.0);
 301     }
 302   }
 303
 304   //QA list
 305   fList = new TList();
 306   fList->SetName("listQA");
 307   fList->SetOwner();
 308
 309   //Balance Function list
 310   fListBF = new TList();
 311   fListBF->SetName("listBF");
 312   fListBF->SetOwner();
 313
 314   if(fRunShuffling) {
 315     fListBFS = new TList();
 316     fListBFS->SetName("listBFShuffled");
 317     fListBFS->SetOwner();
 318   }
 319
 320   if(fRunMixing) {
 321     fListBFM = new TList();
 322     fListBFM->SetName("listBFMixed");
 323     fListBFM->SetOwner();
 324   }
 325
 326   //==============QA================//
 327   //Event stats.
 328   TString gCutName[4] = {"Total","Offline trigger",
 329                          "Vertex","Analyzed"};
 330   fHistEventStats = new TH1F("fHistEventStats",
 331                              "Event statistics;;N_{events}",
 332                              4,0.5,4.5);
 333   for(Int_t i = 1; i <= 4; i++)
 334     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 335   fList->Add(fHistEventStats);
 336
 337   fHistNumberOfAcceptedParticles = new TH1F("fHistNumberOfAcceptedParticles",";N_{acc.};Entries",10000,0,10000);
 338   fList->Add(fHistNumberOfAcceptedParticles);
 339
 340   fHistReactionPlane = new TH1F("fHistReactionPlane","Reaction plane angle;#Psi [rad];Entries",1000,0.,2.*TMath::Pi());
 341   fList->Add(fHistReactionPlane);
 342
 343   //Pseudo-rapidity
 344   fHistEtaTotal = new TH1F("fHistEtaTotal","Pseudo-rapidity (full phase space);#eta;Entries",1000,-15.,15.);
 345   fList->Add(fHistEtaTotal);
 346
 347   fHistEta = new TH1F("fHistEta","Pseudo-rapidity (acceptance);#eta;Entries",1000,-1.5,1.5);
 348   fList->Add(fHistEta);
 349
 350   fHistEtaPhiPos = new TH2F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
 351   fList->Add(fHistEtaPhiPos);
 352   fHistEtaPhiNeg = new TH2F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#varphi (rad)",80,-2.,2.,72,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
 353   fList->Add(fHistEtaPhiNeg);
 354
 355   //Rapidity
 356   fHistRapidity = new TH1F("fHistRapidity","Rapidity (acceptance);y;Entries",1000,-1.5,1.5);
 357   fList->Add(fHistRapidity);
 358   fHistRapidityPions = new TH1F("fHistRapidityPions","Rapidity (acceptance - pions);y;Entries",1000,-1.5,1.5);
 359   fList->Add(fHistRapidityPions);
 360   fHistRapidityKaons = new TH1F("fHistRapidityKaons","Rapidity (acceptance - kaons);y;Entries",1000,-1.5,1.5);
 361   fList->Add(fHistRapidityKaons);
 362   fHistRapidityProtons = new TH1F("fHistRapidityProtons","Rapidity (acceptance - protons);y;Entries",1000,-1.5,1.5);
 363   fList->Add(fHistRapidityProtons);
 364
 365   //Phi
 366   fHistPhi = new TH1F("fHistPhi","Phi (acceptance);#phi (rad);Entries",1000,0.,2*TMath::Pi());
 367   fList->Add(fHistPhi);
 368
 369   fHistPhiPions = new TH1F("fHistPhiPions","Phi (acceptance - pions);#phi (rad);Entries",1000,0.,2*TMath::Pi());
 370   fList->Add(fHistPhiPions);
 371   fHistPhiKaons = new TH1F("fHistPhiKaons","Phi (acceptance - kaons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
 372   fList->Add(fHistPhiKaons);
 373   fHistPhiProtons = new TH1F("fHistPhiProtons","Phi (acceptance - protons);#phi (rad);Entries",1000,0.,2*TMath::Pi());
 374   fList->Add(fHistPhiProtons);
 375
 376   //Pt
 377   fHistPt = new TH1F("fHistPt","Pt (acceptance);p_{t} (GeV/c);Entries",1000,0.,10.);
 378   fList->Add(fHistPt);
 379
 380   fHistPtPions = new TH1F("fHistPtPions","Pt (acceptance - pions);p_{t} (GeV/c);Entries",1000,0.,10.);
 381   fList->Add(fHistPtPions);
 382   fHistPtKaons = new TH1F("fHistPtKaons","Pt (acceptance - kaons);p_{t} (GeV/c);Entries",1000,0.,10.);
 383   fList->Add(fHistPtKaons);
 384   fHistPtProtons = new TH1F("fHistPtProtons","Pt (acceptance - protons);p_{t} (GeV/c);Entries",1000,0.,10.);
 385   fList->Add(fHistPtProtons);
 386
 387   if(fEfficiencyMatrix) fList->Add(fEfficiencyMatrix);
 388
 389   //==============Balance function histograms================//
 390   // Initialize histograms if not done yet
 391   if(!fBalance->GetHistNp()){
 392     AliWarning("Histograms not yet initialized! --> Will be done now");
 393     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 394     fBalance->InitHistograms();
 395   }
 396
 397   if(fRunShuffling) {
 398     if(!fShuffledBalance->GetHistNp()) {
 399       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
 400       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 401       fShuffledBalance->InitHistograms();
 402     }
 403   }
 404
 405   fListBF->Add(fBalance->GetHistNp());
 406   fListBF->Add(fBalance->GetHistNn());
 407   fListBF->Add(fBalance->GetHistNpn());
 408   fListBF->Add(fBalance->GetHistNnn());
 409   fListBF->Add(fBalance->GetHistNpp());
 410   fListBF->Add(fBalance->GetHistNnp());
 411
 412   if(fRunShuffling) {
 413     fListBFS->Add(fShuffledBalance->GetHistNp());
 414     fListBFS->Add(fShuffledBalance->GetHistNn());
 415     fListBFS->Add(fShuffledBalance->GetHistNpn());
 416     fListBFS->Add(fShuffledBalance->GetHistNnn());
 417     fListBFS->Add(fShuffledBalance->GetHistNpp());
 418     fListBFS->Add(fShuffledBalance->GetHistNnp());
 419   }
 420
 421   if(fRunMixing) {
 422     fListBFM->Add(fMixedBalance->GetHistNp());
 423     fListBFM->Add(fMixedBalance->GetHistNn());
 424     fListBFM->Add(fMixedBalance->GetHistNpn());
 425     fListBFM->Add(fMixedBalance->GetHistNnn());
 426     fListBFM->Add(fMixedBalance->GetHistNpp());
 427     fListBFM->Add(fMixedBalance->GetHistNnp());
 428   }
 429
 430   // Event Mixing
 431   if(fRunMixing){
 432     Int_t fMixingTracks = 2000;
 433     Int_t trackDepth = fMixingTracks;
 434     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 435
 436     // centrality bins
 437     // Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
 438     // Double_t* centbins        = centralityBins;
 439     // Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;
 440
 441     // multiplicity bins
 442     Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
 443     Double_t* multbins        = multiplicityBins;
 444     Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;
 445
 446     // Zvtx bins
 447     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
 448     Double_t* vtxbins     = vertexBins;
 449     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;
 450
 451     // Event plane angle (Psi) bins
 452     // Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!
 453     // Double_t* psibins     = psiBins;
 454     // Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;
 455
 456     // // run the event mixing also in bins of event plane (statistics!)
 457     // if(fRunMixingEventPlane){
 458     //   if(fEventClass=="Multiplicity"){
 459     //  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);
 460     //   }
 461     //   else{
 462     //  fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);
 463     //   }
 464     // }
 465     // else{
 466     //if(fEventClass=="Multiplicity"){
 467     fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);
 468     //}
 469     //else{
 470     //fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);
 471     //}
 472   }
 473
 474   TH1::AddDirectory(oldStatus);
 475 }
 476
 477 //________________________________________________________________________
 478 void AliAnalysisTaskToyModel::Run(Int_t nEvents) {
 479   // Main loop
 480   // Called for each event
 481   Short_t vCharge = 0;
 482   Float_t vY = 0.0;
 483   Float_t vEta = 0.0;
 484   Float_t vPhi = 0.0;
 485   Float_t vP[3] = {0.,0.,0.};
 486   Float_t vPt = 0.0;
 487   Float_t vE = 0.0;
 488   Bool_t isPion = kFALSE, isKaon = kFALSE, isProton = kFALSE;
 489
 490   Double_t gDecideCharge = 0.;
 491
 492   if(fUseAllCharges) {
 493     //fPtSpectraAllCharges->SetParameter(0,fParticleMass);
 494     fPtSpectraAllCharges->SetParameter(0,1.05);
 495     fPtSpectraAllCharges->SetParameter(1,fTemperatureAllCharges);
 496
 497     fAzimuthalAngleAllCharges->SetParameter(1,fDirectedFlowAllCharges);
 498     fAzimuthalAngleAllCharges->SetParameter(2,fEllipticFlowAllCharges);
 499     fAzimuthalAngleAllCharges->SetParameter(3,fTriangularFlowAllCharges);
 500     fAzimuthalAngleAllCharges->SetParameter(4,fQuandrangularFlowAllCharges);
 501     fAzimuthalAngleAllCharges->SetParameter(5,fPentangularFlowAllCharges);
 502   }
 503   else {
 504     fPtSpectraPions->SetParameter(0,fPionMass);
 505     fPtSpectraPions->SetParameter(1,fTemperaturePions);
 506     fPtSpectraKaons->SetParameter(0,fKaonMass);
 507     fPtSpectraKaons->SetParameter(1,fTemperatureKaons);
 508     fPtSpectraProtons->SetParameter(0,fProtonMass);
 509     fPtSpectraProtons->SetParameter(1,fTemperatureProtons);
 510
 511     fAzimuthalAnglePions->SetParameter(1,fDirectedFlowPions);
 512     fAzimuthalAnglePions->SetParameter(2,fEllipticFlowPions);
 513     fAzimuthalAnglePions->SetParameter(3,fTriangularFlowPions);
 514     fAzimuthalAnglePions->SetParameter(4,fQuandrangularFlowPions);
 515     fAzimuthalAnglePions->SetParameter(5,fPentangularFlowPions);
 516
 517     fAzimuthalAngleKaons->SetParameter(1,fDirectedFlowKaons);
 518     fAzimuthalAngleKaons->SetParameter(2,fEllipticFlowKaons);
 519     fAzimuthalAngleKaons->SetParameter(3,fTriangularFlowKaons);
 520     fAzimuthalAngleKaons->SetParameter(4,fQuandrangularFlowKaons);
 521     fAzimuthalAngleKaons->SetParameter(5,fPentangularFlowKaons);
 522
 523     fAzimuthalAngleProtons->SetParameter(1,fDirectedFlowProtons);
 524     fAzimuthalAngleProtons->SetParameter(2,fEllipticFlowProtons);
 525     fAzimuthalAngleProtons->SetParameter(3,fTriangularFlowProtons);
 526     fAzimuthalAngleProtons->SetParameter(4,fQuandrangularFlowProtons);
 527     fAzimuthalAngleProtons->SetParameter(5,fPentangularFlowProtons);
 528   }
 529
 530
 531   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
 532     // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
 533     //vector<Double_t> *chargeVectorShuffle[9];   // this will be shuffled
 534     //vector<Double_t> *chargeVector[9];          // original charge
 535     //for(Int_t i = 0; i < 9; i++){
 536     //chargeVectorShuffle[i] = new vector<Double_t>;
 537     //chargeVector[i]        = new vector<Double_t>;
 538     //}
 539
 540     // TObjArray for the accepted particles
 541     // (has to be done here, otherwise mxing with event pool does not work, overwriting pointers!)
 542     TObjArray *tracksMain = new TObjArray();
 543     tracksMain->SetOwner(kTRUE);
 544     TObjArray *tracksMixing = 0x0;
 545     if(fRunMixing) {
 546       tracksMixing = new TObjArray();
 547       tracksMixing->SetOwner(kTRUE);
 548     }
 549
 550     tracksMain->Clear();
 551     if(fRunMixing) tracksMixing->Clear();
 552
 553     fHistEventStats->Fill(1);
 554     fHistEventStats->Fill(2);
 555     fHistEventStats->Fill(3);
 556
 557     if((iEvent%1000) == 0)
 558       cout<<"Event: "<<iEvent<<"/"<<nEvents<<endl;
 559
 560     //Multiplicities
 561     Int_t nMultiplicity = (Int_t)(gRandom->Gaus(fTotalMultiplicityMean,fTotalMultiplicitySigma));
 562     Int_t nNetCharge = (Int_t)(gRandom->Gaus(fNetChargeMean,fNetChargeSigma));
 563
 564     Int_t nGeneratedPositive = (Int_t)((nMultiplicity/2) + nNetCharge);
 565     Int_t nGeneratedNegative = nMultiplicity - nGeneratedPositive;
 566     if(fUseDebug)
 567       Printf("Total multiplicity: %d - Generated positive: %d - Generated negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
 568
 569     //Randomization of the reaction plane
 570     fReactionPlane = 2.0*TMath::Pi()*gRandom->Rndm();
 571     //fReactionPlane = 0.0;
 572     if(fUseAllCharges)
 573       fAzimuthalAngleAllCharges->SetParameter(0,fReactionPlane);
 574     else {
 575       fAzimuthalAnglePions->SetParameter(0,fReactionPlane);
 576       fAzimuthalAngleKaons->SetParameter(0,fReactionPlane);
 577       fAzimuthalAngleProtons->SetParameter(0,fReactionPlane);
 578     }
 579
 580     Int_t gNumberOfAcceptedParticles = 0;
 581     Int_t gNumberOfAcceptedPositiveParticles = 0;
 582     Int_t gNumberOfAcceptedNegativeParticles = 0;
 583     Int_t nGeneratedPions = 0, nGeneratedKaons = 0, nGeneratedProtons = 0;
 584
 585     //Generate particles
 586     for(Int_t iParticleCount = 0; iParticleCount < nMultiplicity; iParticleCount++) {
 587       isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
 588       if(fUseDebug)
 589         Printf("Generating positive: %d(%d)",iParticleCount+1,nGeneratedPositive);
 590
 591       //Pseudo-rapidity sampled from a Gaussian centered @ 0
 592       vEta = gRandom->Gaus(0.0,4.0);
 593
 594       //Fill QA histograms (full phase space)
 595       fHistEtaTotal->Fill(vEta);
 596
 597       //Decide the charge
 598       gDecideCharge = gRandom->Rndm();
 599       if(gDecideCharge <= 1.*nGeneratedPositive/nMultiplicity)
 600         vCharge = 1;
 601       else
 602         vCharge = -1;
 603
 604       //Acceptance
 605       if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
 606
 607       if(!fUseAllCharges) {
 608         //Decide the specie
 609         Double_t randomNumberSpecies = gRandom->Rndm();
 610         if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
 611           nGeneratedPions += 1;
 612           vPt = fPtSpectraPions->GetRandom();
 613           vPhi = fAzimuthalAnglePions->GetRandom();
 614           fParticleMass = fPionMass;
 615           isPion = kTRUE;
 616         }
 617         else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
 618           nGeneratedKaons += 1;
 619           vPt = fPtSpectraKaons->GetRandom();
 620           vPhi = fAzimuthalAngleKaons->GetRandom();
 621           fParticleMass = fKaonMass;
 622           isKaon = kTRUE;
 623         }
 624         else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
 625           nGeneratedProtons += 1;
 626           vPt = fPtSpectraProtons->GetRandom();
 627           vPhi = fAzimuthalAngleProtons->GetRandom();
 628           fParticleMass = fProtonMass;
 629           isProton = kTRUE;
 630         }
 631       }
 632       else {
 633         vPt = fPtSpectraAllCharges->GetRandom();
 634         vPhi = fAzimuthalAngleAllCharges->GetRandom();
 635       }
 636
 637       vP[0] = vPt*TMath::Cos(vPhi);
 638       vP[1] = vPt*TMath::Sin(vPhi);
 639       vP[2] = vPt*TMath::SinH(vEta);
 640       vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
 641                         TMath::Power(vP[0],2) +
 642                         TMath::Power(vP[1],2) +
 643                         TMath::Power(vP[2],2));
 644
 645       vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
 646
 647       //pt coverage
 648       if((vPt < fPtMin) || (vPt > fPtMax)) continue;
 649       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
 650
 651       //acceptance filter
 652       if(fUseAcceptanceParameterization) {
 653         Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 654         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
 655           continue;
 656       }
 657
 658       //Detector effects
 659       if(fSimulateDetectorEffects) {
 660         Double_t randomNumber = gRandom->Rndm();
 661         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(vEta,vPt,vPhi)))
 662           continue;
 663       }
 664
 665       //Fill QA histograms (acceptance);
 666       if(vCharge > 0) {
 667         gNumberOfAcceptedPositiveParticles += 1;
 668         if(vPhi > 3.*TMath::Pi()/2.)
 669           fHistEtaPhiPos->Fill(vEta,vPhi-2.*TMath::Pi());
 670         else
 671           fHistEtaPhiPos->Fill(vEta,vPhi);
 672       }
 673       else {
 674         gNumberOfAcceptedNegativeParticles += 1;
 675         if(vPhi > 3.*TMath::Pi()/2.)
 676           fHistEtaPhiNeg->Fill(vEta,vPhi-2.*TMath::Pi());
 677         else
 678           fHistEtaPhiNeg->Fill(vEta,vPhi);
 679       }
 680
 681       fHistEta->Fill(vEta);
 682       fHistRapidity->Fill(vY);
 683       fHistPhi->Fill(vPhi);
 684       fHistPt->Fill(vPt);
 685       if(isPion) {
 686         fHistRapidityPions->Fill(vY);
 687         fHistPhiPions->Fill(vPhi);
 688         fHistPtPions->Fill(vPt);
 689       }
 690       else if(isKaon) {
 691         fHistRapidityKaons->Fill(vY);
 692         fHistPhiKaons->Fill(vPhi);
 693         fHistPtKaons->Fill(vPt);
 694       }
 695       else if(isProton) {
 696         fHistRapidityProtons->Fill(vY);
 697         fHistPhiProtons->Fill(vPhi);
 698         fHistPtProtons->Fill(vPt);
 699       }
 700       gNumberOfAcceptedParticles += 1;
 701
 702       // add the track to the TObjArray
 703       tracksMain->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
 704       if(fRunMixing)
 705         tracksMixing->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, 1.0));
 706     }//generated positive particle loop
 707
 708     //Jets
 709     if(fUseJets) {
 710       const Int_t nAssociated = 3;
 711
 712       Double_t gPtTrig1 = 0., gPtTrig2 = 0., gPtAssoc = 0.;
 713       Double_t gPhiTrig1 = 0., gPhiTrig2 = 0., gPhiAssoc = 0.;
 714       Double_t gEtaTrig1 = 0., gEtaTrig2 = 0., gEtaAssoc = 0.;
 715       Short_t gChargeTrig1 = 0, gChargeTrig2 = 0, gChargeAssoc = 0;
 716
 717       Double_t gJetCone = 0.2;
 718
 719       //First leading particle
 720       gPtTrig1 = gRandom->Uniform(3.,5.);
 721       gEtaTrig1 = gRandom->Uniform(-0.8,0.8);
 722       gPhiTrig1 = gRandom->Uniform(0.,TMath::TwoPi());
 723
 724       //Decide the charge
 725       gDecideCharge = gRandom->Rndm();
 726       if(gDecideCharge <= 0.5)
 727         gChargeTrig1 = 1;
 728       else
 729         gChargeTrig1 = -1;
 730
 731       //Acceptance
 732       if((gEtaTrig1 < fEtaMin) || (gEtaTrig1 > fEtaMax)) continue;
 733       //pt coverage
 734       if((gPtTrig1 < fPtMin) || (gPtTrig1 > fPtMax)) continue;
 735       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
 736
 737       //acceptance filter
 738       if(fUseAcceptanceParameterization) {
 739         Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 740         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig1))
 741           continue;
 742       }
 743
 744       //Detector effects
 745       if(fSimulateDetectorEffects) {
 746         Double_t randomNumber = gRandom->Rndm();
 747         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig1,gPtTrig1,gPhiTrig1)))
 748           continue;
 749       }
 750
 751       gNumberOfAcceptedParticles += 1;
 752
 753       // add the track to the TObjArray
 754       tracksMain->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
 755       if(fRunMixing)
 756         tracksMixing->Add(new AliBFBasicParticle(gEtaTrig1, gPhiTrig1, gPtTrig1, gChargeTrig1, 1.0));
 757
 758       Int_t iAssociated = 0;
 759       while(iAssociated < nAssociated) {
 760         gPtAssoc = fPtAssoc->GetRandom();
 761         if(gPtAssoc < gPtTrig1) {
 762           gEtaAssoc = gRandom->Uniform(gEtaTrig1 - gJetCone/2.,gEtaTrig1 + gJetCone/2.);
 763           gPhiAssoc = gRandom->Uniform(gPhiTrig1 - gJetCone/2.,gPhiTrig1 + gJetCone/2.);
 764           if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
 765           else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
 766
 767           iAssociated += 1;
 768
 769           gDecideCharge = gRandom->Rndm();
 770           if(gDecideCharge <= 0.5)
 771             gChargeAssoc = 1;
 772           else
 773             gChargeAssoc = -1;
 774
 775           //Acceptance
 776           if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
 777           //pt coverage
 778           if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
 779           //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
 780
 781           //acceptance filter
 782           if(fUseAcceptanceParameterization) {
 783             Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 784             if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
 785               continue;
 786           }
 787
 788           //Detector effects
 789           if(fSimulateDetectorEffects) {
 790             Double_t randomNumber = gRandom->Rndm();
 791             if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
 792               continue;
 793           }
 794
 795           gNumberOfAcceptedParticles += 1;
 796
 797           // add the track to the TObjArray
 798           tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
 799           if(fRunMixing)
 800             tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
 801         }//pt,assoc < pt,trig
 802       }//associated
 803
 804       //back2back
 805       gPtTrig2 = gPtTrig1;
 806       gEtaTrig2 = -gEtaTrig1;
 807       gPhiTrig2 = TMath::Pi() + gPhiTrig1;
 808       if(gPhiTrig2 < 0.) gPhiTrig2 += TMath::TwoPi();
 809       else if(gPhiTrig2 > TMath::TwoPi()) gPhiTrig2 -= TMath::TwoPi();
 810
 811       //Decide the charge
 812       gDecideCharge = gRandom->Rndm();
 813       if(gDecideCharge <= 0.5)
 814         gChargeTrig2 = 1;
 815       else
 816         gChargeTrig2 = -1;
 817
 818       //Acceptance
 819       if((gEtaTrig2 < fEtaMin) || (gEtaTrig2 > fEtaMax)) continue;
 820       //pt coverage
 821       if((gPtTrig2 < fPtMin) || (gPtTrig2 > fPtMax)) continue;
 822       //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
 823
 824       //acceptance filter
 825       if(fUseAcceptanceParameterization) {
 826         Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 827         if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtTrig2))
 828           continue;
 829       }
 830
 831       //Detector effects
 832       if(fSimulateDetectorEffects) {
 833         Double_t randomNumber = gRandom->Rndm();
 834         if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaTrig2,gPtTrig2,gPhiTrig2)))
 835           continue;
 836       }
 837
 838       gNumberOfAcceptedParticles += 1;
 839
 840       // add the track to the TObjArray
 841       tracksMain->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
 842       if(fRunMixing)
 843         tracksMixing->Add(new AliBFBasicParticle(gEtaTrig2, gPhiTrig2, gPtTrig2, gChargeTrig2, 1.0));
 844
 845       iAssociated = 0;
 846       while(iAssociated < nAssociated) {
 847         gPtAssoc = fPtAssoc->GetRandom();
 848         if(gPtAssoc < gPtTrig2) {
 849           gEtaAssoc = gRandom->Uniform(gEtaTrig2 - gJetCone/2.,gEtaTrig2 + gJetCone/2.);
 850           gPhiAssoc = gRandom->Uniform(gPhiTrig2 - gJetCone/2.,gPhiTrig2 + gJetCone/2.);
 851           if(gPhiAssoc < 0.) gPhiAssoc += TMath::TwoPi();
 852           else if(gPhiAssoc > TMath::TwoPi()) gPhiAssoc -= TMath::TwoPi();
 853
 854           iAssociated += 1;
 855
 856           gDecideCharge = gRandom->Rndm();
 857           if(gDecideCharge <= 0.5)
 858             gChargeAssoc = 1;
 859           else
 860             gChargeAssoc = -1;
 861
 862           //Acceptance
 863           if((gEtaAssoc < fEtaMin) || (gEtaAssoc > fEtaMax)) continue;
 864           //pt coverage
 865           if((gPtAssoc < fPtMin) || (gPtAssoc > fPtMax)) continue;
 866           //Printf("pt: %lf - mins: %lf - max: %lf",vPt,fPtMin,fPtMax);
 867
 868           //acceptance filter
 869           if(fUseAcceptanceParameterization) {
 870             Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 871             if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(gPtAssoc))
 872               continue;
 873           }
 874
 875           //Detector effects
 876           if(fSimulateDetectorEffects) {
 877             Double_t randomNumber = gRandom->Rndm();
 878             if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(gEtaAssoc,gPtAssoc,gPhiAssoc)))
 879               continue;
 880           }
 881
 882           gNumberOfAcceptedParticles += 1;
 883
 884           // add the track to the TObjArray
 885           tracksMain->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
 886           if(fRunMixing)
 887             tracksMixing->Add(new AliBFBasicParticle(gEtaAssoc, gPhiAssoc, gPtAssoc, gChargeAssoc, 1.0));
 888         }//pt,assoc < pt,trig
 889       }//associated
 890     }//Jet usage
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900     //Dynamical correlations
 901     Int_t nGeneratedPositiveDynamicalCorrelations = 0;
 902     Int_t nGeneratedNegativeDynamicalCorrelations = 0;
 903     /*Double_t vChargePrime = 0;
 904     Double_t vYPrime = 0.0;
 905     Double_t vEtaPrime = 0.0;
 906     Double_t vPhiPrime = 0.0;
 907     Double_t vPPrime[3] = {0.,0.,0.};
 908     Double_t vPtPrime = 0.0;
 909     Double_t vEPrime = 0.0;
 910     //Generate "correlated" particles
 911     if(fUseDynamicalCorrelations) {
 912       Int_t gNumberOfDynamicalCorrelations = (Int_t)(0.5*gNumberOfAcceptedParticles*fDynamicalCorrelationsPercentage);
 913       for(Int_t iDynamicalCorrelations = 0; iDynamicalCorrelations < gNumberOfDynamicalCorrelations; iDynamicalCorrelations++) {
 914         isPion = kFALSE; isKaon = kFALSE; isProton = kFALSE;
 915
 916         //Pseudo-rapidity sampled from a Gaussian centered @ 0
 917         vEta = gRandom->Gaus(0.0,0.1);
 918         vCharge = 1.0;
 919         nGeneratedPositiveDynamicalCorrelations += 1;
 920
 921         vEtaPrime = -vEta;
 922         vChargePrime = -1.0;
 923         nGeneratedNegativeDynamicalCorrelations += 1;
 924
 925         //Acceptance
 926         if((vEta < fEtaMin) || (vEta > fEtaMax)) continue;
 927         if((vEtaPrime < fEtaMin) || (vEtaPrime > fEtaMax)) continue;
 928
 929         if(!fUseAllCharges) {
 930           //Decide the specie
 931           Double_t randomNumberSpecies = gRandom->Rndm();
 932           if((randomNumberSpecies >= 0.0)&&(randomNumberSpecies < fPionPercentage)) {
 933             nGeneratedPions += 1;
 934             vPt = fPtSpectraPions->GetRandom();
 935             vPhi = fAzimuthalAnglePions->GetRandom();
 936             fParticleMass = fPionMass;
 937             isPion = kTRUE;
 938           }
 939           else if((randomNumberSpecies >= fPionPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage)) {
 940             nGeneratedKaons += 1;
 941             vPt = fPtSpectraKaons->GetRandom();
 942             vPhi = fAzimuthalAngleKaons->GetRandom();
 943             fParticleMass = fKaonMass;
 944             isKaon = kTRUE;
 945           }
 946           else if((randomNumberSpecies >= fPionPercentage + fKaonPercentage)&&(randomNumberSpecies < fPionPercentage + fKaonPercentage + fProtonPercentage)) {
 947             nGeneratedProtons += 1;
 948             vPt = fPtSpectraProtons->GetRandom();
 949             vPtPrime = vPt;
 950             vPhi = fAzimuthalAngleProtons->GetRandom();
 951             fParticleMass = fProtonMass;
 952             isProton = kTRUE;
 953           }
 954         }
 955         else {
 956           vPt = fPtSpectraAllCharges->GetRandom();
 957           vPhi = fAzimuthalAngleAllCharges->GetRandom();
 958         }
 959         vPtPrime = vPt;
 960         vPhiPrime = vPhi;
 961
 962         vP[0] = vPt*TMath::Cos(vPhi);
 963         vP[1] = vPt*TMath::Sin(vPhi);
 964         vP[2] = vPt*TMath::SinH(vEta);
 965         vE = TMath::Sqrt(TMath::Power(fParticleMass,2) +
 966                           TMath::Power(vP[0],2) +
 967                           TMath::Power(vP[1],2) +
 968                           TMath::Power(vP[2],2));
 969
 970         vY = 0.5*TMath::Log((vE + vP[2])/(vE - vP[2]));
 971
 972         vPPrime[0] = vPtPrime*TMath::Cos(vPhiPrime);
 973         vPPrime[1] = vPtPrime*TMath::Sin(vPhiPrime);
 974         vPPrime[2] = vPtPrime*TMath::SinH(vEtaPrime);
 975         vEPrime = TMath::Sqrt(TMath::Power(fParticleMass,2) +
 976                           TMath::Power(vPPrime[0],2) +
 977                           TMath::Power(vPPrime[1],2) +
 978                           TMath::Power(vPPrime[2],2));
 979
 980         vYPrime = 0.5*TMath::Log((vEPrime + vPPrime[2])/(vEPrime - vPPrime[2]));
 981
 982         //pt coverage
 983         if((vPt < fPtMin) || (vPt > fPtMax)) continue;
 984         if((vPtPrime < fPtMin) || (vPtPrime > fPtMax)) continue;
 985
 986         //acceptance filter
 987         if(fUseAcceptanceParameterization) {
 988           Double_t gRandomNumberForAcceptance = gRandom->Rndm();
 989           if(gRandomNumberForAcceptance > fAcceptanceParameterization->Eval(vPt))
 990             continue;
 991
 992           Double_t gRandomNumberForAcceptancePrime = gRandom->Rndm();
 993           if(gRandomNumberForAcceptancePrime > fAcceptanceParameterization->Eval(vPtPrime))
 994             continue;
 995         }
 996
 997         // fill charge vector (positive)
 998         chargeVector[0]->push_back(vCharge);
 999         chargeVector[1]->push_back(vY);
1000         chargeVector[2]->push_back(vEta);
1001         chargeVector[3]->push_back(vPhi);
1002         chargeVector[4]->push_back(vP[0]);
1003         chargeVector[5]->push_back(vP[1]);
1004         chargeVector[6]->push_back(vP[2]);
1005         chargeVector[7]->push_back(vPt);
1006         chargeVector[8]->push_back(vE);
1007
1008         if(fRunShuffling) {
1009           chargeVectorShuffle[0]->push_back(vCharge);
1010           chargeVectorShuffle[1]->push_back(vY);
1011           chargeVectorShuffle[2]->push_back(vEta);
1012           chargeVectorShuffle[3]->push_back(vPhi);
1013           chargeVectorShuffle[4]->push_back(vP[0]);
1014           chargeVectorShuffle[5]->push_back(vP[1]);
1015           chargeVectorShuffle[6]->push_back(vP[2]);
1016           chargeVectorShuffle[7]->push_back(vPt);
1017           chargeVectorShuffle[8]->push_back(vE);
1018         }
1019
1020         // fill charge vector (negative)
1021         chargeVector[0]->push_back(vChargePrime);
1022         chargeVector[1]->push_back(vYPrime);
1023         chargeVector[2]->push_back(vEtaPrime);
1024         chargeVector[3]->push_back(vPhiPrime);
1025         chargeVector[4]->push_back(vPPrime[0]);
1026         chargeVector[5]->push_back(vPPrime[1]);
1027         chargeVector[6]->push_back(vPPrime[2]);
1028         chargeVector[7]->push_back(vPtPrime);
1029         chargeVector[8]->push_back(vEPrime);
1030
1031         if(fRunShuffling) {
1032           chargeVectorShuffle[0]->push_back(vChargePrime);
1033           chargeVectorShuffle[1]->push_back(vYPrime);
1034           chargeVectorShuffle[2]->push_back(vEtaPrime);
1035           chargeVectorShuffle[3]->push_back(vPhiPrime);
1036           chargeVectorShuffle[4]->push_back(vPPrime[0]);
1037           chargeVectorShuffle[5]->push_back(vPPrime[1]);
1038           chargeVectorShuffle[6]->push_back(vPPrime[2]);
1039           chargeVectorShuffle[7]->push_back(vPtPrime);
1040           chargeVectorShuffle[8]->push_back(vEPrime);
1041         }
1042
1043         gNumberOfAcceptedParticles += 2;
1044         }//loop over the dynamical correlations
1045         }*/// usage of dynamical correlations
1046
1047     if(fUseDebug) {
1048       Printf("=======================================================");
1049       Printf("Total: %d - Total positive: %d - Total negative: %d",nMultiplicity,nGeneratedPositive,nGeneratedNegative);
1050       Printf("Accepted positive: %d - Accepted negative: %d",gNumberOfAcceptedPositiveParticles,gNumberOfAcceptedNegativeParticles);
1051       Printf("Correlations: %d - Correlations positive: %d - Correlations negative: %d",nGeneratedPositiveDynamicalCorrelations+nGeneratedNegativeDynamicalCorrelations,nGeneratedPositiveDynamicalCorrelations,nGeneratedNegativeDynamicalCorrelations);
1052       Printf("Number of accepted particles: %d",gNumberOfAcceptedParticles);
1053       if(!fUseAllCharges)
1054         Printf("Pions: %lf - Kaons: %lf - Protons: %lf",1.*nGeneratedPions/nMultiplicity,1.*nGeneratedKaons/nMultiplicity,1.*nGeneratedProtons/nMultiplicity);
1055       //Printf("Calculating the balance function for %d particles",chargeVector[0]->size());
1056     }
1057
1058     fHistEventStats->Fill(4);
1059     fHistNumberOfAcceptedParticles->Fill(gNumberOfAcceptedParticles);
1060     fHistReactionPlane->Fill(fReactionPlane);
1061
1062     // Event mixing
1063     if (fRunMixing)
1064       {
1065         // 1. First get an event pool corresponding in mult (cent) and
1066         //    zvertex to the current event. Once initialized, the pool
1067         //    should contain nMix (reduced) events. This routine does not
1068         //    pre-scan the chain. The first several events of every chain
1069         //    will be skipped until the needed pools are filled to the
1070         //    specified depth. If the pool categories are not too rare, this
1071         //    should not be a problem. If they are rare, you could lose`
1072         //    statistics.
1073
1074         // 2. Collect the whole pool's content of tracks into one TObjArray
1075         //    (bgTracks), which is effectively a single background super-event.
1076
1077         // 3. The reduced and bgTracks arrays must both be passed into
1078         //    FillCorrelations(). Also nMix should be passed in, so a weight
1079         //    of 1./nMix can be applied.
1080
1081         AliEventPool* pool = fPoolMgr->GetEventPool(1., 0.,fReactionPlane);
1082
1083         if (!pool){
1084           AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", 1., 0.,fReactionPlane));
1085         }
1086         else{
1087
1088           //pool->SetDebug(1);
1089
1090           if (pool->IsReady() || pool->NTracksInPool() > 50000 / 10 || pool->GetCurrentNEvents() >= 5){
1091
1092
1093             Int_t nMix = pool->GetCurrentNEvents();
1094             //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << ", tracks in this event = "<<gNumberOfAcceptedParticles <<", event Plane = "<<fReactionPlane<<endl;
1095
1096             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
1097             //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
1098             //if (pool->IsReady())
1099             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
1100
1101             // Fill mixed-event histos here
1102             for (Int_t jMix=0; jMix<nMix; jMix++)
1103               {
1104                 TObjArray* tracksMixed = pool->GetEvent(jMix);
1105                 fMixedBalance->CalculateBalance(fReactionPlane,tracksMain,tracksMixed,1,1.,0.);
1106               }
1107           }
1108
1109           // Update the Event pool
1110           pool->UpdatePool(tracksMain);
1111           //pool->PrintInfo();
1112
1113         }//pool NULL check
1114       }//run mixing
1115
1116     //Calculate the balance function
1117     fBalance->CalculateBalance(fReactionPlane,tracksMain,NULL,1,1.,0.);
1118     //if(fRunMixing)
1119     //  fMixedBalance->CalculateBalance(fReactionPlane,tracksMixing,NULL,1,1.,0.);
1120
1121   }//event loop
1122 }
1123
1124 //________________________________________________________________________
1125 void  AliAnalysisTaskToyModel::FinishOutput() {
1126   //Printf("END BF");
1127
1128   TFile *gOutput = TFile::Open("outputToyModel.root","recreate");
1129   TDirectoryFile *dir = new TDirectoryFile("PWGCFEbyE.outputBalanceFunctionPsiAnalysis","PWGCFEbyE.outputBalanceFunctionPsiAnalysis");
1130
1131   fList->SetName("listQA");
1132   fList->SetOwner(kTRUE);
1133   dir->Add(fList);
1134   //fList->Write("listQA",TObject::kSingleKey);
1135
1136   if (!fBalance) {
1137     AliError("ERROR: fBalance not available");
1138     return;
1139   }
1140   fListBF->SetName("listBF");
1141   fListBF->SetOwner(kTRUE);
1142   dir->Add(fListBF);
1143   //fListBF->Write("listBF",TObject::kSingleKey);
1144
1145   if(fRunShuffling) {
1146     if (!fShuffledBalance) {
1147       AliError("ERROR: fShuffledBalance not available");
1148       return;
1149     }
1150     fListBFS->SetName("listBFShuffled");
1151     fListBFS->SetOwner(kTRUE);
1152     dir->Add(fListBFS);
1153     //fListBFS->Write("listBFShuffled",TObject::kSingleKey);
1154   }
1155
1156   if(fRunMixing) {
1157     if (!fMixedBalance) {
1158       AliError("ERROR: fMixedBalance not available");
1159       return;
1160     }
1161     fListBFM->SetName("listBFMixed");
1162     fListBFM->SetOwner(kTRUE);
1163     dir->Add(fListBFM);
1164     //fListBFM->Write("listBFMixed",TObject::kSingleKey);
1165   }
1166
1167   dir->Write(dir->GetName(),TObject::kSingleKey);
1168   gOutput->Close();
1169 }
1170