]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
added efficiency task for Balance Function (ESD-MC analyis), added V0 analysis for...
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2012 03:53:24 +0000 (03:53 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jul 2012 03:53:24 +0000 (03:53 +0000)
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.cxx [new file with mode: 0644]
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.h [new file with mode: 0644]
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.h
PWGCF/EBYE/macros/AddTaskBalanceEfficiency.C [new file with mode: 0644]
PWGCF/EBYE/macros/AddTaskBalanceTriggered.C

diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.cxx
new file mode 100644 (file)
index 0000000..71951a8
--- /dev/null
@@ -0,0 +1,525 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "TObjArray.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliCentrality.h"
+
+#include "AliAnalysisTaskEfficiencyBF.h"
+
+// ---------------------------------------------------------------------
+//
+// Task for calculating the efficiency of the Balance Function 
+// for single particles and pairs
+//
+// Authors: Panos Christakoglou, Michael Weber
+// 
+// ---------------------------------------------------------------------
+
+ClassImp(AliAnalysisTaskEfficiencyBF)
+
+//________________________________________________________________________
+AliAnalysisTaskEfficiencyBF::AliAnalysisTaskEfficiencyBF(const char *name) 
+  : AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0), 
+  fHistEventStats(0), fHistCentrality(0), fHistNMult(0), 
+  fHistGeneratedEtaPtPlus(0), fHistFindableEtaPtPlus(0), 
+  fHistReconstructedEtaPtPlus(0), fHistSurvivedEtaPtPlus(0),
+  fHistGeneratedEtaPtMinus(0), fHistFindableEtaPtMinus(0), 
+  fHistReconstructedEtaPtMinus(0), fHistSurvivedEtaPtMinus(0),
+  fHistGeneratedEtaPtPlusControl(0), fHistFindableEtaPtPlusControl(0), 
+  fHistReconstructedEtaPtPlusControl(0), fHistSurvivedEtaPtPlusControl(0),
+  fHistGeneratedEtaPtMinusControl(0), fHistFindableEtaPtMinusControl(0), 
+  fHistReconstructedEtaPtMinusControl(0), fHistSurvivedEtaPtMinusControl(0),
+  fHistGeneratedEtaPtPlusPlus(0), fHistFindableEtaPtPlusPlus(0), 
+  fHistReconstructedEtaPtPlusPlus(0), fHistSurvivedEtaPtPlusPlus(0),
+  fHistGeneratedEtaPtMinusMinus(0), fHistFindableEtaPtMinusMinus(0), 
+  fHistReconstructedEtaPtMinusMinus(0), fHistSurvivedEtaPtMinusMinus(0),
+  fHistGeneratedEtaPtPlusMinus(0), fHistFindableEtaPtPlusMinus(0), 
+  fHistReconstructedEtaPtPlusMinus(0), fHistSurvivedEtaPtPlusMinus(0),
+  fESDtrackCuts(0), fAnalysisMode(0), 
+  fCentralityEstimator("V0M"), fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0), 
+  fVxMax(3.0), fVyMax(3.0), fVzMax(10.), 
+  fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), fMaxDCAxy(3.0), fMaxDCAz(3.0),
+  fMinPt(0.3), fMaxPt(1.5), fMaxEta(0.8){  
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  // Output slot #0 id reserved by the base class for AOD
+  // Output slot #1 writes into a TH1 container
+  DefineOutput(1, TList::Class());
+  DefineOutput(2, TList::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEfficiencyBF::UserCreateOutputObjects() {
+  // Create histograms
+  // Called once
+
+  fQAList = new TList();
+  fQAList->SetName("QAList");
+  fQAList->SetOwner();
+
+  fOutputList = new TList();
+  fOutputList->SetName("OutputList");
+  fOutputList->SetOwner();
+
+  //Event stats.
+  TString gCutName[4] = {"Total","Offline trigger",
+                         "Vertex","Analyzed"};
+  fHistEventStats = new TH1F("fHistEventStats",
+                             "Event statistics;;N_{events}",
+                             4,0.5,4.5);
+  for(Int_t i = 1; i <= 4; i++)
+    fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
+  fQAList->Add(fHistEventStats);
+
+  //ESD analysis
+  fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
+                            20,0.5,20.5);
+  fQAList->Add(fHistCentrality);
+  
+  //multiplicity (good MC tracks)
+  TString histName;
+  histName = "fHistNMult";
+  fHistNMult = new TH1F(histName.Data(), 
+                       ";N_{mult.}",
+                       200,0,20000);
+  fQAList->Add(fHistNMult);
+  
+  //eta vs pt for MC positives
+  fHistGeneratedEtaPtPlus = new TH2F("fHistGeneratedEtaPtPlus",
+                                    "Generated positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtPlus);
+  fHistFindableEtaPtPlus = new TH2F("fHistFindableEtaPtPlus",
+                                    "Findable positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtPlus);
+  fHistReconstructedEtaPtPlus = new TH2F("fHistReconstructedEtaPtPlus",
+                                    "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtPlus);
+  fHistSurvivedEtaPtPlus = new TH2F("fHistSurvivedEtaPtPlus",
+                                    "Survived positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtPlus);
+
+  //eta vs pt for MC negatives
+  fHistGeneratedEtaPtMinus = new TH2F("fHistGeneratedEtaPtMinus",
+                                    "Generated positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtMinus);
+  fHistFindableEtaPtMinus = new TH2F("fHistFindableEtaPtMinus",
+                                    "Findable positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtMinus);
+  fHistReconstructedEtaPtMinus = new TH2F("fHistReconstructedEtaPtMinus",
+                                    "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtMinus);
+  fHistSurvivedEtaPtMinus = new TH2F("fHistSurvivedEtaPtMinus",
+                                    "Survived positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtMinus);
+
+  //eta vs pt for MC positives (control)
+  fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
+                                    "Generated positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtPlusControl);
+  fHistFindableEtaPtPlusControl = new TH2F("fHistFindableEtaPtPlusControl",
+                                    "Findable positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtPlusControl);
+  fHistReconstructedEtaPtPlusControl = new TH2F("fHistReconstructedEtaPtPlusControl",
+                                    "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtPlusControl);
+  fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
+                                    "Survived positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtPlusControl);
+
+  //eta vs pt for MC negatives (control)
+  fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
+                                    "Generated positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtMinusControl);
+  fHistFindableEtaPtMinusControl = new TH2F("fHistFindableEtaPtMinusControl",
+                                    "Findable positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtMinusControl);
+  fHistReconstructedEtaPtMinusControl = new TH2F("fHistReconstructedEtaPtMinusControl",
+                                    "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtMinusControl);
+  fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
+                                    "Survived positive primaries;#eta;p_{T} (GeV/c)",
+                                    40,-1.0,1.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtMinusControl);
+
+  //eta vs pt for MC ++
+  fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
+                                    "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
+  fHistFindableEtaPtPlusPlus = new TH2F("fHistFindableEtaPtPlusPlus",
+                                    "Findable ++ primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtPlusPlus);
+  fHistReconstructedEtaPtPlusPlus = new TH2F("fHistReconstructedEtaPtPlusPlus",
+                                    "Reconstructed ++ primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtPlusPlus);
+  fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
+                                    "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
+
+  //eta vs pt for MC --
+  fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
+                                    "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
+  fHistFindableEtaPtMinusMinus = new TH2F("fHistFindableEtaPtMinusMinus",
+                                    "Findable -- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtMinusMinus);
+  fHistReconstructedEtaPtMinusMinus = new TH2F("fHistReconstructedEtaPtMinusMinus",
+                                    "Reconstructed -- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtMinusMinus);
+  fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
+                                    "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
+
+  //eta vs pt for MC +-
+  fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
+                                    "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
+  fHistFindableEtaPtPlusMinus = new TH2F("fHistFindableEtaPtPlusMinus",
+                                    "Findable +- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistFindableEtaPtPlusMinus);
+  fHistReconstructedEtaPtPlusMinus = new TH2F("fHistReconstructedEtaPtPlusMinus",
+                                    "Reconstructed +- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistReconstructedEtaPtPlusMinus);
+  fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
+                                    "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
+                                    40,0.0,2.0,49,0.1,5.0);
+  fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
+
+  fQAList->Print();
+  fOutputList->Print();
+
+  PostData(1, fQAList);
+  PostData(2, fOutputList);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEfficiencyBF::UserExec(Option_t *) {
+  // Main loop
+  // Called for each event
+
+  // Post output data.
+  //ESD analysis
+  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!fESD) {
+    printf("ERROR: fESD not available\n");
+    return;
+  }
+  
+  AliMCEvent* mcEvent = MCEvent();
+  if (!mcEvent) {
+    Printf("ERROR: Could not retrieve MC event");
+    return;
+  }
+  AliStack* stack = mcEvent->Stack();
+  if (!stack) {
+    Printf("ERROR: Could not retrieve MC stack");
+    return;
+  }
+
+  // arrays for 2 particle histograms
+  Int_t nMCLabelCounter         = 0;
+  const Int_t maxMCLabelCounter = 20000;
+
+  Double_t pt[maxMCLabelCounter];
+  Double_t eta[maxMCLabelCounter];
+  Int_t level[maxMCLabelCounter];
+  Int_t charge[maxMCLabelCounter];
+
+
+  //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks()));
+  fHistEventStats->Fill(1); //all events
+    
+  //Centrality stuff
+  AliCentrality *centrality = fESD->GetCentrality();
+  Int_t nCentrality = 0;
+  nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.);
+
+  //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
+
+  if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
+                                         fCentralityPercentileMax,
+                                         fCentralityEstimator.Data())) {
+    fHistEventStats->Fill(2); //triggered + centrality
+    fHistCentrality->Fill(nCentrality+1);
+
+    //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
+  
+    if(fAnalysisMode.CompareTo("TPC") == 0 ) {
+      const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
+      if(vertex) {
+       if(vertex->GetNContributors() > 0) {
+         if(vertex->GetZRes() != 0) {
+           fHistEventStats->Fill(3); //events with a proper vertex
+           if(TMath::Abs(vertex->GetXv()) < fVxMax) {
+             if(TMath::Abs(vertex->GetYv()) < fVyMax) {
+               if(TMath::Abs(vertex->GetZv()) < fVzMax) {
+                 fHistEventStats->Fill(4); //analyzed events
+                 
+                 Int_t nMCParticles = mcEvent->GetNumberOfTracks();
+                 TArrayI labelMCArray(nMCParticles);
+                 
+                 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
+                   AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks);
+                   if (!mcTrack) {
+                     Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
+                     continue;
+                   }
+                   
+                   //exclude particles generated out of the acceptance
+                   Double_t vz = mcTrack->Zv();
+                   if (TMath::Abs(vz) > 50.) continue;
+                   
+                   //acceptance
+                   if(TMath::Abs(mcTrack->Eta()) > 2.5) 
+                     continue;
+                   if((mcTrack->Pt() > 5.0)||(mcTrack->Pt() < 0.1)) 
+                     continue;
+                   
+                   TParticle* particle = mcTrack->Particle();
+                   if(!particle) continue;
+                   if(!stack->IsPhysicalPrimary(iTracks)) continue;
+
+                   if(iTracks <= stack->GetNprimary()) {                     
+                     Short_t gMCCharge = mcTrack->Charge();
+                     
+                     if(gMCCharge > 0)
+                       fHistGeneratedEtaPtPlus->Fill(particle->Eta(),
+                                                     particle->Pt());
+                     else if(gMCCharge < 0)
+                       fHistGeneratedEtaPtMinus->Fill(particle->Eta(),
+                                                      particle->Pt());
+
+                     
+                     // findable tracks --> DOES NOT WORK????
+                     // Loop over Track References
+                     Bool_t labelTPC = kTRUE;
+                     AliTrackReference* trackRef = 0;
+
+                     for (Int_t iTrackRef = 0; iTrackRef < mcTrack->GetNumberOfTrackReferences(); iTrackRef++) {
+                       trackRef = mcTrack->GetTrackReference(iTrackRef);
+                       if(trackRef) {
+                         Int_t detectorId = trackRef->DetectorId();
+                         if (detectorId == AliTrackReference::kTPC) {
+                           labelTPC = kTRUE;
+                           break;
+                         }
+                       }
+                     }//loop over track references
+
+                     if(labelTPC) {
+                       labelMCArray.AddAt(iTracks,nMCLabelCounter);
+
+                       if(nMCLabelCounter >= maxMCLabelCounter){
+                         AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
+                         break;
+                       }
+
+                       // fill the arrays for 2 particle analysis
+                       eta[nMCLabelCounter]    = particle->Eta();
+                       pt[nMCLabelCounter]     = particle->Pt();
+                       charge[nMCLabelCounter] = gMCCharge;
+                       // findable = generated in this case!
+                       level[nMCLabelCounter]  = 1;
+                       
+                       nMCLabelCounter += 1;
+                       
+                       if(gMCCharge > 0)
+                         fHistFindableEtaPtPlus->Fill(particle->Eta(),
+                                                      particle->Pt());
+                       else if(gMCCharge < 0)
+                         fHistFindableEtaPtMinus->Fill(particle->Eta(),
+                                                       particle->Pt());
+                     }
+                   }//primaries
+                 }//loop over MC particles
+               
+                 fHistNMult->Fill(nMCLabelCounter);
+                 
+                 // not used so far
+                 //Float_t dcaXY = 0.0, dcaZ = 0.0;
+
+                 //ESD track loop
+                 Int_t nGoodTracks = fESD->GetNumberOfTracks();
+                 
+                 TArrayI labelArray(nGoodTracks);
+                 Int_t labelCounter = 0;
+                 for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
+                   AliESDtrack* track = fESD->GetTrack(iTracks);
+                   //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks);
+                   if(!track) continue;
+
+                   AliESDtrack *tpcOnlyTrack = new AliESDtrack();
+                   
+                   if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
+                     delete tpcOnlyTrack;
+                     continue;
+                   }
+
+                   Int_t label = TMath::Abs(track->GetTPCLabel());
+                   if(IsLabelUsed(labelArray,label)) continue;
+                   labelArray.AddAt(label,labelCounter);
+                   labelCounter += 1;
+
+                   Bool_t iFound = kFALSE;
+                   Int_t mcGoods = nMCLabelCounter;
+                   for (Int_t k = 0; k < mcGoods; k++) {
+                     Int_t mcLabel = labelMCArray.At(k);
+                     iFound = kFALSE;
+                                             
+                     if (mcLabel != TMath::Abs(label)) continue;
+                     if(mcLabel != label) continue;
+                     if(label > stack->GetNtrack()) continue;
+                     
+                     TParticle *particle = stack->Particle(label);
+                     if(!particle) continue;
+                     
+                     //acceptance
+                     if(TMath::Abs(particle->Eta()) > 2.5) 
+                       continue;
+                     if((particle->Pt() > 5.0)||(particle->Pt() < 0.1)) 
+                       continue;
+                     if(!stack->IsPhysicalPrimary(label)) continue;
+                     
+                     if(label <= stack->GetNprimary()) {
+                       
+                       // reconstructed
+                       level[k]  = 2;
+                       
+                       Short_t gCharge = track->Charge();
+                       if(gCharge > 0)
+                         fHistReconstructedEtaPtPlus->Fill(particle->Eta(),
+                                                           particle->Pt());
+                       else if(gCharge < 0)
+                         fHistReconstructedEtaPtMinus->Fill(particle->Eta(),
+                                                            particle->Pt());
+                       
+                       // track cuts + analysis kinematic cuts
+                       if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
+
+                         // survived
+                         level[k]  = 3;
+
+                         if(gCharge > 0)
+                           fHistSurvivedEtaPtPlus->Fill(particle->Eta(),
+                                                        particle->Pt());
+                         else if(gCharge < 0)
+                           fHistSurvivedEtaPtMinus->Fill(particle->Eta(),
+                                                         particle->Pt());
+                       }//track cuts
+                     }//primary particles
+                   }//findable track loop
+                 }//ESD track loop
+               
+                 labelMCArray.Reset();
+                 labelArray.Reset();
+
+               }//Vz cut
+             }//Vy cut
+           }//Vx cut
+         }//Vz resolution
+       }//number of contributors
+      }//valid vertex
+    }//TPC analysis mode
+  }//centrality  
+      
+
+
+  // Here comes the 2 particle analysis
+  // loop over all good MC particles
+  for (Int_t i = 0; i < nMCLabelCounter ; i++) {
+
+    // control 1D histograms (charge might be different?)
+    if(charge[i] > 0){
+      if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
+      if(level[i] > 1) fHistReconstructedEtaPtPlusControl->Fill(eta[i],pt[i]);
+      if(level[i] > 2) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
+    }
+    else if(charge[i] < 0){
+      if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
+      if(level[i] > 1) fHistReconstructedEtaPtMinusControl->Fill(eta[i],pt[i]);
+      if(level[i] > 2) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
+    }
+
+
+    for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
+      
+      if(charge[i] > 0 && charge[j] > 0 ){
+       if(level[i] > 0 && level[j] > 0) fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 1 && level[j] > 1) fHistReconstructedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 2 && level[j] > 2) fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+      }
+      
+      else if(charge[i] < 0 && charge[j] < 0 ){
+       if(level[i] > 0 && level[j] > 0) fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 1 && level[j] > 1) fHistReconstructedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 2 && level[j] > 2) fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+      }
+      else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
+       if(level[i] > 0 && level[j] > 0) fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 1 && level[j] > 1) fHistReconstructedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+       if(level[i] > 2 && level[j] > 2) fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
+      }
+    }
+  }      
+}
+//________________________________________________________________________
+void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) {
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutputList) {
+    printf("ERROR: Output list not available\n");
+    return;
+  }
+}
+
+//____________________________________________________________________//
+Bool_t AliAnalysisTaskEfficiencyBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
+  //Checks if the label is used already
+  Bool_t status = kFALSE;
+  for(Int_t i = 0; i < labelArray.GetSize(); i++) {
+    if(labelArray.At(i) == label)
+      status = kTRUE;
+  }
+
+  return status;
+}
diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.h b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.h
new file mode 100644 (file)
index 0000000..75a7b5b
--- /dev/null
@@ -0,0 +1,162 @@
+#ifndef ALIANALYSISTASKEFFICIENCYBF_cxx
+#define ALIANALYSISTASKEFFICIENCYBF_cxx
+
+// ---------------------------------------------------------------------
+//
+// Task for calculating the efficiency of the Balance Function 
+// for single particles and pairs
+//
+// Authors: Panos Christakoglou, Michael Weber
+// 
+// ---------------------------------------------------------------------
+
+class TH1F;
+class TH2F;
+class TString;
+class AliESDEvent;
+class AliESDtrackCuts;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskEfficiencyBF : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskEfficiencyBF() : AliAnalysisTaskSE(), 
+    fESD(0), fQAList(0), fOutputList(0), 
+    fHistEventStats(0), fHistCentrality(0), fHistNMult(0), 
+    fHistGeneratedEtaPtPlus(0), fHistFindableEtaPtPlus(0), 
+    fHistReconstructedEtaPtPlus(0), fHistSurvivedEtaPtPlus(0),
+    fHistGeneratedEtaPtMinus(0), fHistFindableEtaPtMinus(0), 
+    fHistReconstructedEtaPtMinus(0), fHistSurvivedEtaPtMinus(0),
+    fHistGeneratedEtaPtPlusControl(0), fHistFindableEtaPtPlusControl(0), 
+    fHistReconstructedEtaPtPlusControl(0), fHistSurvivedEtaPtPlusControl(0),
+    fHistGeneratedEtaPtMinusControl(0), fHistFindableEtaPtMinusControl(0), 
+    fHistReconstructedEtaPtMinusControl(0), fHistSurvivedEtaPtMinusControl(0),
+    fHistGeneratedEtaPtPlusPlus(0), fHistFindableEtaPtPlusPlus(0), 
+    fHistReconstructedEtaPtPlusPlus(0), fHistSurvivedEtaPtPlusPlus(0),
+    fHistGeneratedEtaPtMinusMinus(0), fHistFindableEtaPtMinusMinus(0), 
+    fHistReconstructedEtaPtMinusMinus(0), fHistSurvivedEtaPtMinusMinus(0),
+    fHistGeneratedEtaPtPlusMinus(0), fHistFindableEtaPtPlusMinus(0), 
+    fHistReconstructedEtaPtPlusMinus(0), fHistSurvivedEtaPtPlusMinus(0),
+    fESDtrackCuts(0), fAnalysisMode(0), 
+    fCentralityEstimator("V0M"), fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0), 
+    fVxMax(3.0), fVyMax(3.0), fVzMax(10.), 
+    fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), fMaxDCAxy(3.0), fMaxDCAz(3.0),
+    fMinPt(0.3), fMaxPt(1.5), fMaxEta(0.8) {}
+  AliAnalysisTaskEfficiencyBF(const char *name);
+  virtual ~AliAnalysisTaskEfficiencyBF() {}
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  
+  Bool_t   IsLabelUsed(TArrayI array, Int_t label);
+
+  void SetAnalysisCutObject(AliESDtrackCuts *const trackCuts) {
+    fESDtrackCuts = trackCuts;}
+  void SetVertexDiamond(Double_t vx, Double_t vy, Double_t vz) {
+    fVxMax = vx;
+    fVyMax = vy;
+    fVzMax = vz;
+  }
+  //Centrality
+  void SetCentralityEstimator(const char* centralityEstimator) {
+    fCentralityEstimator = centralityEstimator;}
+  void SetCentralityPercentileRange(Float_t min, Float_t max) { 
+    fCentralityPercentileMin=min;
+    fCentralityPercentileMax=max;
+  }
+
+  void SetAnalysisMode(const char* analysisMode) {
+    fAnalysisMode = analysisMode;}
+
+  //Track cuts
+  void SetMinNumberOfTPCClusters(Double_t min) {
+    fMinNumberOfTPCClusters = min;}
+  void SetMaxChi2PerTPCCluster(Double_t max) {
+    fMaxChi2PerTPCCluster = max;}
+  void SetMaxDCAxy(Double_t max) {
+    fMaxDCAxy = max;}
+  void SetMaxDCAz(Double_t max) {
+    fMaxDCAz = max;}
+  void SetMinPt(Double_t minPt) {
+    fMinPt = minPt;}
+  void SetMaxPt(Double_t maxPt) {
+    fMaxPt = maxPt;}
+  void SetMaxEta(Double_t maxEta) {
+    fMaxEta = maxEta;}
+
+  
+ private:
+  AliESDEvent *fESD;    //! ESD object
+  TList       *fQAList; //! QA list
+  TList       *fOutputList; //! Output list
+  
+  // QA histograms
+  TH1F        *fHistEventStats; //!event stats
+  TH1F        *fHistCentrality; //!centrality
+  TH1F        *fHistNMult; //! nmult   
+
+  // output histograms (single particles)
+  TH2F        *fHistGeneratedEtaPtPlus;//!correction map for positives (generated)
+  TH2F        *fHistFindableEtaPtPlus;//!correction map for positives (findable)
+  TH2F        *fHistReconstructedEtaPtPlus;//!correction map for positives (reconstructed)
+  TH2F        *fHistSurvivedEtaPtPlus;//!correction map positives (survived)
+
+  TH2F        *fHistGeneratedEtaPtMinus;//!correction map for negatives (generated)
+  TH2F        *fHistFindableEtaPtMinus;//!correction map for negatives (findable)
+  TH2F        *fHistReconstructedEtaPtMinus;//!correction map for negatives (reconstructed)
+  TH2F        *fHistSurvivedEtaPtMinus;//!correction map negatives (survived)
+
+  TH2F        *fHistGeneratedEtaPtPlusControl;//!correction map for positives (generated)
+  TH2F        *fHistFindableEtaPtPlusControl;//!correction map for positives (findable)
+  TH2F        *fHistReconstructedEtaPtPlusControl;//!correction map for positives (reconstructed)
+  TH2F        *fHistSurvivedEtaPtPlusControl;//!correction map positives (survived)
+
+  TH2F        *fHistGeneratedEtaPtMinusControl;//!correction map for negatives (generated)
+  TH2F        *fHistFindableEtaPtMinusControl;//!correction map for negatives (findable)
+  TH2F        *fHistReconstructedEtaPtMinusControl;//!correction map for negatives (reconstructed)
+  TH2F        *fHistSurvivedEtaPtMinusControl;//!correction map negatives (survived)
+
+  // output histograms (pairs)
+  TH2F        *fHistGeneratedEtaPtPlusPlus;//!correction map for ++ (generated)
+  TH2F        *fHistFindableEtaPtPlusPlus;//!correction map for ++ (findable)
+  TH2F        *fHistReconstructedEtaPtPlusPlus;//!correction map for ++ (reconstructed)
+  TH2F        *fHistSurvivedEtaPtPlusPlus;//!correction map ++ (survived)
+
+  TH2F        *fHistGeneratedEtaPtMinusMinus;//!correction map for -- (generated)
+  TH2F        *fHistFindableEtaPtMinusMinus;//!correction map for -- (findable)
+  TH2F        *fHistReconstructedEtaPtMinusMinus;//!correction map for -- (reconstructed)
+  TH2F        *fHistSurvivedEtaPtMinusMinus;//!correction map -- (survived)
+
+  TH2F        *fHistGeneratedEtaPtPlusMinus;//!correction map for +- (generated)
+  TH2F        *fHistFindableEtaPtPlusMinus;//!correction map for +- (findable)
+  TH2F        *fHistReconstructedEtaPtPlusMinus;//!correction map for +- (reconstructed)
+  TH2F        *fHistSurvivedEtaPtPlusMinus;//!correction map +- (survived)
+
+
+
+  AliESDtrackCuts *fESDtrackCuts; //ESD track cuts
+
+  TString fAnalysisMode;//"TPC", "Global"
+
+  TString fCentralityEstimator;//"V0M","TRK","TKL","ZDC","FMD"
+  Float_t fCentralityPercentileMin, fCentralityPercentileMax; //min-max centrality percentile
+
+  Double_t fVxMax;//vxmax
+  Double_t fVyMax;//vymax
+  Double_t fVzMax;//vzmax
+
+  Double_t fMinNumberOfTPCClusters; //
+  Double_t fMaxChi2PerTPCCluster; //
+  Double_t fMaxDCAxy, fMaxDCAz;//
+  Double_t fMinPt, fMaxPt;
+  Double_t fMaxEta;
+
+  AliAnalysisTaskEfficiencyBF(const AliAnalysisTaskEfficiencyBF&); // not implemented
+  AliAnalysisTaskEfficiencyBF& operator=(const AliAnalysisTaskEfficiencyBF&); // not implemented
+  
+  ClassDef(AliAnalysisTaskEfficiencyBF, 1); // example of analysis
+};
+
+#endif
index a6c9bf0b72b3dd60d0d44525d77f4355383059f8..f428ada076508fe2ed3e5a9d23b768eefaf4c2cd 100755 (executable)
 // Analysis task for the TriggeredBF code\r
 // Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
 \r
+// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
+//\r
+// For the V0 part:\r
+// --> AliAnalysisTaskExtractV0AOD (by david.chinellato@gmail.com)\r
+//\r
+// +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
+\r
 using std::cout;\r
 using std::endl;\r
 using std::vector;\r
@@ -56,11 +63,15 @@ AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name)
   fMixingTracks(50000),\r
   fMixedBalance(0),\r
   fPoolMgr(0),\r
+  fRunV0(kFALSE),\r
+  fPIDResponse(0),\r
+  fPIDCombined(0),\r
   fList(0),\r
   fListTriggeredBF(0),\r
   fListTriggeredBFS(0),\r
   fListTriggeredBFM(0),\r
   fHistListPIDQA(0),\r
+  fHistListV0(0),\r
   fHistEventStats(0),\r
   fHistCentStats(0),\r
   fHistTriggerStats(0),\r
@@ -78,6 +89,24 @@ AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name)
   fHistPhiAfter(0),\r
   fHistV0M(0),\r
   fHistRefTracks(0),\r
+  fHistV0MultiplicityBeforeTrigSel(0),\r
+  fHistV0MultiplicityForTrigEvt(0),\r
+  fHistV0MultiplicityForSelEvt(0),\r
+  fHistV0MultiplicityForSelEvtNoTPCOnly(0),\r
+  fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),\r
+  fHistMultiplicityBeforeTrigSel(0),\r
+  fHistMultiplicityForTrigEvt(0),\r
+  fHistMultiplicity(0),\r
+  fHistMultiplicityNoTPCOnly(0),\r
+  fHistMultiplicityNoTPCOnlyNoPileup(0),\r
+  fHistV0InvMassK0(0),\r
+  fHistV0InvMassLambda(0),\r
+  fHistV0InvMassAntiLambda(0),\r
+  fHistV0Armenteros(0),\r
+  fHistV0SelInvMassK0(0),\r
+  fHistV0SelInvMassLambda(0),\r
+  fHistV0SelInvMassAntiLambda(0),\r
+  fHistV0SelArmenteros(0),\r
   fCentralityEstimator("V0M"),\r
   fUseCentrality(kFALSE),\r
   fCentralityPercentileMin(0.), \r
@@ -111,6 +140,7 @@ AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name)
   DefineOutput(2, TList::Class());\r
   DefineOutput(3, TList::Class());\r
   DefineOutput(4, TList::Class());\r
+  DefineOutput(5, TList::Class());\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -223,8 +253,140 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
   fList->Add(fHistRefTracks);\r
 \r
+  //------------------------------------------------\r
+  // V0 Multiplicity Histograms\r
+  //------------------------------------------------\r
+  if(fRunV0){\r
+    fHistListV0 = new TList();\r
+    fHistListV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner\r
+    \r
+    if(! fHistV0MultiplicityBeforeTrigSel) {\r
+      fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", \r
+                                                 "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", \r
+                                                 25, 0, 25);\r
+      fHistListV0->Add(fHistV0MultiplicityBeforeTrigSel);\r
+    }\r
+    \r
+    if(! fHistV0MultiplicityForTrigEvt) {\r
+      fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", \r
+                                              "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", \r
+                                              25, 0, 25);\r
+      fHistListV0->Add(fHistV0MultiplicityForTrigEvt);\r
+    }\r
+    \r
+    if(! fHistV0MultiplicityForSelEvt) {\r
+      fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", \r
+                                             "V0s per event;Nbr of V0s/Evt;Events", \r
+                                             25, 0, 25);\r
+      fHistListV0->Add(fHistV0MultiplicityForSelEvt);\r
+    }\r
+    \r
+    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {\r
+    fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", \r
+                                                    "V0s per event;Nbr of V0s/Evt;Events", \r
+                                                    25, 0, 25);\r
+    fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);\r
+    }\r
+    \r
+    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {\r
+      fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", \r
+                                                            "V0s per event;Nbr of V0s/Evt;Events", \r
+                                                              25, 0, 25);\r
+      fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);\r
+    }\r
+    \r
+    //------------------------------------------------\r
+    // Track Multiplicity Histograms\r
+    //------------------------------------------------\r
+    \r
+    if(! fHistMultiplicityBeforeTrigSel) {\r
+      fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", \r
+                                             "Tracks per event;Nbr of Tracks;Events", \r
+                                               200, 0, 200);           \r
+      fHistListV0->Add(fHistMultiplicityBeforeTrigSel);\r
+    }\r
+    if(! fHistMultiplicityForTrigEvt) {\r
+      fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", \r
+                                            "Tracks per event;Nbr of Tracks;Events", \r
+                                            200, 0, 200);              \r
+      fHistListV0->Add(fHistMultiplicityForTrigEvt);\r
+    }\r
+    if(! fHistMultiplicity) {\r
+      fHistMultiplicity = new TH1F("fHistMultiplicity", \r
+                                  "Tracks per event;Nbr of Tracks;Events", \r
+                                  200, 0, 200);                \r
+      fHistListV0->Add(fHistMultiplicity);\r
+    }\r
+    if(! fHistMultiplicityNoTPCOnly) {\r
+      fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", \r
+                                           "Tracks per event;Nbr of Tracks;Events", \r
+                                           200, 0, 200);               \r
+    fHistListV0->Add(fHistMultiplicityNoTPCOnly);\r
+    }\r
+    if(! fHistMultiplicityNoTPCOnlyNoPileup) {\r
+      fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", \r
+                                                   "Tracks per event;Nbr of Tracks;Events", \r
+                                                   200, 0, 200);               \r
+      fHistListV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);\r
+    }\r
 \r
-\r
+    //------------------------------------------------\r
+    // V0 selection Histograms (before)\r
+    //------------------------------------------------\r
+    if(!fHistV0InvMassK0) {\r
+      fHistV0InvMassK0 = new TH1F("fHistV0InvMassK0",\r
+                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0InvMassK0);\r
+    }\r
+    if(!fHistV0InvMassLambda) {\r
+      fHistV0InvMassLambda = new TH1F("fHistV0InvMassLambda",\r
+                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0InvMassLambda);\r
+    }\r
+    if(!fHistV0InvMassAntiLambda) {\r
+      fHistV0InvMassAntiLambda = new TH1F("fHistV0InvMassAntiLambda",\r
+                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0InvMassAntiLambda);\r
+    }\r
+    if(!fHistV0Armenteros) {\r
+      fHistV0Armenteros = new TH2F("fHistV0Armenteros",\r
+                                 "Armenteros plot;#alpha;q_{t}",\r
+                                  200,-1,1,200,0,0.5);\r
+      fHistListV0->Add(fHistV0Armenteros);\r
+    }\r
+    \r
+    //------------------------------------------------\r
+    // V0 selection Histograms (after)\r
+    //------------------------------------------------\r
+    if(!fHistV0SelInvMassK0) {\r
+      fHistV0SelInvMassK0 = new TH1F("fHistV0SelInvMassK0",\r
+                                 "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0SelInvMassK0);\r
+    }\r
+    if(!fHistV0SelInvMassLambda) {\r
+      fHistV0SelInvMassLambda = new TH1F("fHistV0SelInvMassLambda",\r
+                                 "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0SelInvMassLambda);\r
+    }\r
+    if(!fHistV0SelInvMassAntiLambda) {\r
+      fHistV0SelInvMassAntiLambda = new TH1F("fHistV0SelInvMassAntiLambda",\r
+                                 "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
+                                 200,0,2);\r
+      fHistListV0->Add(fHistV0SelInvMassAntiLambda);\r
+    }\r
+    if(!fHistV0SelArmenteros) {\r
+      fHistV0SelArmenteros = new TH2F("fHistV0SelArmenteros",\r
+                                 "Armenteros plot;#alpha;q_{t}",\r
+                                  200,-1,1,200,0,0.5);\r
+      fHistListV0->Add(fHistV0SelArmenteros);\r
+    }\r
+  }//V0\r
+    \r
   // Balance function histograms\r
   // Initialize histograms if not done yet\r
   if(!fBalance->GetHistNp()){\r
@@ -274,6 +436,11 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
     fListTriggeredBFM->Add(fMixedBalance->GetHistNnp());\r
   }  \r
 \r
+  // PID Response task active?\r
+  if(fRunV0) {\r
+    fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
+    if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
+  }\r
 \r
   // Event Mixing\r
   Int_t trackDepth = fMixingTracks; \r
@@ -296,6 +463,7 @@ void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {
   PostData(2, fListTriggeredBF);\r
   if(fRunShuffling) PostData(3, fListTriggeredBFS);\r
   if(fRunMixing) PostData(4, fListTriggeredBFM);\r
+  if(fRunV0) PostData(5,fHistListV0);\r
 }\r
 \r
 //________________________________________________________________________\r
@@ -321,7 +489,9 @@ void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {
     }\r
     \r
     // get the accepted tracks in main event\r
-    TObjArray *tracksMain = GetAcceptedTracks(eventMain);\r
+    TObjArray *tracksMain = NULL;\r
+    if(fRunV0) tracksMain = GetAcceptedV0s(eventMain);\r
+    else       tracksMain = GetAcceptedTracks(eventMain);\r
 \r
     // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
     TObjArray* tracksShuffled = NULL;\r
@@ -407,11 +577,18 @@ Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){
 \r
   Bool_t isSelectedMain = kTRUE;\r
   Float_t fCentrality = -1.;\r
+  Int_t nV0s          = event->GetNumberOfV0s();\r
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
   \r
   if(fUseOfflineTrigger)\r
     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
   \r
+  //V0 QA histograms (before trigger selection)\r
+  if(fRunV0){\r
+    fHistMultiplicityBeforeTrigSel->Fill ( -1 );\r
+    fHistV0MultiplicityBeforeTrigSel->Fill ( nV0s );\r
+  }\r
+  \r
   if(isSelectedMain) {\r
     fHistEventStats->Fill(2); //triggered events\r
     \r
@@ -445,6 +622,12 @@ Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){
        fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
        fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
        fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
+\r
+       //V0 QA histograms (after trigger selection)\r
+       if(fRunV0){\r
+         fHistMultiplicityForTrigEvt->Fill ( fCentrality );\r
+         fHistV0MultiplicityForTrigEvt->Fill ( nV0s );\r
+       }\r
       }\r
     }\r
     \r
@@ -465,6 +648,31 @@ Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){
                fHistVy->Fill(vertex->GetY());\r
                fHistVz->Fill(vertex->GetZ());\r
 \r
+\r
+\r
+               //V0 QA histograms (vertex Z check)\r
+               if(fRunV0){\r
+                 fHistV0MultiplicityForSelEvt ->Fill( nV0s );\r
+                 fHistMultiplicity->Fill(fCentrality);\r
+\r
+                 //V0 QA histograms (Only look at events with well-established PV)\r
+                 const AliAODVertex *lPrimarySPDVtx = ((AliAODEvent*)event)->GetPrimaryVertexSPD();\r
+                 if(lPrimarySPDVtx){\r
+                   fHistMultiplicityNoTPCOnly->Fill ( fCentrality );\r
+                   fHistV0MultiplicityForSelEvtNoTPCOnly->Fill ( nV0s );\r
+                   \r
+                   //V0 QA histograms (Pileup Rejection)\r
+                   // FIXME : quality selection regarding pile-up rejection \r
+                   fHistMultiplicityNoTPCOnlyNoPileup->Fill(fCentrality);\r
+                   fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( nV0s );\r
+\r
+                 }\r
+                 else{\r
+                   return -1;\r
+                 }\r
+               }\r
+               \r
+               \r
                // take only events inside centrality class\r
                if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){\r
                  return fCentrality;           \r
@@ -553,6 +761,289 @@ TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedTracks(AliVEvent *event){
   return tracksAccepted;\r
 }\r
 \r
+//________________________________________________________________________\r
+TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedV0s(AliVEvent *event){\r
+  // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
+  // Fills QA histograms\r
+\r
+  //output TObjArray holding all good tracks\r
+  TObjArray* tracksAccepted = new TObjArray;\r
+  tracksAccepted->SetOwner(kTRUE);\r
+\r
+  Double_t v_charge;\r
+  Double_t v_eta;\r
+  Double_t v_phi;\r
+  Double_t v_pt;\r
+  \r
+  //------------------------------------------------\r
+  // MAIN LAMBDA LOOP STARTS HERE (basically a copy of AliAnalysisTaskExtractV0AOD)\r
+  //------------------------------------------------\r
+\r
+  // parameters (for the time being hard coded here) --> from David for EbyE Lambdas\r
+  Bool_t fkUseOnTheFly = kFALSE;\r
+  Double_t fRapidityBoundary  = 0.5; \r
+  Double_t fCutDaughterEta    = 0.8;\r
+  Double_t fCutV0Radius       = 0.9;\r
+  Double_t fCutDCANegToPV     = 0.1;\r
+  Double_t fCutDCAPosToPV     = 0.1;\r
+  Double_t fCutDCAV0Daughters = 1.0;\r
+  Double_t fCutV0CosPA        = 0.9995;\r
+  Double_t fMassLambda        = 1.115683;\r
+  Double_t fCutMassLambda     = 0.007;\r
+  Double_t fCutProperLifetime = 3*7.9;\r
+  Double_t fCutLeastNumberOfCrossedRows = 70;\r
+  Double_t fCutLeastNumberOfCrossedRowsOverFindable = 0.8;\r
+  Double_t fCutTPCPIDNSigmasProton  = 3.0;\r
+  Double_t fCutTPCPIDNSigmasPion    = 5.0;\r
+\r
+\r
+  //Variable definition\r
+  Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;\r
+  Double_t lChi2V0 = 0;\r
+  Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;\r
+  Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;\r
+  Double_t lV0CosineOfPointingAngle = 0;\r
+  Double_t lV0Radius = 0, lPt = 0;\r
+  Double_t lEta = 0, lPhi = 0;\r
+  Double_t lRap = 0, lRapK0Short = 0, lRapLambda = 0;\r
+  Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;\r
+  Double_t lAlphaV0 = 0, lPtArmV0 = 0;\r
+  \r
+  Double_t fMinV0Pt = 0; \r
+  Double_t fMaxV0Pt = 100; \r
+  \r
+\r
+  \r
+  // some event observables\r
+  Int_t nv0s = event->GetNumberOfV0s();\r
+  Double_t tPrimaryVtxPosition[3];\r
+  const AliVVertex *primaryVtx = event->GetPrimaryVertex();\r
+  tPrimaryVtxPosition[0] = primaryVtx->GetX();\r
+  tPrimaryVtxPosition[1] = primaryVtx->GetY();\r
+  tPrimaryVtxPosition[2] = primaryVtx->GetZ();\r
+\r
+\r
+  //loop over V0s  \r
+  for (Int_t iV0 = 0; iV0 < nv0s; iV0++) \r
+    {// This is the begining of the V0 loop\r
+      AliAODv0 *v0 = ((AliAODEvent*)event)->GetV0(iV0);\r
+      if (!v0) continue;\r
+\r
+      //Obsolete at AOD level... \r
+      //---> Fix On-the-Fly candidates, count how many swapped\r
+      //if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){\r
+      //  fHistSwappedV0Counter -> Fill( 1 );\r
+      //}else{\r
+      //  fHistSwappedV0Counter -> Fill( 0 ); \r
+      //}\r
+      //if ( fkUseOnTheFly ) CheckChargeV0(v0); \r
+      \r
+      Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0); \r
+      Double_t tV0mom[3];\r
+      v0->GetPxPyPz( tV0mom ); \r
+      Double_t lV0TotalMomentum = TMath::Sqrt(\r
+                                             tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );\r
+      \r
+      lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);\r
+      lPt = v0->Pt();\r
+      lEta = v0->Eta();\r
+      lPhi = v0->Phi()*TMath::RadToDeg();\r
+      lRapK0Short = v0->RapK0Short();\r
+      lRapLambda  = v0->RapLambda();\r
+      lRap        = lRapLambda;//v0->Y(); //FIXME!!!\r
+      if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;\r
+      \r
+      //UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPosID());\r
+      //UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetPosID());\r
+\r
+      Double_t lMomPos[3]; //v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);\r
+      Double_t lMomNeg[3]; //v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);\r
+      lMomPos[0] = v0->MomPosX();\r
+      lMomPos[1] = v0->MomPosY();\r
+      lMomPos[2] = v0->MomPosZ();\r
+      lMomNeg[0] = v0->MomNegX();\r
+      lMomNeg[1] = v0->MomNegY();\r
+      lMomNeg[2] = v0->MomNegZ();\r
+      \r
+      AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter\r
+      AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter\r
+      if (!pTrack || !nTrack) {\r
+       Printf("ERROR: Could not retreive one of the daughter track");\r
+       continue;\r
+      }\r
+\r
+      //Daughter Eta for Eta selection, afterwards\r
+      Double_t lNegEta = nTrack->Eta();\r
+      Double_t lPosEta = pTrack->Eta();\r
+      \r
+      // Filter like-sign V0 (next: add counter and distribution)\r
+      if ( pTrack->Charge() == nTrack->Charge()){\r
+       continue;\r
+      } \r
+      \r
+      //Quick test this far! \r
+      \r
+\r
+      //________________________________________________________________________\r
+      // Track quality cuts \r
+      Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);\r
+      Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);\r
+      Float_t lLeastNbrCrossedRows =  (lPosTrackCrossedRows>lNegTrackCrossedRows) ? lNegTrackCrossedRows : lPosTrackCrossedRows;\r
+\r
+      // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)\r
+      if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
+      if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
+      \r
+      if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;\r
+      \r
+      //Findable clusters > 0 condition\r
+      if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;\r
+      \r
+      //Compute ratio Crossed Rows / Findable clusters\r
+      //Note: above test avoids division by zero! \r
+      Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); \r
+      Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); \r
+      Float_t lLeastNbrCrossedRowsOverFindable = (lPosTrackCrossedRowsOverFindable>lNegTrackCrossedRowsOverFindable) ? lNegTrackCrossedRowsOverFindable : lPosTrackCrossedRowsOverFindable;\r
+\r
+      //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here\r
+      if ( lLeastNbrCrossedRowsOverFindable < 0.8) continue;\r
+      \r
+      //End track Quality Cuts\r
+      //________________________________________________________________________\r
+      \r
+      \r
+      lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();\r
+      lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();\r
+          \r
+      lOnFlyStatus = v0->GetOnFlyStatus();\r
+      lChi2V0 = v0->Chi2V0();\r
+      lDcaV0Daughters = v0->DcaV0Daughters();\r
+      lDcaV0ToPrimVertex = v0->DcaV0ToPrimVertex();\r
+      lV0CosineOfPointingAngle = v0->CosPointingAngle(tPrimaryVtxPosition);\r
+      \r
+      // Distance over total momentum\r
+      Double_t lDistOverTotMom = TMath::Sqrt(\r
+                                   TMath::Power( tDecayVertexV0[0] - tPrimaryVtxPosition[0] , 2) +\r
+                                   TMath::Power( tDecayVertexV0[1] - tPrimaryVtxPosition[1] , 2) +\r
+                                   TMath::Power( tDecayVertexV0[2] - tPrimaryVtxPosition[2] , 2)\r
+                                   );\r
+      lDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure\r
+      \r
+      \r
+      // Getting invariant mass infos directly from ESD\r
+      lInvMassK0s        = v0->MassK0Short();\r
+      lInvMassLambda     = v0->MassLambda();\r
+      lInvMassAntiLambda = v0->MassAntiLambda();\r
+      lAlphaV0 = v0->AlphaV0();\r
+      lPtArmV0 = v0->PtArmV0();\r
+\r
+      //Official means of acquiring N-sigmas \r
+      Double_t lNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );\r
+      Double_t lNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );\r
+      Double_t lNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );\r
+      Double_t lNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );\r
+\r
+      //V0 QA histograms (before V0 selection)\r
+      fHistV0InvMassK0->Fill(lInvMassK0s);\r
+      fHistV0InvMassLambda->Fill(lInvMassLambda);\r
+      fHistV0InvMassAntiLambda->Fill(lInvMassAntiLambda);\r
+      fHistV0Armenteros->Fill(lAlphaV0,lPtArmV0);\r
+      \r
+      \r
+      //First Selection: Reject OnFly\r
+      if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){\r
+       \r
+\r
+       //Second Selection: rough 20-sigma band, parametric.    \r
+       //K0Short: Enough to parametrize peak broadening with linear function.    \r
+       Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*lPt; \r
+       Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*lPt;\r
+       \r
+       //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)\r
+       //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)\r
+       Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*lPt); \r
+       Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*lPt);\r
+       \r
+       //Do Selection      \r
+       if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda     ) || \r
+           (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda     ) || \r
+           (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short    ) ){\r
+\r
+\r
+      //         //Pre-selection in case this is AA...\r
+      //         //if( fkIsNuclear == kFALSE ) fTree->Fill();\r
+      //         //if( fkIsNuclear == kTRUE){ \r
+      //         //If this is a nuclear collision___________________\r
+      //         // ... pre-filter with TPC, daughter eta selection\r
+\r
+         \r
+         if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda \r
+              && TMath::Abs(lNSigmasPosProton) < 6.0 && TMath::Abs(lNSigmasNegPion) < 6.0 ) || \r
+             (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda \r
+              && TMath::Abs(lNSigmasNegProton) < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ||  \r
+             (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short \r
+              && TMath::Abs(lNSigmasNegPion)   < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ){\r
+           \r
+           //insane test\r
+           if ( TMath::Abs(lNegEta)<0.8 && TMath::Abs(lPosEta)<0.8 ){\r
+\r
+             // start the fine selection (usually done in post processing, but we don't have time to waste) --> Lambdas!\r
+             if(\r
+                TMath::Abs(lRap)<fRapidityBoundary &&\r
+                TMath::Abs(lNegEta)       <= fCutDaughterEta               &&                   \r
+                TMath::Abs(lPosEta)       <= fCutDaughterEta               &&\r
+                lV0Radius                 >= fCutV0Radius                  &&\r
+                lDcaNegToPrimVertex       >= fCutDCANegToPV                &&\r
+                lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&\r
+                lDcaV0Daughters           <= fCutDCAV0Daughters            &&\r
+                lV0CosineOfPointingAngle  >= fCutV0CosPA                   && \r
+                fMassLambda*lDistOverTotMom    <= fCutProperLifetime       &&\r
+                lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&\r
+                lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&\r
+                lPtArmV0 * 5 < TMath::Abs(lAlphaV0)                        && \r
+                ((TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmasPion     &&\r
+                 TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmasProton) ||\r
+                 (TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmasPion     &&\r
+                  TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmasProton))            \r
+                )\r
+               {\r
+\r
+                 //V0 QA histograms (after V0 selection)\r
+                 fHistV0SelInvMassK0->Fill(lInvMassK0s);\r
+                 fHistV0SelInvMassLambda->Fill(lInvMassLambda);\r
+                 fHistV0SelInvMassAntiLambda->Fill(lInvMassAntiLambda);\r
+\r
+                 // this means a V0 candidate is found\r
+                 if(TMath::Abs(lInvMassLambda-fMassLambda) < fCutMassLambda ||\r
+                    TMath::Abs(lInvMassAntiLambda-fMassLambda) < fCutMassLambda){\r
+\r
+                   fHistV0SelArmenteros->Fill(lAlphaV0,lPtArmV0);                \r
+\r
+                   v_eta    = lEta;\r
+                   v_phi    = lPhi;\r
+                   v_pt     = lPt;\r
+                   if(lAlphaV0 > 0) v_charge = 1;\r
+                   if(lAlphaV0 < 0) v_charge = -1;\r
+\r
+                   // fill QA histograms\r
+                   fHistPt->Fill(v_pt);\r
+                   fHistEta->Fill(v_eta);\r
+                   fHistPhi->Fill(v_phi);\r
+                   \r
+                   // add the track to the TObjArray\r
+                   tracksAccepted->Add(new AliBFBasicParticle(v_eta, v_phi, v_pt, v_charge));\r
+                 }\r
+               }\r
+             }\r
+         }\r
+         //}//end nuclear_____________________________________\r
+       }\r
+      }\r
+    }//V0 loop\r
+  \r
+  return tracksAccepted;\r
+}\r
+\r
 //________________________________________________________________________\r
 TObjArray* AliAnalysisTaskTriggeredBF::GetShuffledTracks(TObjArray *tracks){\r
   // Clones TObjArray and returns it with tracks after shuffling the charges\r
index a335d0db1f64688f99b8c9bcfb45adc87bc7ddd6..db5a6a81dcf476027db757872b5db301b511dfac 100755 (executable)
@@ -7,6 +7,8 @@
 #include "AliLog.h"\r
 #include "AliAnalysisTaskSE.h"\r
 #include "AliBalanceTriggered.h"\r
+#include "AliPIDResponse.h"\r
+#include "AliPIDCombined.h"\r
 \r
 class TList;\r
 class TH1F;\r
@@ -41,6 +43,7 @@ class AliAnalysisTaskTriggeredBF : public AliAnalysisTaskSE {
   }\r
   void SetMixingTracks(Int_t tracks) { fMixingTracks = tracks; }\r
 \r
+  void SetRunV0(Bool_t runV0 = kTRUE) { fRunV0 = runV0; }\r
  \r
   void SetVertexDiamond(Double_t vx, Double_t vy, Double_t vz) {\r
     fVxMax = vx;\r
@@ -97,6 +100,7 @@ class AliAnalysisTaskTriggeredBF : public AliAnalysisTaskSE {
  private:\r
   Float_t    IsEventAccepted(AliVEvent* event);\r
   TObjArray* GetAcceptedTracks(AliVEvent* event);\r
+  TObjArray* GetAcceptedV0s(AliVEvent* event);\r
   TObjArray* GetShuffledTracks(TObjArray* tracks);\r
 \r
   AliBalanceTriggered *fBalance; //TriggeredBF object\r
@@ -106,13 +110,17 @@ class AliAnalysisTaskTriggeredBF : public AliAnalysisTaskSE {
   Int_t  fMixingTracks;\r
   AliBalanceTriggered *fMixedBalance; //TriggeredBF object (mixed)\r
   AliEventPoolManager*     fPoolMgr;         //! event pool manager\r
-    \r
+  Bool_t fRunV0;\r
+\r
+  AliPIDResponse       *fPIDResponse;     //! PID response object\r
+  AliPIDCombined       *fPIDCombined;     //! combined PID object\r
 \r
   TList *fList; //fList object\r
   TList *fListTriggeredBF; //fList object\r
   TList *fListTriggeredBFS; //fList object (shuffling)\r
   TList *fListTriggeredBFM; //fList object (mixing)\r
   TList *fHistListPIDQA;  //! list of histograms\r
+  TList *fHistListV0; // list of V0 histograms\r
 \r
   TH1F *fHistEventStats; //event stats\r
   TH2F *fHistCentStats; //centrality stats\r
@@ -133,6 +141,30 @@ class AliAnalysisTaskTriggeredBF : public AliAnalysisTaskSE {
   TH2F *fHistV0M;//\r
   TH2F *fHistRefTracks;//\r
 \r
+  // V0 histograms\r
+  TH1F    *fHistV0MultiplicityBeforeTrigSel;             //! V0 multiplicity distribution\r
+  TH1F    *fHistV0MultiplicityForTrigEvt;                //! V0 multiplicity distribution\r
+  TH1F    *fHistV0MultiplicityForSelEvt;                 //! V0 multiplicity distribution\r
+  TH1F    *fHistV0MultiplicityForSelEvtNoTPCOnly;        //! V0 multiplicity distribution\r
+  TH1F    *fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup;//! V0 multiplicity distribution\r
+  \r
+  TH1F    *fHistMultiplicityBeforeTrigSel;             //! multiplicity distribution    \r
+  TH1F    *fHistMultiplicityForTrigEvt;                        //! multiplicity distribution\r
+  TH1F    *fHistMultiplicity;                                                  //! multiplicity distribution\r
+  TH1F    *fHistMultiplicityNoTPCOnly;                         //! multiplicity distribution\r
+  TH1F    *fHistMultiplicityNoTPCOnlyNoPileup;                 //! multiplicity distribution\r
+\r
+  //before selection\r
+  TH1F* fHistV0InvMassK0;                           // Invariant mass K0\r
+  TH1F* fHistV0InvMassLambda;                       // Invariant mass Lambda\r
+  TH1F* fHistV0InvMassAntiLambda;                   // Invariant mass AntiLambda\r
+  TH2F* fHistV0Armenteros;                          // Armenteros plot\r
+\r
+  //after selection\r
+  TH1F* fHistV0SelInvMassK0;                           // Invariant mass K0\r
+  TH1F* fHistV0SelInvMassLambda;                       // Invariant mass Lambda\r
+  TH1F* fHistV0SelInvMassAntiLambda;                   // Invariant mass AntiLambda\r
+  TH2F* fHistV0SelArmenteros;                          // Armenteros plot\r
  \r
   TString fCentralityEstimator;      //"V0M","TRK","TKL","ZDC","FMD"\r
   Bool_t fUseCentrality;//use the centrality (PbPb) or not (pp)\r
diff --git a/PWGCF/EBYE/macros/AddTaskBalanceEfficiency.C b/PWGCF/EBYE/macros/AddTaskBalanceEfficiency.C
new file mode 100644 (file)
index 0000000..eab0075
--- /dev/null
@@ -0,0 +1,82 @@
+//_________________________________________________________//\r
+AliAnalysisTaskEfficiencyBF *AddTaskBalanceEfficiency(\r
+                                                     TString  centralityEstimator="V0M",\r
+                                                     Double_t centrMin=0.,\r
+                                                     Double_t centrMax=80.,\r
+                                                     Double_t vertexZ=10.,\r
+                                                     TString fileNameBase="AnalysisResults"\r
+                                                    ) {\r
+\r
+  // Creates a balance function analysis task and adds it to the analysis manager.\r
+  // Get the pointer to the existing analysis manager via the static access method.\r
+  TString centralityName("");\r
+\r
+  centralityName+=Form("%s_%.0f-%.0f_%.0f",centralityEstimator.Data(),centrMin,centrMax,vertexZ);\r
+  TString outputFileName(fileNameBase);\r
+  outputFileName.Append(".root");\r
+  \r
+  //===========================================================================\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskTriggeredBF", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }\r
+\r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //===========================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskTriggeredBF", "This task requires an input event handler");\r
+    return NULL;\r
+  }\r
+  TString analysisType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"\r
+  if(dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())) analysisType = "MC";\r
+\r
\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+\r
+\r
+  AliAnalysisTaskEfficiencyBF *taskEfficiencyBF = new AliAnalysisTaskEfficiencyBF("TaskEfficiencyBF");\r
+\r
+  // analysis mode\r
+  taskEfficiencyBF->SetAnalysisMode("TPC");\r
+\r
+  // centrality\r
+  taskEfficiencyBF->SetCentralityEstimator(centralityEstimator);\r
+  taskEfficiencyBF->SetCentralityPercentileRange(centrMin,centrMax);\r
+  taskEfficiencyBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+\r
+  // vertex\r
+  taskEfficiencyBF->SetVertexDiamond(.3,.3,vertexZ);\r
+\r
+  // track cuts\r
+  AliESDtrackCuts *cuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();\r
+  cuts->SetRequireTPCStandAlone(kTRUE); // TPC only cuts!  \r
+  taskEfficiencyBF->SetAnalysisCutObject(cuts);\r
+\r
+  // analysis kinematic cuts\r
+  taskEfficiencyBF->SetMinPt(0.3);\r
+  taskEfficiencyBF->SetMaxPt(1.5);\r
+  taskEfficiencyBF->SetMaxEta(0.8);\r
+\r
+\r
+  // ADD the task\r
+  //===========================================================================\r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskEfficiencyBF);\r
+\r
+  \r
+  // Create ONLY the output containers for the data produced by the task.\r
+  // Get and connect other common input/output containers via the manager as below\r
+  //==============================================================================\r
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();\r
+  outputFileName += ":PWGCFEbyE.outputBalanceFunctionEfficiencyAnalysis";\r
+  AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQA_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutEfficiencyBF = mgr->CreateContainer(Form("listEfficiencyBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
\r
+  mgr->ConnectInput(taskEfficiencyBF, 0, mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(taskEfficiencyBF, 1, coutQA);\r
+  mgr->ConnectOutput(taskEfficiencyBF, 2, coutEfficiencyBF);\r
+\r
+  return taskEfficiencyBF;\r
+}\r
index 5d2070d5583f90d22a9b8cedb5b2fbbbe64381d6..903f73c905b5aecb897ee3a852d31dc61a0c64df 100644 (file)
@@ -11,6 +11,7 @@ AliAnalysisTaskTriggeredBF *AddTaskBalanceTriggered(Double_t centrMin=0.,
                                                 Double_t centrMax=80.,\r
                                                 Bool_t gRunShuffling=kFALSE,\r
                                                 Bool_t gRunMixing=kFALSE,\r
+                                                Bool_t gRunV0=kFALSE,\r
                                                 TString centralityEstimator="V0M",\r
                                                 Double_t vertexZ=10.,\r
                                                 Double_t DCAxy=-1,\r
@@ -55,6 +56,7 @@ AliAnalysisTaskTriggeredBF *AddTaskBalanceTriggered(Double_t centrMin=0.,
   centralityName+="_Bit";\r
   centralityName+=Form("%d",AODfilterBit);\r
   if(bCentralTrigger)   centralityName+="_withCentralTrigger";\r
+  if(gRunV0)   centralityName+="_V0";\r
 \r
 \r
 \r
@@ -113,6 +115,9 @@ AliAnalysisTaskTriggeredBF *AddTaskBalanceTriggered(Double_t centrMin=0.,
     taskTriggeredBF->SetMixingObject(bfm);\r
     taskTriggeredBF->SetMixingTracks(50000);\r
   }\r
+  if(gRunV0){\r
+    taskTriggeredBF->SetRunV0(kTRUE);\r
+  }\r
 \r
   taskTriggeredBF->SetCentralityPercentileRange(centrMin,centrMax);\r
   if(analysisType == "AOD") {\r
@@ -153,12 +158,14 @@ AliAnalysisTaskTriggeredBF *AddTaskBalanceTriggered(Double_t centrMin=0.,
   AliAnalysisDataContainer *coutTriggeredBF = mgr->CreateContainer(Form("listTriggeredBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
   if(gRunShuffling) AliAnalysisDataContainer *coutTriggeredBFS = mgr->CreateContainer(Form("listTriggeredBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
   if(gRunMixing) AliAnalysisDataContainer *coutTriggeredBFM = mgr->CreateContainer(Form("listTriggeredBFMixed_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunV0) AliAnalysisDataContainer *coutQAV0 = mgr->CreateContainer(Form("listQAV0_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
 \r
   mgr->ConnectInput(taskTriggeredBF, 0, mgr->GetCommonInputContainer());\r
   mgr->ConnectOutput(taskTriggeredBF, 1, coutQA);\r
   mgr->ConnectOutput(taskTriggeredBF, 2, coutTriggeredBF);\r
   if(gRunShuffling) mgr->ConnectOutput(taskTriggeredBF, 3, coutTriggeredBFS);\r
   if(gRunMixing) mgr->ConnectOutput(taskTriggeredBF, 4, coutTriggeredBFM);\r
+  if(gRunV0) mgr->ConnectOutput(taskTriggeredBF, 5, coutQAV0);\r
 \r
   return taskTriggeredBF;\r
 }\r