PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEventMixingBF.cxx

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