--- /dev/null
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3D.h"
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "TObjArray.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliTHn.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliCentrality.h"
+#include "AliGenHijingEventHeader.h"
+
+#include "AliAnalysisTaskEfficiencyBFPsi.h"
+
+// ---------------------------------------------------------------------
+//
+// Task for calculating the efficiency of the Balance Function
+// for single particles and pairs
+//
+// Authors: Panos Christakoglou, Michael Weber
+//
+// ---------------------------------------------------------------------
+
+ClassImp(AliAnalysisTaskEfficiencyBFPsi)
+
+//________________________________________________________________________
+ AliAnalysisTaskEfficiencyBFPsi::AliAnalysisTaskEfficiencyBFPsi(const char *name):
+ AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0),
+ fHistEventStats(0), fHistCentrality(0), fHistNMult(0),
+ fHistGeneratedPlus(0), fHistSurvivedPlus(0),
+ fHistGeneratedMinus(0), fHistSurvivedMinus(0),
+ fHistGeneratedPlusPlus(0), fHistSurvivedPlusPlus(0),
+ fHistGeneratedMinusMinus(0), fHistSurvivedMinusMinus(0),
+ fHistGeneratedPlusMinus(0), fHistSurvivedPlusMinus(0),
+ fHistGeneratedMinusPlus(0), fHistSurvivedMinusPlus(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), fEtaRangeMax(0.8),
+ fPtRangeMin(0.1), fPtRangeMax(5.0), fPhiRangeMin(0.0),fPhiRangeMax(6.28) {
+ // 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 AliAnalysisTaskEfficiencyBFPsi::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);
+
+ //Setting up the AliTHn
+ Int_t anaSteps = 1; // analysis steps
+ Int_t iBinSingle[kVariablesSingle]; // binning for track variables
+ Double_t* dBinsSingle[kVariablesSingle]; // bins for track variables
+ TString axisTitleSingle[kVariablesSingle]; // axis titles for track variables
+
+ // two particle histograms
+ Int_t iBinPair[kVariablesPair]; // binning for track variables
+ Double_t* dBinsPair[kVariablesPair]; // bins for track variables
+ TString axisTitlePair[kVariablesPair]; // axis titles for track variables
+
+ //Psi_2: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (all)
+ const Int_t kNPsi2Bins = 4;
+ Double_t psi2Bins[kNPsi2Bins+1] = {-0.5,0.5,1.5,2.5,3.5};
+ iBinSingle[0] = kNPsi2Bins;
+ dBinsSingle[0] = psi2Bins;
+ axisTitleSingle[0] = "#phi - #Psi_{2} (a.u.)";
+ iBinPair[0] = kNPsi2Bins;
+ dBinsPair[0] = psi2Bins;
+ axisTitlePair[0] = "#phi - #Psi_{2} (a.u.)";
+
+ // delta eta
+ const Int_t kNDeltaEtaBins = 80;
+ Double_t deltaEtaBins[kNDeltaEtaBins+1];
+ for(Int_t i = 0; i < kNDeltaEtaBins+1; i++)
+ deltaEtaBins[i] = -2.0 + i * 0.05;
+ iBinPair[1] = kNDeltaEtaBins;
+ dBinsPair[1] = deltaEtaBins;
+ axisTitlePair[1] = "#Delta #eta";
+
+ // delta phi
+ const Int_t kNDeltaPhiBins = 72;
+ Double_t deltaPhiBins[kNDeltaPhiBins+1];
+ for(Int_t i = 0; i < kNDeltaPhiBins+1; i++){
+ deltaPhiBins[i] = -180.0 + i * 5.;
+ }
+ iBinPair[2] = kNDeltaPhiBins;
+ dBinsPair[2] = deltaPhiBins;
+ axisTitlePair[2] = "#Delta #phi (#circ)";
+
+ // pt(trigger-associated)
+ const Int_t kNPtBins = 16;
+ Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
+ iBinSingle[1] = kNPtBins;
+ dBinsSingle[1] = ptBins;
+ axisTitleSingle[1] = "p_{t}^{trig.} (GeV/c)";
+
+ iBinPair[3] = kNPtBins;
+ dBinsPair[3] = ptBins;
+ axisTitlePair[3] = "p_{t}^{trig.} (GeV/c)";
+
+ iBinPair[4] = kNPtBins;
+ dBinsPair[4] = ptBins;
+ axisTitlePair[4] = "p_{t}^{assoc.} (GeV/c)";
+
+ //=============================//
+ //Generated: Single particle distributions
+ fHistGeneratedPlus = new AliTHn("fHistGeneratedPlus",
+ "Generated positive primaries",
+ anaSteps,kVariablesSingle,iBinSingle);
+ for (Int_t j = 0; j < kVariablesSingle; j++) {
+ fHistGeneratedPlus->SetBinLimits(j, dBinsSingle[j]);
+ fHistGeneratedPlus->SetVarTitle(j, axisTitleSingle[j]);
+ }
+ fOutputList->Add(fHistGeneratedPlus);
+
+ fHistGeneratedMinus = new AliTHn("fHistGeneratedMinus",
+ "Generated positive primaries",
+ anaSteps,kVariablesSingle,iBinSingle);
+ for (Int_t j = 0; j < kVariablesSingle; j++) {
+ fHistGeneratedMinus->SetBinLimits(j, dBinsSingle[j]);
+ fHistGeneratedMinus->SetVarTitle(j, axisTitleSingle[j]);
+ }
+ fOutputList->Add(fHistGeneratedMinus);
+
+ //Survived: Single particle distributions
+ fHistSurvivedPlus = new AliTHn("fHistSurvivedPlus",
+ "Survived positive primaries",
+ anaSteps,kVariablesSingle,iBinSingle);
+ for (Int_t j = 0; j < kVariablesSingle; j++) {
+ fHistSurvivedPlus->SetBinLimits(j, dBinsSingle[j]);
+ fHistSurvivedPlus->SetVarTitle(j, axisTitleSingle[j]);
+ }
+ fOutputList->Add(fHistSurvivedPlus);
+
+ fHistSurvivedMinus = new AliTHn("fHistSurvivedMinus",
+ "Survived positive primaries",
+ anaSteps,kVariablesSingle,iBinSingle);
+ for (Int_t j = 0; j < kVariablesSingle; j++) {
+ fHistSurvivedMinus->SetBinLimits(j, dBinsSingle[j]);
+ fHistSurvivedMinus->SetVarTitle(j, axisTitleSingle[j]);
+ }
+ fOutputList->Add(fHistSurvivedMinus);
+
+ //=============================//
+ //Generated-Survived: Particle pairs +-
+ fHistGeneratedPlusMinus = new AliTHn("fHistGeneratedPlusMinus",
+ "Generated +- primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistGeneratedPlusMinus->SetBinLimits(j, dBinsPair[j]);
+ fHistGeneratedPlusMinus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistGeneratedPlusMinus);
+
+ fHistSurvivedPlusMinus = new AliTHn("fHistSurvivedPlusMinus",
+ "Survived +- primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistSurvivedPlusMinus->SetBinLimits(j, dBinsPair[j]);
+ fHistSurvivedPlusMinus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistSurvivedPlusMinus);
+
+ //Generated-Survived: Particle pairs -+
+ fHistGeneratedMinusPlus = new AliTHn("fHistGeneratedMinusPlus",
+ "Generated -+ primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistGeneratedMinusPlus->SetBinLimits(j, dBinsPair[j]);
+ fHistGeneratedMinusPlus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistGeneratedMinusPlus);
+
+ fHistSurvivedMinusPlus = new AliTHn("fHistSurvivedMinusPlus",
+ "Survived -+ primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistSurvivedMinusPlus->SetBinLimits(j, dBinsPair[j]);
+ fHistSurvivedMinusPlus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistSurvivedMinusPlus);
+
+ //Generated-Survived: Particle pairs ++
+ fHistGeneratedPlusPlus = new AliTHn("fHistGeneratedPlusPlus",
+ "Generated ++ primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistGeneratedPlusPlus->SetBinLimits(j, dBinsPair[j]);
+ fHistGeneratedPlusPlus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistGeneratedPlusPlus);
+
+ fHistSurvivedPlusPlus = new AliTHn("fHistSurvivedPlusPlus",
+ "Survived ++ primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistSurvivedPlusPlus->SetBinLimits(j, dBinsPair[j]);
+ fHistSurvivedPlusPlus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistSurvivedPlusPlus);
+
+ //Generated-Survived: Particle pairs --
+ fHistGeneratedMinusMinus = new AliTHn("fHistGeneratedMinusMinus",
+ "Generated -- primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistGeneratedMinusMinus->SetBinLimits(j, dBinsPair[j]);
+ fHistGeneratedMinusMinus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistGeneratedMinusMinus);
+
+ fHistSurvivedMinusMinus = new AliTHn("fHistSurvivedMinusMinus",
+ "Survived -- primaries",
+ anaSteps,kVariablesPair,iBinPair);
+ for (Int_t j = 0; j < kVariablesPair; j++) {
+ fHistSurvivedMinusMinus->SetBinLimits(j, dBinsPair[j]);
+ fHistSurvivedMinusMinus->SetVarTitle(j, axisTitlePair[j]);
+ }
+ fOutputList->Add(fHistSurvivedMinusMinus);
+ //=============================//
+
+ fQAList->Print();
+ fOutputList->Print();
+
+ PostData(1, fQAList);
+ PostData(2, fOutputList);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEfficiencyBFPsi::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 phiMinusPsi[maxMCLabelCounter];
+ Double_t eta[maxMCLabelCounter];
+ Double_t pt[maxMCLabelCounter];
+ Double_t phi[maxMCLabelCounter];
+ Int_t level[maxMCLabelCounter];
+ Int_t charge[maxMCLabelCounter];
+
+ Double_t trackVariablesSingle[kVariablesSingle];
+ Double_t trackVariablesPair[kVariablesPair];
+
+ Double_t gReactionPlane = 0.;
+ AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(dynamic_cast<AliMCEvent*>(mcEvent)->GenEventHeader());
+ if (headerH) {
+ gReactionPlane = headerH->ReactionPlaneAngle();
+ gReactionPlane *= TMath::RadToDeg();
+ }
+
+ //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()) > fEtaRangeMax)
+ continue;
+ if((mcTrack->Pt() > fPtRangeMax)||(mcTrack->Pt() < fPtRangeMin))
+ continue;
+ if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin))
+ continue;
+
+ TParticle* particle = mcTrack->Particle();
+ if(!particle) continue;
+ if(!stack->IsPhysicalPrimary(iTracks)) continue;
+
+ if(iTracks <= stack->GetNprimary()) {
+ Short_t gMCCharge = mcTrack->Charge();
+ Float_t firstPhi = mcTrack->Phi()*TMath::RadToDeg();
+ Float_t firstPt = mcTrack->Pt();
+
+ // Event plane (determine psi bin)
+ Double_t gPsiMinusPhi = 0.;
+ Double_t gPsiMinusPhiBin = -10.;
+ gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
+ //in-plane
+ if((gPsiMinusPhi <= 7.5)||
+ ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
+ gPsiMinusPhiBin = 0.0;
+ //intermediate
+ else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
+ ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
+ ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
+ ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
+ gPsiMinusPhiBin = 1.0;
+ //out of plane
+ else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
+ ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
+ gPsiMinusPhiBin = 2.0;
+ //everything else
+ else
+ gPsiMinusPhiBin = 3.0;
+
+ trackVariablesSingle[0] = gPsiMinusPhiBin;
+ trackVariablesSingle[1] = firstPt;
+
+ if(gMCCharge > 0)
+ fHistGeneratedPlus->Fill(trackVariablesSingle,0,1.);
+ else if(gMCCharge < 0)
+ fHistGeneratedMinus->Fill(trackVariablesSingle,0,1.);
+
+ 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
+ phiMinusPsi[nMCLabelCounter] = gPsiMinusPhiBin;
+ eta[nMCLabelCounter] = particle->Eta();
+ pt[nMCLabelCounter] = particle->Pt();
+ phi[nMCLabelCounter] = particle->Phi();
+ charge[nMCLabelCounter] = gMCCharge;
+ // findable = generated in this case!
+
+ level[nMCLabelCounter] = 1;
+ nMCLabelCounter += 1;
+ }//primaries
+ }//loop over MC particles
+
+ fHistNMult->Fill(nMCLabelCounter);
+
+ //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()) > fEtaRangeMax)
+ continue;
+ if((particle->Pt() > fPtRangeMax)||(particle->Pt() < fPtRangeMin))
+ continue;
+ if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin))
+ continue;
+
+ if(!stack->IsPhysicalPrimary(label)) continue;
+
+ if(label <= stack->GetNprimary()) {
+ Short_t gCharge = track->Charge();
+ // track cuts + analysis kinematic cuts
+ if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
+ // survived
+ level[k] = 2;
+
+ Float_t firstPhi = particle->Phi()*TMath::RadToDeg();
+ Float_t firstPt = particle->Pt();
+
+ // Event plane (determine psi bin)
+ Double_t gPsiMinusPhi = 0.;
+ Double_t gPsiMinusPhiBin = -10.;
+ gPsiMinusPhi = TMath::Abs(firstPhi - gReactionPlane);
+ //in-plane
+ if((gPsiMinusPhi <= 7.5)||
+ ((172.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5)))
+ gPsiMinusPhiBin = 0.0;
+ //intermediate
+ else if(((37.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5))||
+ ((127.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5))||
+ ((217.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5))||
+ ((307.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5)))
+ gPsiMinusPhiBin = 1.0;
+ //out of plane
+ else if(((82.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5))||
+ ((262.5 <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5)))
+ gPsiMinusPhiBin = 2.0;
+ //everything else
+ else
+ gPsiMinusPhiBin = 3.0;
+
+ trackVariablesSingle[0] = gPsiMinusPhiBin;
+ trackVariablesSingle[1] = firstPt;
+
+ if(gCharge > 0)
+ fHistSurvivedPlus->Fill(trackVariablesSingle,0,1.);
+ else if(gCharge < 0)
+ fHistSurvivedMinus->Fill(trackVariablesSingle,0,1.);
+ }//track cuts
+ }//primary particles
+ }//loop over the generated
+ }//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++) {
+ Float_t firstEta = eta[i];
+ Float_t firstPhi = phi[i];
+ Float_t firstPt = pt[i];
+ Float_t gPhisMinusPsiBin = phiMinusPsi[i];
+ for (Int_t j = 0; j < nMCLabelCounter ; j++) {
+ if(i == j) continue;
+
+ Float_t secondEta = eta[j];
+ Float_t secondPhi = phi[j];
+ Float_t secondPt = pt[j];
+
+ trackVariablesPair[0] = gPhisMinusPsiBin;
+ trackVariablesPair[1] = firstEta - secondEta; // delta eta
+ trackVariablesPair[2] = firstPhi - secondPhi; // delta phi
+ if (trackVariablesPair[2] > 180.) // delta phi between -180 and 180
+ trackVariablesPair[2] -= 360.;
+ if (trackVariablesPair[2] < - 180.)
+ trackVariablesPair[2] += 360.;
+ trackVariablesPair[3] = firstPt; // pt trigger
+ trackVariablesPair[4] = secondPt; // pt
+
+ //++ pairs
+ if(charge[i] > 0 && charge[j] > 0 ) {
+ if(level[i] > 0 && level[j] > 0)
+ fHistGeneratedPlusPlus->Fill(trackVariablesPair,0,1.);
+
+ if(level[i] > 1 && level[j] > 1)
+ fHistSurvivedPlusPlus->Fill(trackVariablesPair,0,1.);
+ }
+
+ //-- pairs
+ else if(charge[i] < 0 && charge[j] < 0 ) {
+ if(level[i] > 0 && level[j] > 0)
+ fHistGeneratedMinusMinus->Fill(trackVariablesPair,0,1.);
+
+ if(level[i] > 1 && level[j] > 1)
+ fHistSurvivedMinusMinus->Fill(trackVariablesPair,0,1.);
+ }
+
+ //+- pairs
+ else if(charge[i] > 0 && charge[j] < 0 ) {
+ if(level[i] > 0 && level[j] > 0)
+ fHistGeneratedPlusMinus->Fill(trackVariablesPair,0,1.);
+
+ if(level[i] > 1 && level[j] > 1)
+ fHistSurvivedPlusMinus->Fill(trackVariablesPair,0,1.);
+ }
+
+ //-+ pairs
+ else if(charge[i] < 0 && charge[j] > 0 ) {
+ if(level[i] > 0 && level[j] > 0)
+ fHistGeneratedMinusPlus->Fill(trackVariablesPair,0,1.);
+
+ if(level[i] > 1 && level[j] > 1)
+ fHistSurvivedMinusPlus->Fill(trackVariablesPair,0,1.);
+ }
+ }//second particle loop
+ }//first particle loop
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEfficiencyBFPsi::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 AliAnalysisTaskEfficiencyBFPsi::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;
+}
--- /dev/null
+#ifndef ALIANALYSISTASKEFFICIENCYBFPSI_cxx
+#define ALIANALYSISTASKEFFICIENCYBFPSI_cxx
+
+// ---------------------------------------------------------------------
+//
+// Task for calculating the efficiency of the Balance Function
+// for single particles and pairs: multi-D analysis
+//
+// Authors: Panos Christakoglou, Michael Weber
+//
+// ---------------------------------------------------------------------
+
+class TH1F;
+class TH3D;
+class TH2F;
+class TString;
+class AliESDEvent;
+class AliESDtrackCuts;
+#include "AliTHn.h"
+
+#include "AliAnalysisTaskSE.h"
+
+const Int_t kVariablesSingle = 2; // track variables in histogram (phi-Psi2, pTtrig)
+const Int_t kVariablesPair = 5; // track variables in histogram (phi-Psi2, dEta, dPhi, pTtrig, ptAssociated)
+
+class AliAnalysisTaskEfficiencyBFPsi : public AliAnalysisTaskSE {
+ public:
+ AliAnalysisTaskEfficiencyBFPsi() : AliAnalysisTaskSE(),
+ fESD(0), fQAList(0), fOutputList(0),
+ fHistEventStats(0), fHistCentrality(0), fHistNMult(0),
+ fHistGeneratedPlus(0), fHistSurvivedPlus(0),
+ fHistGeneratedMinus(0),fHistSurvivedMinus(0),
+ fHistGeneratedPlusPlus(0), fHistSurvivedPlusPlus(0),
+ fHistGeneratedMinusMinus(0), fHistSurvivedMinusMinus(0),
+ fHistGeneratedPlusMinus(0), fHistSurvivedPlusMinus(0),
+ fHistGeneratedMinusPlus(0), fHistSurvivedMinusPlus(0),
+ fESDtrackCuts(0), fAnalysisMode(0),
+ fCentralityEstimator("V0M"),
+ fCentralityPercentileMin(0.0), fCentralityPercentileMax(10.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), fEtaRangeMax(0.8),
+ fPtRangeMin(0.1), fPtRangeMax(5.0), fPhiRangeMin(0.0),fPhiRangeMax(6.28){}
+ AliAnalysisTaskEfficiencyBFPsi(const char *name);
+ virtual ~AliAnalysisTaskEfficiencyBFPsi() {}
+
+ 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;}
+
+ void SetEtaRangeMax(Double_t maxRangeEta){
+ fEtaRangeMax = maxRangeEta;}//
+ void SetPtRangeMin(Double_t minRangePt){
+ fPtRangeMin = minRangePt;} //
+ void SetPtRangeMax(Double_t maxRangePt){
+ fPtRangeMax = maxRangePt;} //
+ void SetPhiRangeMin(Double_t minRangePhi){
+ fPhiRangeMin = minRangePhi;} //
+ void SetPhiRangeMax(Double_t maxRangePhi){
+ fPhiRangeMax = maxRangePhi;} //
+
+
+ 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)
+ AliTHn *fHistGeneratedPlus;//!correction map for positives (generated)
+ AliTHn *fHistSurvivedPlus;//!correction map positives (survived)
+ AliTHn *fHistGeneratedMinus;//!correction map for negatives (generated)
+ AliTHn *fHistSurvivedMinus;//!correction map negatives (survived)
+
+ // output histograms (pairs)
+ AliTHn *fHistGeneratedPlusPlus;//!correction map for ++ (generated)
+ AliTHn *fHistSurvivedPlusPlus;//!correction map ++ (survived)
+ AliTHn *fHistGeneratedMinusMinus;//!correction map for -- (generated)
+ AliTHn *fHistSurvivedMinusMinus;//!correction map -- (survived)
+ AliTHn *fHistGeneratedPlusMinus;//!correction map for +- (generated)
+ AliTHn *fHistSurvivedPlusMinus;//!correction map +- (survived)
+ AliTHn *fHistGeneratedMinusPlus;//!correction map for -+ (generated)
+ AliTHn *fHistSurvivedMinusPlus;//!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;
+ Double_t fEtaRangeMax; // acceptance cuts
+ Double_t fPtRangeMin; //acceptance cuts
+ Double_t fPtRangeMax; //acceptance cuts
+ Double_t fPhiRangeMin; //acceptance cuts
+ Double_t fPhiRangeMax; // acceptance cuts
+
+ AliAnalysisTaskEfficiencyBFPsi(const AliAnalysisTaskEfficiencyBFPsi&); // not implemented
+ AliAnalysisTaskEfficiencyBFPsi& operator=(const AliAnalysisTaskEfficiencyBFPsi&); // not implemented
+
+ ClassDef(AliAnalysisTaskEfficiencyBFPsi, 1); // example of analysis
+};
+
+#endif
--- /dev/null
+enum anaModes {mLocal,mLocalPAR,mPROOF,mGrid,mGridPAR};
+//mLocal: Analyze locally files in your computer using aliroot
+//mLocalPAR: Analyze locally files in your computer using root + PAR files
+//mPROOF: Analyze CAF files with PROOF
+//mGrid: Analyze files on Grid via AliEn plug-in and using precompiled FLOW libraries
+//mGridPAR: Analyze files on Grid via AliEn plug-in and using par files for FLOW package
+
+//Analysis modes
+const TString analysisMode = "TPC"; //"TPC", "Global"
+
+//Centrality stuff
+Int_t binfirst = 0; //where do we start numbering bins
+Int_t binlast = 8; //where do we stop numbering bins
+const Int_t numberOfCentralityBins = 9;
+Double_t centralityArray[numberOfCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; // in centrality percentile
+Bool_t kUsePID = kFALSE;
+const TString centralityEstimator = "V0M";
+Double_t vertexZ = 10.;
+
+//output file
+TString commonOutputFileName = "outputEfficiency";
+
+//void runEfficiencyTaskCentralityTrain(Int_t mode = mPROOF,
+//Int_t nRuns = 600000,
+//Bool_t DATA = kFALSE,
+//const Char_t* dataDir="/alice/data/LHC10h_000137161_p1_4plus#esdTree", Int_t offset=0) {
+void runBalancePsiEfficiencyTaskCentralityTrain(Int_t mode = mLocal, Bool_t DATA = kFALSE) {
+//void runEfficiencyTaskCentralityTrain(const char* runListFileName = "listOfRuns.txt",Int_t mode = mGrid, Bool_t DATA = kFALSE) {
+ // Time:
+ TStopwatch timer;
+ timer.Start();
+
+ // Load needed libraries:
+ LoadLibraries(mode);
+
+ // Create and configure the AliEn plug-in:
+ if(mode == mGrid || mode == mGridPAR) {
+ gROOT->LoadMacro("CreateAlienHandler.C");
+ AliAnalysisGrid *alienHandler = CreateAlienHandler(runListFileName);
+ if (!alienHandler) return;
+ gROOT->LoadMacro("AliAnalysisTaskEfficiencyBF.cxx++"); //gROOT->LoadMacro("AliAnalysisEfficiencyTaskCentralityTrain.cxx++");
+ }
+ // Chains:
+ if(mode==mLocal || mode == mLocalPAR) {
+ //gROOT->LoadMacro("AliAnalysisEfficiencyTaskCentralityTrain.cxx++");
+ //gROOT->LoadMacro("AliAnalysisTaskEfficiencyBF.cxx++");
+ TChain* chain = new TChain("esdTree");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set1/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set2/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set3/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set4/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set5/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set6/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set7/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set8/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set9/AliESDs.root");
+ chain->Add("/glusterfs/alice1/alice2/pchrist/HeavyIons/MC/Set10/AliESDs.root");
+ }
+ //Proof
+ if(mode == mPROOF) {
+ gROOT->ProcessLine(Form(".include %s/include", gSystem->ExpandPathName("$ALICE_ROOT")));
+ gProof->Load("AliAnalysisTaskEfficiencyBF.cxx++");
+ }
+
+ // Create analysis manager:
+ AliAnalysisManager *mgr = new AliAnalysisManager("bfPsiCorrectionManager");
+ // Connect plug-in to the analysis manager:
+ if(mode == mGrid || mode == mGridPAR) {
+ mgr->SetGridHandler(alienHandler);
+ }
+ // Event handlers:
+ AliVEventHandler* esdH = new AliESDInputHandler;
+ mgr->SetInputEventHandler(esdH);
+ AliMCEventHandler *mc = new AliMCEventHandler();
+ //mc->SetReadTR(kFALSE);
+ mgr->SetMCtruthEventHandler(mc);
+
+ // Task to check the offline trigger:
+ //if(mode == mLocal || mode == mGrid || mode == mGridPAR)
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ AliPhysicsSelectionTask* physicsSelTask = AddTaskPhysicsSelection(!DATA);
+ physicsSelTask->GetPhysicsSelection()->SetAnalyzeMC();
+ // Enable debug printouts:
+ mgr->SetDebugLevel(2);
+
+ //Add the centrality determination task
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+ AliCentralitySelectionTask *centralityTask = AddTaskCentrality();
+ centralityTask->SetMCInput();
+ //centralityTask->SetPass(2);
+ //AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
+ //taskCentrality->SelectCollisionCandidates(AliVEvent::kMB);
+
+ if(kUsePID) {
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+ AddTaskPIDResponse(1);
+ }
+
+ // Load the analysis task:
+ gROOT->LoadMacro("$ALICE_ROOT/PWGCF/EBYE/macros/AddTaskBalancePsiEfficiency.C");
+
+ //for (Int_t i=binfirst; i<binlast+1; i++) {
+ for (Int_t i = 0; i < 9; i++) {
+ Double_t lowCentralityBinEdge = centralityArray[i];
+ Double_t highCentralityBinEdge = centralityArray[i+1];
+ Printf("\nWagon for centrality bin %i: %.0f-%.0f",i,lowCentralityBinEdge,highCentralityBinEdge);
+ AddTaskBalancePsiEfficiency(centralityEstimator.Data(), lowCentralityBinEdge, highCentralityBinEdge, vertexZ, commonOutputFileName.Data());
+ } // end of for (Int_t i=0; i<numberOfCentralityBins; i++)
+
+ // Run the analysis:
+ if(!mgr->InitAnalysis()){return;}
+ mgr->PrintStatus();
+ if(mode == mLocal || mode == mLocalPAR)
+ mgr->StartAnalysis("local",chain);
+ else if(mode == mPROOF)
+ mgr->StartAnalysis("proof",dataDir,nRuns,offset);
+ else if(mode == mGrid || mode == mGridPAR)
+ mgr->StartAnalysis("grid");
+
+ // Print real and CPU time used for analysis:
+ timer.Stop();
+ timer.Print();
+
+} // end of void runTaskFluctuations(...)
+
+//=============================================================//
+void LoadLibraries(const anaModes mode) {
+ //--------------------------------------
+ // Load the needed libraries most of them already loaded by aliroot
+ //--------------------------------------
+ gSystem->Load("libTree");
+ gSystem->Load("libGeom");
+ gSystem->Load("libVMC");
+ gSystem->Load("libXMLIO");
+ gSystem->Load("libPhysics");
+
+ //----------------------------------------------------------
+ // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
+ //----------------------------------------------------------
+ if (mode==mLocal || mode==mGrid || mode == mGridPAR) {
+ //--------------------------------------------------------
+ // If you want to use already compiled libraries
+ // in the aliroot distribution
+ //--------------------------------------------------------
+ gSystem->Load("libSTEERBase");
+ gSystem->Load("libESD");
+ gSystem->Load("libAOD");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libEventMixing.so");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libPWGTools.so");
+ gSystem->Load("libPWGCFebye.so");
+
+ // Use AliRoot includes to compile our task
+ gROOT->ProcessLine(".include $ALICE_ROOT/include");
+ }
+
+ else if (mode == mLocalPAR) {
+ //--------------------------------------------------------
+ //If you want to use root and par files from aliroot
+ //--------------------------------------------------------
+ SetupPar("STEERBase");
+ SetupPar("ESD");
+ SetupPar("AOD");
+ SetupPar("ANALYSIS");
+ SetupPar("ANALYSISalice");
+}
+
+ //---------------------------------------------------------
+ // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
+ //---------------------------------------------------------
+ else if (mode==mPROOF) {
+ // Connect to proof
+ printf("*** Connect to PROOF ***\n");
+ gEnv->SetValue("XSec.GSI.DelegProxy","2");
+ // Put appropriate username here
+ TProof::Open("alice-caf.cern.ch");
+ //TProof::Open("skaf.saske.sk");
+ //TProof::Open("prf000-iep-grid.saske.sk");
+
+ gProof->EnablePackage("VO_ALICE@AliRoot::v4-21-12-AN");
+ }
+
+} // end of void LoadLibraries(const anaModes mode)
+
+//===============================================================================================
+
+void SetupPar(char* pararchivename) {
+ //Load par files, create analysis libraries
+ //For testing, if par file already decompressed and modified
+ //classes then do not decompress.
+
+ TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
+ TString parpar(Form("%s.par", pararchivename)) ;
+ if ( gSystem->AccessPathName(parpar.Data()) ) {
+ gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
+ TString processline(Form(".! make %s", parpar.Data())) ;
+ gROOT->ProcessLine(processline.Data()) ;
+ gSystem->ChangeDirectory(cdir) ;
+ processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
+ gROOT->ProcessLine(processline.Data()) ;
+ }
+ if ( gSystem->AccessPathName(pararchivename) ) {
+ TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
+ gROOT->ProcessLine(processline.Data());
+ }
+
+ TString ocwd = gSystem->WorkingDirectory();
+ gSystem->ChangeDirectory(pararchivename);
+
+ // check for BUILD.sh and execute
+ if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+ printf("*******************************\n");
+ printf("*** Building PAR archive ***\n");
+ cout<<pararchivename<<endl;
+ printf("*******************************\n");
+ if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+ Error("runProcess","Cannot Build the PAR Archive! - Abort!");
+ return -1;
+ }
+ }
+ // check for SETUP.C and execute
+ if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+ printf("*******************************\n");
+ printf("*** Setup PAR archive ***\n");
+ cout<<pararchivename<<endl;
+ printf("*******************************\n");
+ gROOT->Macro("PROOF-INF/SETUP.C");
+ }
+
+ gSystem->ChangeDirectory(ocwd.Data());
+ printf("Current dir: %s\n", ocwd.Data());
+
+} // end of void SetupPar(char* pararchivename)
+
+//===============================================================================================
+
+// Helper macros for creating chains
+// from: CreateESDChain.C,v 1.10 jgrosseo Exp
+
+TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+ // creates chain of files in a given directory or file containing a list.
+ // In case of directory the structure is expected as:
+ // <aDataDir>/<dir0>/AliESDs.root
+ // <aDataDir>/<dir1>/AliESDs.root
+ // ...
+
+ if (!aDataDir)
+ return 0;
+
+ Long_t id, size, flags, modtime;
+ if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+ {
+ printf("%s not found.\n", aDataDir);
+ return 0;
+ }
+
+ TChain* chain = new TChain("esdTree");
+ TChain* chaingAlice = 0;
+
+ if (flags & 2)
+ {
+ TString execDir(gSystem->pwd());
+ TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+ TList* dirList = baseDir->GetListOfFiles();
+ Int_t nDirs = dirList->GetEntries();
+ gSystem->cd(execDir);
+
+ Int_t count = 0;
+
+ for (Int_t iDir=0; iDir<nDirs; ++iDir)
+ {
+ TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+ if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+ continue;
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ TString presentDirName(aDataDir);
+ presentDirName += "/";
+ presentDirName += presentDir->GetName();
+ chain->Add(presentDirName + "/AliESDs.root/esdTree");
+ // cerr<<presentDirName<<endl;
+ }
+
+ }
+ else
+ {
+ // Open the input stream
+ ifstream in;
+ in.open(aDataDir);
+
+ Int_t count = 0;
+
+ // Read the input list of files and add them to the chain
+ TString esdfile;
+ while(in.good()) {
+ in >> esdfile;
+ if (!esdfile.Contains("root")) continue; // protection
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ // add esd file
+ chain->Add(esdfile);
+ }
+
+ in.close();
+ }
+
+ return chain;
+
+} // end of TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+
+//===============================================================================================
+
+TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+ // creates chain of files in a given directory or file containing a list.
+ // In case of directory the structure is expected as:
+ // <aDataDir>/<dir0>/AliAOD.root
+ // <aDataDir>/<dir1>/AliAOD.root
+ // ...
+
+ if (!aDataDir)
+ return 0;
+
+ Long_t id, size, flags, modtime;
+ if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+ {
+ printf("%s not found.\n", aDataDir);
+ return 0;
+ }
+
+ TChain* chain = new TChain("aodTree");
+ TChain* chaingAlice = 0;
+
+ if (flags & 2)
+ {
+ TString execDir(gSystem->pwd());
+ TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+ TList* dirList = baseDir->GetListOfFiles();
+ Int_t nDirs = dirList->GetEntries();
+ gSystem->cd(execDir);
+
+ Int_t count = 0;
+
+ for (Int_t iDir=0; iDir<nDirs; ++iDir)
+ {
+ TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+ if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+ continue;
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ TString presentDirName(aDataDir);
+ presentDirName += "/";
+ presentDirName += presentDir->GetName();
+ chain->Add(presentDirName + "/AliAOD.root/aodTree");
+ // cerr<<presentDirName<<endl;
+ }
+
+ }
+ else
+ {
+ // Open the input stream
+ ifstream in;
+ in.open(aDataDir);
+
+ Int_t count = 0;
+
+ // Read the input list of files and add them to the chain
+ TString aodfile;
+ while(in.good()) {
+ in >> aodfile;
+ if (!aodfile.Contains("root")) continue; // protection
+
+ if (offset > 0)
+ {
+ --offset;
+ continue;
+ }
+
+ if (count++ == aRuns)
+ break;
+
+ // add aod file
+ chain->Add(aodfile);
+ }
+
+ in.close();
+ }
+
+ return chain;
+
+} // end of TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff