PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx

   1 #include "TChain.h"
   2 #include "TList.h"
   3 #include "TCanvas.h"
   4 #include "TLorentzVector.h"
   5 #include "TGraphErrors.h"
   6 #include "TH1F.h"
   7 #include "TH2F.h"
   8 #include "TArrayF.h"
   9 #include "TF1.h"
  10 #include "TRandom.h"
  11
  12 #include "AliAnalysisTaskSE.h"
  13 #include "AliAnalysisManager.h"
  14
  15 #include "AliESDVertex.h"
  16 #include "AliESDEvent.h"
  17 #include "AliESDInputHandler.h"
  18 #include "AliAODEvent.h"
  19 #include "AliAODTrack.h"
  20 #include "AliAODInputHandler.h"
  21 #include "AliGenEventHeader.h"
  22 #include "AliGenHijingEventHeader.h"
  23 #include "AliMCEventHandler.h"
  24 #include "AliMCEvent.h"
  25 #include "AliMixInputEventHandler.h"
  26 #include "AliStack.h"
  27 #include "AliESDtrackCuts.h"
  28
  29 #include "TH2D.h"
  30 #include "AliPID.h"
  31 #include "AliPIDResponse.h"
  32 #include "AliPIDCombined.h"
  33
  34 #include "AliAnalysisTaskEventMixingBF.h"
  35 #include "AliBalanceEventMixing.h"
  36
  37
  38 // Analysis task for the EventMixingBF code
  39 // Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch
  40
  41 ClassImp(AliAnalysisTaskEventMixingBF)
  42
  43 //________________________________________________________________________
  44 AliAnalysisTaskEventMixingBF::AliAnalysisTaskEventMixingBF(const char *name)
  45 : AliAnalysisTaskSE(name),
  46   fBalance(0),
  47   fRunShuffling(kFALSE),
  48   fShuffledBalance(0),
  49   fList(0),
  50   fListEventMixingBF(0),
  51   fListEventMixingBFS(0),
  52   fHistListPIDQA(0),
  53   fHistEventStats(0),
  54   fHistCentStats(0),
  55   fHistTriggerStats(0),
  56   fHistTrackStats(0),
  57   fHistVx(0),
  58   fHistVy(0),
  59   fHistVz(0),
  60   fHistClus(0),
  61   fHistDCA(0),
  62   fHistChi2(0),
  63   fHistPt(0),
  64   fHistEta(0),
  65   fHistPhi(0),
  66   fHistPhiBefore(0),
  67   fHistPhiAfter(0),
  68   fHistV0M(0),
  69   fHistRefTracks(0),
  70   fHistdEdxVsPTPCbeforePID(NULL),
  71   fHistBetavsPTOFbeforePID(NULL),
  72   fHistProbTPCvsPtbeforePID(NULL),
  73   fHistProbTOFvsPtbeforePID(NULL),
  74   fHistProbTPCTOFvsPtbeforePID(NULL),
  75   fHistNSigmaTPCvsPtbeforePID(NULL),
  76   fHistNSigmaTOFvsPtbeforePID(NULL),
  77   fHistdEdxVsPTPCafterPID(NULL),
  78   fHistBetavsPTOFafterPID(NULL),
  79   fHistProbTPCvsPtafterPID(NULL),
  80   fHistProbTOFvsPtafterPID(NULL),
  81   fHistProbTPCTOFvsPtafterPID(NULL),
  82   fHistNSigmaTPCvsPtafterPID(NULL),
  83   fHistNSigmaTOFvsPtafterPID(NULL),
  84   fPIDResponse(0x0),
  85   fPIDCombined(0x0),
  86   fParticleOfInterest(kPion),
  87   fPidDetectorConfig(kTPCTOF),
  88   fUsePID(kFALSE),
  89   fUsePIDnSigma(kTRUE),
  90   fUsePIDPropabilities(kFALSE),
  91   fPIDNSigma(3.),
  92   fMinAcceptedPIDProbability(0.8),
  93   fESDtrackCuts(0),
  94   fCentralityEstimator("V0M"),
  95   fUseCentrality(kFALSE),
  96   fCentralityPercentileMin(0.),
  97   fCentralityPercentileMax(5.),
  98   fImpactParameterMin(0.),
  99   fImpactParameterMax(20.),
 100   fUseMultiplicity(kFALSE),
 101   fNumberOfAcceptedTracksMin(0),
 102   fNumberOfAcceptedTracksMax(10000),
 103   fHistNumberOfAcceptedTracks(0),
 104   fUseOfflineTrigger(kFALSE),
 105   fVxMax(0.3),
 106   fVyMax(0.3),
 107   fVzMax(10.),
 108   nAODtrackCutBit(128),
 109   fPtMin(0.3),
 110   fPtMax(1.5),
 111   fEtaMin(-0.8),
 112   fEtaMax(-0.8),
 113   fDCAxyCut(-1),
 114   fDCAzCut(-1),
 115   fTPCchi2Cut(-1),
 116   fNClustersTPCCut(-1),
 117   fAcceptanceParameterization(0),
 118   fDifferentialV2(0),
 119   fUseFlowAfterBurner(kFALSE),
 120   fExcludeResonancesInMC(kFALSE),
 121   fUseMCPdgCode(kFALSE),
 122   fPDGCodeToBeAnalyzed(-1),
 123   fMainEvent(0x0),
 124   fMixEvent(0x0)
 125 {
 126   // Constructor
 127   // Define input and output slots here
 128   // Input slot #0 works with a TChain
 129   DefineInput(0, TChain::Class());
 130   // Output slot #0 writes into a TH1 container
 131   DefineOutput(1, TList::Class());
 132   DefineOutput(2, TList::Class());
 133   DefineOutput(3, TList::Class());
 134   DefineOutput(4, TList::Class());
 135 }
 136
 137 //________________________________________________________________________
 138 AliAnalysisTaskEventMixingBF::~AliAnalysisTaskEventMixingBF() {
 139
 140   // delete fBalance;
 141   // delete fShuffledBalance;
 142   // delete fList;
 143   // delete fListEventMixingBF;
 144   // delete fListEventMixingBFS;
 145
 146   // delete fHistEventStats;
 147   // delete fHistTrackStats;
 148   // delete fHistVx;
 149   // delete fHistVy;
 150   // delete fHistVz;
 151
 152   // delete fHistClus;
 153   // delete fHistDCA;
 154   // delete fHistChi2;
 155   // delete fHistPt;
 156   // delete fHistEta;
 157   // delete fHistPhi;
 158   // delete fHistV0M;
 159 }
 160
 161 //________________________________________________________________________
 162 void AliAnalysisTaskEventMixingBF::UserCreateOutputObjects() {
 163   // Create histograms
 164   // Called once
 165   if(!fBalance) {
 166     fBalance = new AliBalanceEventMixing();
 167     fBalance->SetAnalysisLevel("ESD");
 168     //fBalance->SetNumberOfBins(-1,16);
 169     fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
 170   }
 171   if(fRunShuffling) {
 172     if(!fShuffledBalance) {
 173       fShuffledBalance = new AliBalanceEventMixing();
 174       fShuffledBalance->SetAnalysisLevel("ESD");
 175       //fShuffledBalance->SetNumberOfBins(-1,16);
 176       fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);
 177     }
 178   }
 179
 180   //QA list
 181   fList = new TList();
 182   fList->SetName("listQA");
 183   fList->SetOwner();
 184
 185   //Balance Function list
 186   fListEventMixingBF = new TList();
 187   fListEventMixingBF->SetName("listEventMixingBF");
 188   fListEventMixingBF->SetOwner();
 189
 190   if(fRunShuffling) {
 191     fListEventMixingBFS = new TList();
 192     fListEventMixingBFS->SetName("listEventMixingBFShuffled");
 193     fListEventMixingBFS->SetOwner();
 194   }
 195
 196   //PID QA list
 197   if(fUsePID) {
 198     fHistListPIDQA = new TList();
 199     fHistListPIDQA->SetName("listQAPID");
 200     fHistListPIDQA->SetOwner();
 201   }
 202
 203   //Event stats.
 204   TString gCutName[4] = {"Total","Offline trigger",
 205                          "Vertex","Analyzed"};
 206   fHistEventStats = new TH1F("fHistEventStats",
 207                              "Event statistics;;N_{events}",
 208                              4,0.5,4.5);
 209   for(Int_t i = 1; i <= 4; i++)
 210     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
 211   fList->Add(fHistEventStats);
 212
 213   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
 214   fHistCentStats = new TH2F("fHistCentStats",
 215                              "Centrality statistics;;Cent percentile",
 216                             9,-0.5,8.5,220,-5,105);
 217   for(Int_t i = 1; i <= 9; i++)
 218     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
 219   fList->Add(fHistCentStats);
 220
 221   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
 222   fList->Add(fHistTriggerStats);
 223
 224   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);
 225   fList->Add(fHistTrackStats);
 226
 227   fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);
 228   fList->Add(fHistNumberOfAcceptedTracks);
 229
 230   // Vertex distributions
 231   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
 232   fList->Add(fHistVx);
 233   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
 234   fList->Add(fHistVy);
 235   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
 236   fList->Add(fHistVz);
 237
 238   // QA histograms
 239   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
 240   fList->Add(fHistClus);
 241   fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
 242   fList->Add(fHistChi2);
 243   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
 244   fList->Add(fHistDCA);
 245   fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);
 246   fList->Add(fHistPt);
 247   fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);
 248   fList->Add(fHistEta);
 249   fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);
 250   fList->Add(fHistPhi);
 251   fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());
 252   fList->Add(fHistPhiBefore);
 253   fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());
 254   fList->Add(fHistPhiAfter);
 255   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
 256   fList->Add(fHistV0M);
 257   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
 258   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);
 259   for(Int_t i = 1; i <= 6; i++)
 260     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
 261   fList->Add(fHistRefTracks);
 262
 263   // Balance function histograms
 264   // Initialize histograms if not done yet
 265   if(!fBalance->GetHistNp(0)){
 266     AliWarning("Histograms not yet initialized! --> Will be done now");
 267     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 268     fBalance->InitHistograms();
 269   }
 270
 271   if(fRunShuffling) {
 272     if(!fShuffledBalance->GetHistNp(0)) {
 273       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");
 274       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
 275       fShuffledBalance->InitHistograms();
 276     }
 277   }
 278
 279   for(Int_t a = 0; a < ANALYSIS_TYPES; a++){
 280     fListEventMixingBF->Add(fBalance->GetHistNp(a));
 281     fListEventMixingBF->Add(fBalance->GetHistNn(a));
 282     fListEventMixingBF->Add(fBalance->GetHistNpn(a));
 283     fListEventMixingBF->Add(fBalance->GetHistNnn(a));
 284     fListEventMixingBF->Add(fBalance->GetHistNpp(a));
 285     fListEventMixingBF->Add(fBalance->GetHistNnp(a));
 286
 287     if(fRunShuffling) {
 288       fListEventMixingBFS->Add(fShuffledBalance->GetHistNp(a));
 289       fListEventMixingBFS->Add(fShuffledBalance->GetHistNn(a));
 290       fListEventMixingBFS->Add(fShuffledBalance->GetHistNpn(a));
 291       fListEventMixingBFS->Add(fShuffledBalance->GetHistNnn(a));
 292       fListEventMixingBFS->Add(fShuffledBalance->GetHistNpp(a));
 293       fListEventMixingBFS->Add(fShuffledBalance->GetHistNnp(a));
 294     }
 295   }
 296
 297   if(fESDtrackCuts) fList->Add(fESDtrackCuts);
 298
 299   //====================PID========================//
 300   if(fUsePID) {
 301     fPIDCombined = new AliPIDCombined();
 302     fPIDCombined->SetDefaultTPCPriors();
 303
 304     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000);
 305     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition
 306
 307     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2);
 308     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition
 309
 310     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0);
 311     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition
 312
 313     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
 314     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition
 315
 316     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0);
 317     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition
 318
 319     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500);
 320     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition
 321
 322     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500);
 323     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition
 324
 325     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000);
 326     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition
 327
 328     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2);
 329     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition
 330
 331     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2);
 332     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition
 333
 334     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2);
 335     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition
 336
 337     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0);
 338     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition
 339
 340     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500);
 341     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition
 342
 343     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500);
 344     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition
 345   }
 346   //====================PID========================//
 347
 348   // Post output data.
 349   PostData(1, fList);
 350   PostData(2, fListEventMixingBF);
 351   if(fRunShuffling) PostData(3, fListEventMixingBFS);
 352   if(fUsePID) PostData(4, fHistListPIDQA);       //PID
 353 }
 354
 355 //________________________________________________________________________
 356 void AliAnalysisTaskEventMixingBF::UserExec(Option_t *) {
 357   // Main loop
 358   // Called for each event
 359   // NOTHING TO DO for event mixing!
 360 }
 361
 362 //________________________________________________________________________
 363 void  AliAnalysisTaskEventMixingBF::FinishTaskOutput(){
 364   //Printf("END EventMixingBF");
 365
 366   if (!fBalance) {
 367     AliError("ERROR: fBalance not available");
 368     return;
 369   }
 370   if(fRunShuffling) {
 371     if (!fShuffledBalance) {
 372       AliError("ERROR: fShuffledBalance not available");
 373       return;
 374     }
 375   }
 376
 377 }
 378
 379 //________________________________________________________________________
 380 void AliAnalysisTaskEventMixingBF::Terminate(Option_t *) {
 381   // Draw result to the screen
 382   // Called once at the end of the query
 383
 384   // not implemented ...
 385
 386 }
 387
 388 void AliAnalysisTaskEventMixingBF::UserExecMix(Option_t *)
 389 {
 390   // Main loop for event mixing
 391
 392   TString gAnalysisLevel = fBalance->GetAnalysisLevel();
 393
 394   AliMixInputEventHandler *mixIEH = SetupEventsForMixing();
 395
 396   Float_t fCentrality           = 0.;
 397
 398   // for HBT like cuts need magnetic field sign
 399   Float_t bSign = 0; // only used in AOD so far
 400
 401   // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E)
 402   vector<Double_t> *chargeVector[9];          // original charge
 403   for(Int_t i = 0; i < 9; i++){
 404     chargeVector[i]        = new vector<Double_t>;
 405   }
 406
 407   Double_t vCharge;
 408   Double_t vY;
 409   Double_t vEta;
 410   Double_t vPhi;
 411   Double_t vP[3];
 412   Double_t vPt;
 413   Double_t vE;
 414
 415   Int_t iMainTrackUsed = -1;
 416
 417   // -------------------------------------------------------------
 418   // At the moment MIXING only for AODs
 419   if(mixIEH){
 420
 421     //AOD analysis (vertex and track cuts also here!!!!)
 422     if(gAnalysisLevel == "AOD") {
 423       AliAODEvent* aodEventMain = dynamic_cast<AliAODEvent*>(fMainEvent);
 424       if(!aodEventMain) {
 425         AliError("ERROR: aodEventMain not available");
 426         return;
 427       }
 428       AliAODEvent *aodEventMix  = dynamic_cast<AliAODEvent *>(fMixEvent);
 429      if(!aodEventMix) {
 430         AliError("ERROR: aodEventMix not available");
 431         return;
 432       }
 433
 434      // for HBT like cuts need magnetic field sign
 435      bSign = (aodEventMain->GetMagneticField() > 0) ? 1 : -1;
 436
 437      AliAODHeader *aodHeaderMain = dynamic_cast<AliAODHeader*>(aodEventMain->GetHeader());
 438      if(!aodHeaderMain) AliFatal("Not a standard AOD");
 439
 440       // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
 441       fHistEventStats->Fill(1); //all events
 442
 443       // this is not needed (checked in mixing handler!)
 444       Bool_t isSelectedMain = kTRUE;
 445       Bool_t isSelectedMix = kTRUE;
 446
 447       if(fUseOfflineTrigger){
 448         isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 449         isSelectedMix = ((AliInputEventHandler*)((AliMultiInputEventHandler *)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetFirstMultiInputHandler())->IsEventSelected();
 450       }
 451
 452       if(isSelectedMain && isSelectedMix) {
 453         fHistEventStats->Fill(2); //triggered events
 454
 455         //Centrality stuff (centrality in AOD header)
 456         if(fUseCentrality) {
 457           fCentrality = aodHeaderMain->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 458
 459           // QA for centrality estimators
 460           fHistCentStats->Fill(0.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0M"));
 461           fHistCentStats->Fill(1.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("FMD"));
 462           fHistCentStats->Fill(2.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TRK"));
 463           fHistCentStats->Fill(3.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKL"));
 464           fHistCentStats->Fill(4.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL0"));
 465           fHistCentStats->Fill(5.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("CL1"));
 466           fHistCentStats->Fill(6.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
 467           fHistCentStats->Fill(7.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
 468           fHistCentStats->Fill(8.,aodHeaderMain->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
 469
 470           // take only events inside centrality class
 471           if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
 472             return;
 473
 474           // centrality QA (V0M)
 475           fHistV0M->Fill(aodEventMain->GetVZEROData()->GetMTotV0A(), aodEventMain->GetVZEROData()->GetMTotV0C());
 476
 477           // centrality QA (reference tracks)
 478           fHistRefTracks->Fill(0.,aodHeaderMain->GetRefMultiplicity());
 479           fHistRefTracks->Fill(1.,aodHeaderMain->GetRefMultiplicityPos());
 480           fHistRefTracks->Fill(2.,aodHeaderMain->GetRefMultiplicityNeg());
 481           fHistRefTracks->Fill(3.,aodHeaderMain->GetTPConlyRefMultiplicity());
 482           fHistRefTracks->Fill(4.,aodHeaderMain->GetNumberOfITSClusters(0));
 483           fHistRefTracks->Fill(5.,aodHeaderMain->GetNumberOfITSClusters(1));
 484           fHistRefTracks->Fill(6.,aodHeaderMain->GetNumberOfITSClusters(2));
 485           fHistRefTracks->Fill(7.,aodHeaderMain->GetNumberOfITSClusters(3));
 486           fHistRefTracks->Fill(8.,aodHeaderMain->GetNumberOfITSClusters(4));
 487         }
 488
 489         // // this is crashing (Bug in ROOT (to be solved)) but not needed (checked in mixing handler)
 490         // const AliAODVertex *vertexMain = aodEventMain->GetPrimaryVertex();
 491         // const AliAODVertex *vertexMix  = aodEventMix->GetPrimaryVertex();
 492
 493         // if(vertexMain && vertexMix) {
 494         //    Double32_t fCovMain[6];
 495         //    Double32_t fCovMix[6];
 496         //    vertexMain->GetCovarianceMatrix(fCovMain);
 497         //    vertexMix->GetCovarianceMatrix(fCovMix);
 498
 499         //    if(vertexMain->GetNContributors() > 0 && vertexMix->GetNContributors() > 0) {
 500         //     if(fCovMain[5] != 0 && fCovMix[5] != 0) {
 501         //       fHistEventStats->Fill(3); //events with a proper vertex
 502         //       if(TMath::Abs(vertexMain->GetX()) < fVxMax && TMath::Abs(vertexMix->GetX()) < fVxMax ) {
 503         //      if(TMath::Abs(vertexMain->GetY()) < fVyMax && TMath::Abs(vertexMix->GetY()) < fVyMax) {
 504         //        if(TMath::Abs(vertexMain->GetZ()) < fVzMax && TMath::Abs(vertexMix->GetZ()) < fVzMax) {
 505         //          fHistEventStats->Fill(4); //analyzed events
 506         //          fHistVx->Fill(vertexMain->GetX());
 507         //          fHistVy->Fill(vertexMain->GetY());
 508         //          fHistVz->Fill(vertexMain->GetZ());
 509
 510                     // Loop over tracks in main event
 511                     for (Int_t iTracksMain = 0; iTracksMain < aodEventMain->GetNumberOfTracks(); iTracksMain++) {
 512                       AliAODTrack* aodTrackMain = dynamic_cast<AliAODTrack *>(aodEventMain->GetTrack(iTracksMain));
 513                       if (!aodTrackMain) {
 514                         AliError(Form("ERROR: Could not receive track %d", iTracksMain));
 515                         continue;
 516                       }
 517
 518                       // AOD track cuts
 519
 520                       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
 521                       // take only TPC only tracks
 522                       fHistTrackStats->Fill(aodTrackMain->GetFilterMap());
 523                       if(!aodTrackMain->TestFilterBit(nAODtrackCutBit)) continue;
 524
 525                       vCharge = aodTrackMain->Charge();
 526                       vY      = aodTrackMain->Y();
 527                       vEta    = aodTrackMain->Eta();
 528                       vPhi    = aodTrackMain->Phi() * TMath::RadToDeg();
 529                       vE      = aodTrackMain->E();
 530                       vPt     = aodTrackMain->Pt();
 531                       aodTrackMain->PxPyPz(vP);
 532
 533                       Float_t dcaXYMain = aodTrackMain->DCA();      // this is the DCA from global track (not exactly what is cut on)
 534                       Float_t dcaZMain  = aodTrackMain->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
 535
 536
 537                       // Kinematics cuts from ESD track cuts
 538                       if( vPt < fPtMin || vPt > fPtMax)      continue;
 539                       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
 540
 541                       // Extra DCA cuts (for systematic studies [!= -1])
 542                       if( fDCAxyCut != -1 && fDCAzCut != -1){
 543                         if(TMath::Sqrt((dcaXYMain*dcaXYMain)/(fDCAxyCut*fDCAxyCut)+(dcaZMain*dcaZMain)/(fDCAzCut*fDCAzCut)) > 1 ){
 544                           continue;  // 2D cut
 545                         }
 546                       }
 547
 548                       // Extra TPC cuts (for systematic studies [!= -1])
 549                       if( fTPCchi2Cut != -1 && aodTrackMain->Chi2perNDF() > fTPCchi2Cut){
 550                         continue;
 551                       }
 552                       if( fNClustersTPCCut != -1 && aodTrackMain->GetTPCNcls() < fNClustersTPCCut){
 553                         continue;
 554                       }
 555
 556                       // fill QA histograms
 557                       fHistClus->Fill(aodTrackMain->GetITSNcls(),aodTrackMain->GetTPCNcls());
 558                       fHistDCA->Fill(dcaZMain,dcaXYMain);
 559                       fHistChi2->Fill(aodTrackMain->Chi2perNDF());
 560                       fHistPt->Fill(vPt);
 561                       fHistEta->Fill(vEta);
 562                       fHistPhi->Fill(vPhi);
 563
 564                       // fill charge vector
 565                       chargeVector[0]->push_back(vCharge);
 566                       chargeVector[1]->push_back(vY);
 567                       chargeVector[2]->push_back(vEta);
 568                       chargeVector[3]->push_back(vPhi);
 569                       chargeVector[4]->push_back(vP[0]);
 570                       chargeVector[5]->push_back(vP[1]);
 571                       chargeVector[6]->push_back(vP[2]);
 572                       chargeVector[7]->push_back(vPt);
 573                       chargeVector[8]->push_back(vE);
 574
 575                       // -------------------------------------------------------------
 576                       // for each track in main event loop over all tracks in mix event
 577                       for (Int_t iTracksMix = 0; iTracksMix < aodEventMix->GetNumberOfTracks(); iTracksMix++) {
 578
 579                         AliAODTrack* aodTrackMix = dynamic_cast<AliAODTrack *>(aodEventMix->GetTrack(iTracksMix));
 580                         if (!aodTrackMix) {
 581                           AliError(Form("ERROR: Could not receive track %d", iTracksMix));
 582                           continue;
 583                         }
 584
 585                         // AOD track cuts
 586
 587                         // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
 588                         // take only TPC only tracks
 589                         fHistTrackStats->Fill(aodTrackMix->GetFilterMap());
 590                         if(!aodTrackMix->TestFilterBit(nAODtrackCutBit)) continue;
 591
 592                         vCharge = aodTrackMix->Charge();
 593                         vY      = aodTrackMix->Y();
 594                         vEta    = aodTrackMix->Eta();
 595                         vPhi    = aodTrackMix->Phi() * TMath::RadToDeg();
 596                         vE      = aodTrackMix->E();
 597                         vPt     = aodTrackMix->Pt();
 598                         aodTrackMix->PxPyPz(vP);
 599
 600                         Float_t dcaXYMix = aodTrackMix->DCA();      // this is the DCA from global track (not exactly what is cut on)
 601                         Float_t dcaZMix  = aodTrackMix->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
 602
 603
 604                         // Kinematics cuts from ESD track cuts
 605                         if( vPt < fPtMin || vPt > fPtMax)      continue;
 606                         if( vEta < fEtaMin || vEta > fEtaMax)  continue;
 607
 608                         // Extra DCA cuts (for systematic studies [!= -1])
 609                         if( fDCAxyCut != -1 && fDCAxyCut != -1){
 610                           if(TMath::Sqrt((dcaXYMix*dcaXYMix)/(fDCAxyCut*fDCAxyCut)+(dcaZMix*dcaZMix)/(fDCAzCut*fDCAzCut)) > 1 ){
 611                             continue;  // 2D cut
 612                           }
 613                         }
 614
 615                         // Extra TPC cuts (for systematic studies [!= -1])
 616                         if( fTPCchi2Cut != -1 && aodTrackMix->Chi2perNDF() > fTPCchi2Cut){
 617                           continue;
 618                         }
 619                         if( fNClustersTPCCut != -1 && aodTrackMix->GetTPCNcls() < fNClustersTPCCut){
 620                           continue;
 621                         }
 622
 623                         // fill QA histograms
 624                         fHistClus->Fill(aodTrackMix->GetITSNcls(),aodTrackMix->GetTPCNcls());
 625                         fHistDCA->Fill(dcaZMix,dcaXYMix);
 626                         fHistChi2->Fill(aodTrackMix->Chi2perNDF());
 627                         fHistPt->Fill(vPt);
 628                         fHistEta->Fill(vEta);
 629                         fHistPhi->Fill(vPhi);
 630
 631                         // fill charge vector
 632                         chargeVector[0]->push_back(vCharge);
 633                         chargeVector[1]->push_back(vY);
 634                         chargeVector[2]->push_back(vEta);
 635                         chargeVector[3]->push_back(vPhi);
 636                         chargeVector[4]->push_back(vP[0]);
 637                         chargeVector[5]->push_back(vP[1]);
 638                         chargeVector[6]->push_back(vP[2]);
 639                         chargeVector[7]->push_back(vPt);
 640                         chargeVector[8]->push_back(vE);
 641
 642
 643
 644                       } //mix track loop
 645
 646                       // calculate balance function for each track in main event
 647                       iMainTrackUsed++; // is needed to do no double counting in Balance Function calculation
 648                       if(iMainTrackUsed >= (Int_t)chargeVector[0]->size()) break; //do not allow more tracks than in mixed event!
 649                       fBalance->CalculateBalance(fCentrality,chargeVector,iMainTrackUsed,bSign);
 650                       // clean charge vector afterwards
 651                       for(Int_t i = 0; i < 9; i++){
 652                         chargeVector[i]->clear();
 653                       }
 654
 655
 656                     } //main track loop
 657       //                  }//Vz cut
 658       //                }//Vy cut
 659       //              }//Vx cut
 660       //            }//proper vertexresolution
 661       //    }//proper number of contributors
 662       // }//vertex object valid
 663       }//triggered event
 664     }//AOD analysis
 665   }
 666 }
 667
 668 AliMixInputEventHandler *AliAnalysisTaskEventMixingBF::SetupEventsForMixing() {
 669   //sets the input handlers for event mixing
 670
 671   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 672   AliMultiInputEventHandler *inEvHMain = dynamic_cast<AliMultiInputEventHandler *>(mgr->GetInputEventHandler());
 673   if (inEvHMain) {
 674
 675       AliMixInputEventHandler *mixEH = dynamic_cast<AliMixInputEventHandler *>(inEvHMain->GetFirstMultiInputHandler());
 676       if (!mixEH) return nullptr;
 677       if (mixEH->CurrentBinIndex() < 0) {
 678          AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1");
 679          return nullptr;
 680       }
 681       AliDebug(AliLog::kDebug, Form("Mixing %lld %d [%lld,%lld] %d", mixEH->CurrentEntry(), mixEH->CurrentBinIndex(), mixEH->CurrentEntryMain(), mixEH->CurrentEntryMix(), mixEH->NumberMixed()));
 682
 683       AliInputEventHandler      *ihMainCurrent     = inEvHMain->GetFirstInputEventHandler();
 684       fMainEvent = ihMainCurrent->GetEvent();
 685
 686       AliMultiInputEventHandler *inEvHMixedCurrent = mixEH->GetFirstMultiInputHandler(); // for buffer = 1
 687       AliInputEventHandler      *ihMixedCurrent    = inEvHMixedCurrent->GetFirstInputEventHandler();
 688       fMixEvent                                    = ihMixedCurrent->GetEvent();
 689
 690       return mixEH;
 691   }
 692   return NULL;
 693 }