[u/mrichter/AliRoot.git] / PWG2 / EBYE / MeanPtFluctuations / AliAnalysisTaskPtFlucPbPb.cxx

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TProfile.h"
#include "TList.h"
#include "iostream"

#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"

#include "AliESD.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrackCuts.h"
#include "AliCentrality.h"

#include "AliLog.h"

#include "AliAnalysisTaskPtFlucPbPb.h"


using namespace std;

// Analysis of Pt FLuctuations (PbPb)
// Author: Stefan Heckel
// Version of PbPb task: 5.0, 18.04.2011


ClassImp(AliAnalysisTaskPtFlucPbPb)

//________________________________________________________________________
AliAnalysisTaskPtFlucPbPb::AliAnalysisTaskPtFlucPbPb(const char *name)
  :AliAnalysisTaskSE(name),
  fESD(0),
  fOutputList(0),
  fPtSpec(0),
  fMult(0),
  fMultSum(0),
  fMultSumPt(0),
  fMultNrPairs(0),
  fCent(0),
  fCentSum(0),
  fCentSumPt(0),
  fCentNrPairs(0),
  fEta(0),
  fEtaPhiPlus(0),
  fEtaPhiMinus(0),
  fVtxZ(0),
  fVtxZCut(0),
  fVtxZCont(0),
  fVtxZTrackCuts(0),
  fEventMeanPt(0),
  fEventMeanPtSq(0),
  fMultEventMeanPt(0),
  fMultEventMeanPtSq(0),
  fCentEventMeanPt(0),
  fCentEventMeanPtSq(0),
  fTwoPartCorrEv(0),
  fTwoPartCorrEvSq(0),
  fTwoPartCorrEvCent(0),
  fTwoPartCorrEvCentSq(0),
  fESDTrackCuts(0),
  fMaxVertexZ(0),
  fNContributors(0),
  fUseCentrality(0),
  fMC(0)
{
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
AliAnalysisTaskPtFlucPbPb::~AliAnalysisTaskPtFlucPbPb()
{
  if(fOutputList) delete fOutputList;  fOutputList =0;
}

//________________________________________________________________________
void AliAnalysisTaskPtFlucPbPb::UserCreateOutputObjects()
{
  // Create histograms
  // Called once

  OpenFile(1, "RECREATE");
  fOutputList = new TList();
  fOutputList->SetOwner();


  fPtSpec = new TH1F("fPtSpec","Pt spectrum",50,0,2.5);

  fMult = new TH1F("fMult","Multiplicity distribution",30,0,3000);

  fMultSum = new TH1F("fMultSum","Sum of nrtracks of all events in multiplicity bins",30,0,3000);

  fMultSumPt = new TH1F("fMultSumPt","Sum of pTs of all events in multiplicity bins",30,0,3000);

  fMultNrPairs = new TH1F("fMultNrPairs","Sum of number of pairs in multiplicity bins",30,0,3000);

  fCent = new TH1F("fCent","Centrality distribution",21,0,105);

  fCentSum = new TH1F("fCentSum","Sum of nrtracks of all events in centrality bins",21,0,105);

  fCentSumPt = new TH1F("fCentSumPt","Sum of pTs of all events in centrality bins",21,0,105);

  fCentNrPairs = new TH1F("fCentNrPairs","Sum of number of pairs in centrality bins",21,0,105);


  fEta = new TH1F("fEta","Eta distribution",80,-2,2);

  fEtaPhiPlus = new TH1F("fEtaPhiPlus","Phi distribution for positive eta",62,0,6.2);

  fEtaPhiMinus = new TH1F("fEtaPhiMinus","Phi distribution for negative eta",62,0,6.2);

  fVtxZ = new TH1F("fVtxZ","Vertex Z distribution before cuts",100,-20,20);

  fVtxZCut = new TH1F("fVtxZCut","Vertex Z distribution after vtxZ cut",110,-11,11);

  fVtxZCont = new TH1F("fVtxZCont","Vertex Z distribution after nCont cut",110,-11,11);

  fVtxZTrackCuts = new TH1F("fVtxZTrackCuts","Vertex Z distribution after track cuts",110,-11,11);


  fEventMeanPt = new TH1F("fEventMeanPt","Mean-Pt event by event",50,0,2.5);

  fEventMeanPtSq = new TH1F("fEventMeanPtSq","Mean-Pt event by event squared",100,0,5);

  fMultEventMeanPt = new TH1F("fMultEventMeanPt","Mean-Pt event by event vs. multiplicity",30,0,3000);

  fMultEventMeanPtSq = new TH1F("fMultEventMeanPtSq","Mean-Pt event by event squared vs. multiplicity",30,0,3000);

  fCentEventMeanPt = new TH1F("fCentEventMeanPt","Mean-Pt event by event vs. centrality",21,0,105);

  fCentEventMeanPtSq = new TH1F("fCentEventMeanPtSq","Mean-Pt event by event squared vs. centrality",21,0,105);


  fTwoPartCorrEv = new TH1F("fTwoPartCorrEv","Two-particle correlator vs. multiplicity",30,0,3000);

  fTwoPartCorrEvSq = new TH1F("fTwoPartCorrEvSq","Two-particle correlator squared vs. multiplicity",30,0,3000);

  fTwoPartCorrEvCent = new TH1F("fTwoPartCorrEvCent","Two-particle correlator vs. centrality",21,0,105);

  fTwoPartCorrEvCentSq = new TH1F("fTwoPartCorrEvCentSq","Two-particle correlator squared vs. centrality",21,0,105);


  // Add histograms to the output list
  fOutputList->Add(fPtSpec);
  fOutputList->Add(fMult);
  fOutputList->Add(fMultSum);
  fOutputList->Add(fMultSumPt);
  fOutputList->Add(fMultNrPairs);
  fOutputList->Add(fCent);
  fOutputList->Add(fCentSum);
  fOutputList->Add(fCentSumPt);
  fOutputList->Add(fCentNrPairs);
  fOutputList->Add(fEta);
  fOutputList->Add(fEtaPhiPlus);
  fOutputList->Add(fEtaPhiMinus);
  fOutputList->Add(fVtxZ);
  fOutputList->Add(fVtxZCut);
  fOutputList->Add(fVtxZCont);
  fOutputList->Add(fVtxZTrackCuts);
  fOutputList->Add(fEventMeanPt);
  fOutputList->Add(fEventMeanPtSq);
  fOutputList->Add(fMultEventMeanPt);
  fOutputList->Add(fMultEventMeanPtSq);
  fOutputList->Add(fCentEventMeanPt);
  fOutputList->Add(fCentEventMeanPtSq);
  fOutputList->Add(fTwoPartCorrEv);
  fOutputList->Add(fTwoPartCorrEvSq);
  fOutputList->Add(fTwoPartCorrEvCent);
  fOutputList->Add(fTwoPartCorrEvCentSq);
}

//________________________________________________________________________
void AliAnalysisTaskPtFlucPbPb::UserExec(Option_t *)
{
  // Main loop
  // Called for each event


  // --- ESD event handler ---

  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

  if (!esdH) {
    Printf("ERROR: Could not get ESDInputHandler");
  }
  else {
    fESD = esdH->GetEvent();
  }

  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }


  // --- End event handler ---


  // --- Vertex cuts ---

//   // TPC+ITS vertex
//   const AliESDVertex* vtxESD = fESD->GetPrimaryVertexTracks();
//   Double_t vtxZ = vtxESD->GetZv();
//   Double_t vtxNCont = vtxESD->GetNContributors();

  // TPConly vertex
  const AliESDVertex* vtxESDTPC = fESD->GetPrimaryVertexTPC();
  Double_t vtxZ = vtxESDTPC->GetZv();
  Double_t vtxNCont = vtxESDTPC->GetNContributors();

  fVtxZ->Fill(vtxZ); // VtxZ before cuts

  // Event cut on the z-Position of the vertex
  if(vtxZ > fMaxVertexZ || vtxZ < (-1.*fMaxVertexZ)) {
//     Printf("VertexZ out of range, Zv = %f",vtxZ);
    return;
  }

  fVtxZCut->Fill(vtxZ); // VtxZ after cut on vtxZ

  // Event cut on the number of contributors
  if(vtxNCont < fNContributors) {
//     Printf("Vertex has no contributors");
    return;
  }
//   else {
//     Printf("Vertex nContributors = %i",vtxNCont);
//   }

  fVtxZCont->Fill(vtxZ); // VtxZ after cut on nContributors

  // --- End vertex cuts ---


  Int_t nrTracks = 0;		// count for number of tracks which pass the track cuts
  Int_t nrESDTracks = 0;	// number of tracks in the ESD file
  Double_t sumPt = 0.;		// sum of track Pts
  Double_t twoPartCorrEv = 0.;	// Two-particle correlator of one event for multiplicity bin analysis
  Double_t twoPartCorrEvCent = 0.; // Two-particle correlator of one event for centrality bin analysis

  Double_t tracks[3000] = {0.};	// array of track Pts, needed for the two-particle correlator

  Double_t trackPt = 0., trackEta = 0., trackPhi = 0.;
  Double_t eventMeanPt = 0., eventMeanPtSq = 0., evMptMult = 0.;
  Double_t twoPartCorrPair = 0., twoPartCorrEvSq = 0.;
  Double_t twoPartCorrPairCent = 0., twoPartCorrEvCentSq = 0.;
  Double_t nrPairs = 0.;

  Double_t *nbins = 0x0; // Mean pT values for multiplicity bin analysis
  Double_t *centbins = 0x0; // Mean pT values for centrality bin analysis

  Double_t centralityVZERO = 0.;
  Int_t centralityVZEROBin = 0;

  Int_t centBin = 0;
  Double_t evMptCent = 0.;


//   Printf("MC: %d",fMC);


// --- Mean pT values ---

// !!!!! Have to be calculated in a first run - for each sample, which should be analysed, separately! !!!!!


// -- Mean pT values for multiplicity bins (first bin is for nrTracks = 0 and has to be 0) --

  // MC PbPb 2.76 ATeV 11a10a (Hijing) (Pt-Range: 0.15 - 2)
  Double_t nbinsMC276[32] = {0.000, 0.525, 0.535, 0.539, 0.541, 0.543, 0.544, 0.545, 0.545, 0.546, 0.547, 0.548, 0.548, 0.549, 0.549, 0.550, 0.551, 0.551, 0.552, 0.552, 0.553, 0.554, 0.554, 0.555, 0.556, 0.556, 0.557, 0.558, 0.549, 0.000, 0.000, 0.000};

  // Data PbPb 2.76 ATeV 10h.pass1 (Pt-Range: 0.15 - 2)
  Double_t nbinsData276[32] = {0.000, 0.580, 0.605, 0.618, 0.626, 0.632, 0.636, 0.640, 0.642, 0.645, 0.647, 0.648, 0.649, 0.650, 0.651, 0.652, 0.653, 0.653, 0.654, 0.654, 0.654, 0.654, 0.654, 0.654, 0.655, 0.656, 0.658, 0.660, 0.661, 0.659, 0.660, 0.656};

// -- End mean pT values for multiplicity bins --


// -- Mean pT values for centrality bins --

  // MC PbPb 2.76 ATeV 11a10a (Hijing) (Pt-Range: 0.15 - 2)
  Double_t centbinsMC276[21] = {0.555, 0.553, 0.551, 0.549, 0.548, 0.547, 0.545, 0.544, 0.543, 0.541, 0.539, 0.537, 0.536, 0.533, 0.530, 0.528, 0.524, 0.518, 0.505, 0.000, 0.534};

  // Data PbPb 2.76 ATeV 10h.pass1 (Pt-Range: 0.15 - 2)
  Double_t centbinsData276[21] = {0.654, 0.654, 0.652, 0.650, 0.647, 0.644, 0.640, 0.635, 0.629, 0.623, 0.616, 0.609, 0.601, 0.594, 0.586, 0.579, 0.568, 0.551, 0.535, 0.000, 0.621};

// -- End mean pT values for centrality bins --


// -- Selection of MC/Data; whole sample values --

if (fMC) { // - MC -

//   Printf(" -- MC, 2.76 ATeV -- ");

  nbins = nbinsMC276;
  centbins = centbinsMC276;

} // - End MC -
else { // - Data -

//   Printf(" -- Data, 2.76 ATeV -- ");

  nbins = nbinsData276;
  centbins = centbinsData276;

} // - End data -

// -- End selection of MC/Data; whole sample values --

// --- End mean pT values ---


  nrESDTracks = fESD->GetNumberOfTracks();
//   if( nrESDTracks ) { printf("Found event with %i tracks.\n",nrESDTracks); }


  if(fUseCentrality != 0) {

    AliCentrality *esdCentrality = fESD->GetCentrality();

    //V0
    if (fUseCentrality == 1){

      centralityVZERO = esdCentrality->GetCentralityPercentile("V0M");

      if      ( centralityVZERO >   0. && centralityVZERO <   5.) centralityVZEROBin =  0;
      else if ( centralityVZERO >=  5. && centralityVZERO <  10.) centralityVZEROBin =  5;
      else if ( centralityVZERO >= 10. && centralityVZERO <  15.) centralityVZEROBin = 10;
      else if ( centralityVZERO >= 15. && centralityVZERO <  20.) centralityVZEROBin = 15;
      else if ( centralityVZERO >= 20. && centralityVZERO <  25.) centralityVZEROBin = 20;
      else if ( centralityVZERO >= 25. && centralityVZERO <  30.) centralityVZEROBin = 25;
      else if ( centralityVZERO >= 30. && centralityVZERO <  35.) centralityVZEROBin = 30;
      else if ( centralityVZERO >= 35. && centralityVZERO <  40.) centralityVZEROBin = 35;
      else if ( centralityVZERO >= 40. && centralityVZERO <  45.) centralityVZEROBin = 40;
      else if ( centralityVZERO >= 45. && centralityVZERO <  50.) centralityVZEROBin = 45;
      else if ( centralityVZERO >= 50. && centralityVZERO <  55.) centralityVZEROBin = 50;
      else if ( centralityVZERO >= 55. && centralityVZERO <  60.) centralityVZEROBin = 55;
      else if ( centralityVZERO >= 60. && centralityVZERO <  65.) centralityVZEROBin = 60;
      else if ( centralityVZERO >= 65. && centralityVZERO <  70.) centralityVZEROBin = 65;
      else if ( centralityVZERO >= 70. && centralityVZERO <  75.) centralityVZEROBin = 70;
      else if ( centralityVZERO >= 75. && centralityVZERO <  80.) centralityVZEROBin = 75;
      else if ( centralityVZERO >= 80. && centralityVZERO <  85.) centralityVZEROBin = 80;
      else if ( centralityVZERO >= 85. && centralityVZERO <  90.) centralityVZEROBin = 85;
      else if ( centralityVZERO >= 90. && centralityVZERO <  95.) centralityVZEROBin = 90;
      else if ( centralityVZERO >= 95. && centralityVZERO <  99.) centralityVZEROBin = 95;
      else if ( centralityVZERO >= 99. ) centralityVZEROBin = 100;
      else if ( centralityVZERO <= 0.  ) centralityVZEROBin = 100;

      centBin = centralityVZEROBin / 5;
      evMptCent = centbins[centBin];
    }

//     Printf("\nRaw mult: %i   CentVZERO: %f   CentVZERObin: %i   centBin: %i   evMptCent: %.3f ",nrESDTracks,centralityVZERO,centralityVZEROBin,centBin,evMptCent);

  }


  // --- Loop over all tracks of one event ---

  for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
    AliESDtrack* track = fESD->GetTrack(iTracks);
    if (!track) {
      Printf("ERROR: Could not receive track %d\n", iTracks);
      continue;
    }
    if(!fESDTrackCuts) {
      Printf("ERROR: No esd track cut");
      continue;
    }
    else {
      if(!fESDTrackCuts->AcceptTrack(track))
	continue;
    }
	trackPt = track->Pt();
	fPtSpec->Fill(trackPt);
	tracks[nrTracks] = trackPt;

	trackEta = track->Eta();
	fEta->Fill(trackEta);

	trackPhi = track->Phi();

	if (trackEta > 0.) {
	  fEtaPhiPlus->Fill(trackPhi);
	}
	else if (trackEta < 0.) {
	  fEtaPhiMinus->Fill(trackPhi);
	}

	sumPt = sumPt + trackPt;
	nrTracks++;

// 	printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);

  } // --- End track loop ---


  // --- Calculation of various values and filling of histograms (for remaining events with N_acc > 0) ---

  if(nrTracks != 0) {

    // VtxZ after all track cuts
    fVtxZTrackCuts->Fill(vtxZ);

    // Multiplicity distributions
    fMult->Fill(nrTracks);
    fMultSum->Fill(nrTracks,nrTracks);

    // Number of pairs in event
    nrPairs = 0.5 * nrTracks * (nrTracks-1);
    fMultNrPairs->Fill(nrTracks,nrPairs);

    // Calculation of mean Pt and mean Pt Squared
    eventMeanPt = sumPt / nrTracks;
    eventMeanPtSq = eventMeanPt * eventMeanPt;

    // Mean-Pt and Mean-Pt squared
    fEventMeanPt->Fill(eventMeanPt);
    fEventMeanPtSq->Fill(eventMeanPtSq);

    // Mean-Pt and Mean-Pt squared depending on multiplicity
    fMultEventMeanPt->Fill(nrTracks,eventMeanPt);
    fMultEventMeanPtSq->Fill(nrTracks,eventMeanPtSq);
    fMultSumPt->Fill(nrTracks,sumPt);

//     printf("nrTracks: %i   sumPt: %.8f   meanPt: %.8f   meanPtSq: %.8f\n",nrTracks,sumPt,eventMeanPt,eventMeanPtSq);


    // Centrality V0
    if (fUseCentrality == 1) {

      // Centrality distributions
      fCent->Fill(centralityVZEROBin);
      fCentSum->Fill(centralityVZEROBin,nrTracks);

      // Number of pairs in event
      fCentNrPairs->Fill(centralityVZEROBin,nrPairs);

      // Mean-Pt and Mean-Pt squared depending on centrality
      fCentEventMeanPt->Fill(centralityVZEROBin,eventMeanPt);
      fCentEventMeanPtSq->Fill(centralityVZEROBin,eventMeanPtSq);
      fCentSumPt->Fill(centralityVZEROBin,sumPt);

//       printf("CentV0bin: %i   NrTracks: %i   NrPairs: %.0f   SumPt: %.1f\n",centralityVZEROBin,nrTracks,nrPairs,sumPt);
    }


    // --- Two-particle correlator ---

    // -- Analysis in multiplicity bins --

    for (int k=1; k<32; k++) {
      if (nrTracks > 100 * (k-1) && nrTracks <= 100 * k) {
	evMptMult = nbins[k];
	break;
      }
    }
//     printf("nrTracks = %3i   evMptMult = %.3f \n",nrTracks,evMptMult);

    for (int i=0; i<nrTracks; i++) {
      for (int j=i+1; j<nrTracks; j++) {
	twoPartCorrPair = (tracks[i] - evMptMult) * (tracks[j] - evMptMult);
	twoPartCorrEv = twoPartCorrEv + twoPartCorrPair;
      }
    }

    twoPartCorrEvSq = twoPartCorrEv * twoPartCorrEv;
    fTwoPartCorrEv->Fill(nrTracks,twoPartCorrEv);
    fTwoPartCorrEvSq->Fill(nrTracks,twoPartCorrEvSq);

//     printf("twoPartCorrEv = %.3f   twoPartCorrEvSq = %.3f \n\n",twoPartCorrEv,twoPartCorrEvSq);

    // -- End analysis in multiplicity bins --


    // -- Analysis in centrality bins --

    // Centrality V0
    if (fUseCentrality == 1) {

//       printf("CentV0bin: %i  ",centralityVZEROBin);

      if (centralityVZEROBin < 91) {

// 	printf("CentV0bin: %i   ",centralityVZEROBin);

	for (int i=0; i<nrTracks; i++) {
	  for (int j=i+1; j<nrTracks; j++) {
	    twoPartCorrPairCent = (tracks[i] - evMptCent) * (tracks[j] - evMptCent);
	    twoPartCorrEvCent = twoPartCorrEvCent + twoPartCorrPairCent;
	  }
	}

	twoPartCorrEvCentSq = twoPartCorrEvCent * twoPartCorrEvCent;
	fTwoPartCorrEvCent->Fill(centralityVZEROBin,twoPartCorrEvCent);
	fTwoPartCorrEvCentSq->Fill(centralityVZEROBin,twoPartCorrEvCentSq);

// 	printf("CentV0bin: %i   evMptCent: %.3f   CorrEvCent: %.2f   CorrEvCentSq: %.2f\n",centralityVZEROBin,evMptCent,twoPartCorrEvCent,twoPartCorrEvCentSq);
      }
    }

    // -- End analysis in centrality bins --

    // --- End two-particle correlator ---


  } // --- End calculation of various values and filling of histograms ---


  // Post output data
  PostData(1, fOutputList);
}

void AliAnalysisTaskPtFlucPbPb::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: fOutputList not available\n");
    return;
  }

}
Commit	Line	Data
befb009c	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TCanvas.h"
	6	#include "TF1.h"
	7	#include "TProfile.h"
	8	#include "TList.h"
	9	#include "iostream"
	10
	11	#include "AliAnalysisTaskSE.h"
	12	#include "AliAnalysisManager.h"
	13
	14	#include "AliESD.h"
	15	#include "AliESDEvent.h"
	16	#include "AliESDInputHandler.h"
	17	#include "AliESDtrackCuts.h"
	18	#include "AliCentrality.h"
	19
	20	#include "AliLog.h"
	21
	22	#include "AliAnalysisTaskPtFlucPbPb.h"
	23
	24
	25	using namespace std;
	26
	27	// Analysis of Pt FLuctuations (PbPb)
	28	// Author: Stefan Heckel
	29	// Version of PbPb task: 5.0, 18.04.2011
	30
	31
	32	ClassImp(AliAnalysisTaskPtFlucPbPb)
	33
	34	//________________________________________________________________________
	35	AliAnalysisTaskPtFlucPbPb::AliAnalysisTaskPtFlucPbPb(const char *name)
	36	:AliAnalysisTaskSE(name),
	37	fESD(0),
	38	fOutputList(0),
	39	fPtSpec(0),
	40	fMult(0),
	41	fMultSum(0),
	42	fMultSumPt(0),
	43	fMultNrPairs(0),
	44	fCent(0),
	45	fCentSum(0),
	46	fCentSumPt(0),
	47	fCentNrPairs(0),
	48	fEta(0),
	49	fEtaPhiPlus(0),
	50	fEtaPhiMinus(0),
	51	fVtxZ(0),
	52	fVtxZCut(0),
	53	fVtxZCont(0),
	54	fVtxZTrackCuts(0),
	55	fEventMeanPt(0),
	56	fEventMeanPtSq(0),
	57	fMultEventMeanPt(0),
	58	fMultEventMeanPtSq(0),
	59	fCentEventMeanPt(0),
	60	fCentEventMeanPtSq(0),
	61	fTwoPartCorrEv(0),
	62	fTwoPartCorrEvSq(0),
	63	fTwoPartCorrEvCent(0),
	64	fTwoPartCorrEvCentSq(0),
65	fESDTrackCuts(0),
66	fMaxVertexZ(0),
67	fNContributors(0),
68	fUseCentrality(0),
69	fMC(0)
70	{
71	DefineOutput(1, TList::Class());
72	}
73
74	//________________________________________________________________________
75	AliAnalysisTaskPtFlucPbPb::~AliAnalysisTaskPtFlucPbPb()
76	{
77	if(fOutputList) delete fOutputList; fOutputList =0;
78	}
79
80	//________________________________________________________________________
81	void AliAnalysisTaskPtFlucPbPb::UserCreateOutputObjects()
82	{
83	// Create histograms
84	// Called once
85
86	OpenFile(1, "RECREATE");
87	fOutputList = new TList();
88	fOutputList->SetOwner();
89
90
91	fPtSpec = new TH1F("fPtSpec","Pt spectrum",50,0,2.5);
92
93	fMult = new TH1F("fMult","Multiplicity distribution",30,0,3000);
94
95	fMultSum = new TH1F("fMultSum","Sum of nrtracks of all events in multiplicity bins",30,0,3000);
96
97	fMultSumPt = new TH1F("fMultSumPt","Sum of pTs of all events in multiplicity bins",30,0,3000);
98
99	fMultNrPairs = new TH1F("fMultNrPairs","Sum of number of pairs in multiplicity bins",30,0,3000);
100
101	fCent = new TH1F("fCent","Centrality distribution",21,0,105);
102
103	fCentSum = new TH1F("fCentSum","Sum of nrtracks of all events in centrality bins",21,0,105);
104
105	fCentSumPt = new TH1F("fCentSumPt","Sum of pTs of all events in centrality bins",21,0,105);
106
107	fCentNrPairs = new TH1F("fCentNrPairs","Sum of number of pairs in centrality bins",21,0,105);
108
109
110	fEta = new TH1F("fEta","Eta distribution",80,-2,2);
111
112	fEtaPhiPlus = new TH1F("fEtaPhiPlus","Phi distribution for positive eta",62,0,6.2);
113
114	fEtaPhiMinus = new TH1F("fEtaPhiMinus","Phi distribution for negative eta",62,0,6.2);
115
116	fVtxZ = new TH1F("fVtxZ","Vertex Z distribution before cuts",100,-20,20);
117
118	fVtxZCut = new TH1F("fVtxZCut","Vertex Z distribution after vtxZ cut",110,-11,11);
119
120	fVtxZCont = new TH1F("fVtxZCont","Vertex Z distribution after nCont cut",110,-11,11);
121
122	fVtxZTrackCuts = new TH1F("fVtxZTrackCuts","Vertex Z distribution after track cuts",110,-11,11);
123
124
125	fEventMeanPt = new TH1F("fEventMeanPt","Mean-Pt event by event",50,0,2.5);
126
127	fEventMeanPtSq = new TH1F("fEventMeanPtSq","Mean-Pt event by event squared",100,0,5);
128
129	fMultEventMeanPt = new TH1F("fMultEventMeanPt","Mean-Pt event by event vs. multiplicity",30,0,3000);
130
131	fMultEventMeanPtSq = new TH1F("fMultEventMeanPtSq","Mean-Pt event by event squared vs. multiplicity",30,0,3000);
132
133	fCentEventMeanPt = new TH1F("fCentEventMeanPt","Mean-Pt event by event vs. centrality",21,0,105);
134
135	fCentEventMeanPtSq = new TH1F("fCentEventMeanPtSq","Mean-Pt event by event squared vs. centrality",21,0,105);
136
137
138	fTwoPartCorrEv = new TH1F("fTwoPartCorrEv","Two-particle correlator vs. multiplicity",30,0,3000);
139
140	fTwoPartCorrEvSq = new TH1F("fTwoPartCorrEvSq","Two-particle correlator squared vs. multiplicity",30,0,3000);
141
142	fTwoPartCorrEvCent = new TH1F("fTwoPartCorrEvCent","Two-particle correlator vs. centrality",21,0,105);
143
144	fTwoPartCorrEvCentSq = new TH1F("fTwoPartCorrEvCentSq","Two-particle correlator squared vs. centrality",21,0,105);
145
146
147	// Add histograms to the output list
148	fOutputList->Add(fPtSpec);
149	fOutputList->Add(fMult);
150	fOutputList->Add(fMultSum);
151	fOutputList->Add(fMultSumPt);
152	fOutputList->Add(fMultNrPairs);
153	fOutputList->Add(fCent);
154	fOutputList->Add(fCentSum);
155	fOutputList->Add(fCentSumPt);
156	fOutputList->Add(fCentNrPairs);
157	fOutputList->Add(fEta);
158	fOutputList->Add(fEtaPhiPlus);
159	fOutputList->Add(fEtaPhiMinus);
160	fOutputList->Add(fVtxZ);
161	fOutputList->Add(fVtxZCut);
162	fOutputList->Add(fVtxZCont);
163	fOutputList->Add(fVtxZTrackCuts);
164	fOutputList->Add(fEventMeanPt);
165	fOutputList->Add(fEventMeanPtSq);
166	fOutputList->Add(fMultEventMeanPt);
167	fOutputList->Add(fMultEventMeanPtSq);
168	fOutputList->Add(fCentEventMeanPt);
169	fOutputList->Add(fCentEventMeanPtSq);
170	fOutputList->Add(fTwoPartCorrEv);
171	fOutputList->Add(fTwoPartCorrEvSq);
172	fOutputList->Add(fTwoPartCorrEvCent);
173	fOutputList->Add(fTwoPartCorrEvCentSq);
174	}
175
176	//________________________________________________________________________
177	void AliAnalysisTaskPtFlucPbPb::UserExec(Option_t *)
178	{
179	// Main loop
180	// Called for each event
181
182
183	// --- ESD event handler ---
184
185	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
186
187	if (!esdH) {
188	Printf("ERROR: Could not get ESDInputHandler");
189	}
190	else {
191	fESD = esdH->GetEvent();
192	}
193
194	if (!fESD) {
195	Printf("ERROR: fESD not available");
196	return;
197	}
198
befb009c	199
	200	// --- End event handler ---
	201
	202
	203	// --- Vertex cuts ---
	204
	205	// // TPC+ITS vertex
	206	// const AliESDVertex* vtxESD = fESD->GetPrimaryVertexTracks();
	207	// Double_t vtxZ = vtxESD->GetZv();
	208	// Double_t vtxNCont = vtxESD->GetNContributors();
	209
	210	// TPConly vertex
	211	const AliESDVertex* vtxESDTPC = fESD->GetPrimaryVertexTPC();
	212	Double_t vtxZ = vtxESDTPC->GetZv();
	213	Double_t vtxNCont = vtxESDTPC->GetNContributors();
	214
	215	fVtxZ->Fill(vtxZ); // VtxZ before cuts
	216
	217	// Event cut on the z-Position of the vertex
	218	if(vtxZ > fMaxVertexZ \|\| vtxZ < (-1.*fMaxVertexZ)) {
	219	// Printf("VertexZ out of range, Zv = %f",vtxZ);
	220	return;
	221	}
	222
	223	fVtxZCut->Fill(vtxZ); // VtxZ after cut on vtxZ
	224
	225	// Event cut on the number of contributors
	226	if(vtxNCont < fNContributors) {
	227	// Printf("Vertex has no contributors");
	228	return;
	229	}
	230	// else {
	231	// Printf("Vertex nContributors = %i",vtxNCont);
	232	// }
	233
	234	fVtxZCont->Fill(vtxZ); // VtxZ after cut on nContributors
	235
	236	// --- End vertex cuts ---
	237
	238
	239	Int_t nrTracks = 0; // count for number of tracks which pass the track cuts
	240	Int_t nrESDTracks = 0; // number of tracks in the ESD file
	241	Double_t sumPt = 0.; // sum of track Pts
	242	Double_t twoPartCorrEv = 0.; // Two-particle correlator of one event for multiplicity bin analysis
	243	Double_t twoPartCorrEvCent = 0.; // Two-particle correlator of one event for centrality bin analysis
	244
	245	Double_t tracks[3000] = {0.}; // array of track Pts, needed for the two-particle correlator
	246
	247	Double_t trackPt = 0., trackEta = 0., trackPhi = 0.;
	248	Double_t eventMeanPt = 0., eventMeanPtSq = 0., evMptMult = 0.;
	249	Double_t twoPartCorrPair = 0., twoPartCorrEvSq = 0.;
	250	Double_t twoPartCorrPairCent = 0., twoPartCorrEvCentSq = 0.;
	251	Double_t nrPairs = 0.;
	252
	253	Double_t *nbins = 0x0; // Mean pT values for multiplicity bin analysis
	254	Double_t *centbins = 0x0; // Mean pT values for centrality bin analysis
	255
	256	Double_t centralityVZERO = 0.;
e6a162b1	257	Int_t centralityVZEROBin = 0;
befb009c	258
e6a162b1	259	Int_t centBin = 0;
befb009c	260	Double_t evMptCent = 0.;
	261
	262
	263	// Printf("MC: %d",fMC);
	264
	265
	266	// --- Mean pT values ---
	267
	268	// !!!!! Have to be calculated in a first run - for each sample, which should be analysed, separately! !!!!!
	269
	270
	271	// -- Mean pT values for multiplicity bins (first bin is for nrTracks = 0 and has to be 0) --
	272
	273	// MC PbPb 2.76 ATeV 11a10a (Hijing) (Pt-Range: 0.15 - 2)
	274	Double_t nbinsMC276[32] = {0.000, 0.525, 0.535, 0.539, 0.541, 0.543, 0.544, 0.545, 0.545, 0.546, 0.547, 0.548, 0.548, 0.549, 0.549, 0.550, 0.551, 0.551, 0.552, 0.552, 0.553, 0.554, 0.554, 0.555, 0.556, 0.556, 0.557, 0.558, 0.549, 0.000, 0.000, 0.000};
	275
	276	// Data PbPb 2.76 ATeV 10h.pass1 (Pt-Range: 0.15 - 2)
	277	Double_t nbinsData276[32] = {0.000, 0.580, 0.605, 0.618, 0.626, 0.632, 0.636, 0.640, 0.642, 0.645, 0.647, 0.648, 0.649, 0.650, 0.651, 0.652, 0.653, 0.653, 0.654, 0.654, 0.654, 0.654, 0.654, 0.654, 0.655, 0.656, 0.658, 0.660, 0.661, 0.659, 0.660, 0.656};
	278
	279	// -- End mean pT values for multiplicity bins --
	280
	281
	282	// -- Mean pT values for centrality bins --
	283
	284	// MC PbPb 2.76 ATeV 11a10a (Hijing) (Pt-Range: 0.15 - 2)
	285	Double_t centbinsMC276[21] = {0.555, 0.553, 0.551, 0.549, 0.548, 0.547, 0.545, 0.544, 0.543, 0.541, 0.539, 0.537, 0.536, 0.533, 0.530, 0.528, 0.524, 0.518, 0.505, 0.000, 0.534};
	286
	287	// Data PbPb 2.76 ATeV 10h.pass1 (Pt-Range: 0.15 - 2)
	288	Double_t centbinsData276[21] = {0.654, 0.654, 0.652, 0.650, 0.647, 0.644, 0.640, 0.635, 0.629, 0.623, 0.616, 0.609, 0.601, 0.594, 0.586, 0.579, 0.568, 0.551, 0.535, 0.000, 0.621};
	289
	290	// -- End mean pT values for centrality bins --
	291
	292
	293	// -- Selection of MC/Data; whole sample values --
	294
	295	if (fMC) { // - MC -
	296
	297	// Printf(" -- MC, 2.76 ATeV -- ");
	298
	299	nbins = nbinsMC276;
	300	centbins = centbinsMC276;
	301
	302	} // - End MC -
	303	else { // - Data -
	304
	305	// Printf(" -- Data, 2.76 ATeV -- ");
	306
	307	nbins = nbinsData276;
	308	centbins = centbinsData276;
	309
	310	} // - End data -
	311
	312	// -- End selection of MC/Data; whole sample values --
	313
	314	// --- End mean pT values ---
	315
	316
	317	nrESDTracks = fESD->GetNumberOfTracks();
	318	// if( nrESDTracks ) { printf("Found event with %i tracks.\n",nrESDTracks); }
	319
	320
	321	if(fUseCentrality != 0) {
	322
	323	AliCentrality *esdCentrality = fESD->GetCentrality();
324
325	//V0
326	if (fUseCentrality == 1){
327
328	centralityVZERO = esdCentrality->GetCentralityPercentile("V0M");
329
330	if ( centralityVZERO > 0. && centralityVZERO < 5.) centralityVZEROBin = 0;
331	else if ( centralityVZERO >= 5. && centralityVZERO < 10.) centralityVZEROBin = 5;
332	else if ( centralityVZERO >= 10. && centralityVZERO < 15.) centralityVZEROBin = 10;
333	else if ( centralityVZERO >= 15. && centralityVZERO < 20.) centralityVZEROBin = 15;
334	else if ( centralityVZERO >= 20. && centralityVZERO < 25.) centralityVZEROBin = 20;
335	else if ( centralityVZERO >= 25. && centralityVZERO < 30.) centralityVZEROBin = 25;
336	else if ( centralityVZERO >= 30. && centralityVZERO < 35.) centralityVZEROBin = 30;
337	else if ( centralityVZERO >= 35. && centralityVZERO < 40.) centralityVZEROBin = 35;
338	else if ( centralityVZERO >= 40. && centralityVZERO < 45.) centralityVZEROBin = 40;
339	else if ( centralityVZERO >= 45. && centralityVZERO < 50.) centralityVZEROBin = 45;
340	else if ( centralityVZERO >= 50. && centralityVZERO < 55.) centralityVZEROBin = 50;
341	else if ( centralityVZERO >= 55. && centralityVZERO < 60.) centralityVZEROBin = 55;
342	else if ( centralityVZERO >= 60. && centralityVZERO < 65.) centralityVZEROBin = 60;
343	else if ( centralityVZERO >= 65. && centralityVZERO < 70.) centralityVZEROBin = 65;
344	else if ( centralityVZERO >= 70. && centralityVZERO < 75.) centralityVZEROBin = 70;
345	else if ( centralityVZERO >= 75. && centralityVZERO < 80.) centralityVZEROBin = 75;
346	else if ( centralityVZERO >= 80. && centralityVZERO < 85.) centralityVZEROBin = 80;
347	else if ( centralityVZERO >= 85. && centralityVZERO < 90.) centralityVZEROBin = 85;
348	else if ( centralityVZERO >= 90. && centralityVZERO < 95.) centralityVZEROBin = 90;
349	else if ( centralityVZERO >= 95. && centralityVZERO < 99.) centralityVZEROBin = 95;
350	else if ( centralityVZERO >= 99. ) centralityVZEROBin = 100;
351	else if ( centralityVZERO <= 0. ) centralityVZEROBin = 100;
352
353	centBin = centralityVZEROBin / 5;
354	evMptCent = centbins[centBin];
355	}
356
357	// Printf("\nRaw mult: %i CentVZERO: %f CentVZERObin: %i centBin: %i evMptCent: %.3f ",nrESDTracks,centralityVZERO,centralityVZEROBin,centBin,evMptCent);
358
359	}
360
361
362	// --- Loop over all tracks of one event ---
363
364	for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
365	AliESDtrack* track = fESD->GetTrack(iTracks);
366	if (!track) {
367	Printf("ERROR: Could not receive track %d\n", iTracks);
368	continue;
369	}
e6a162b1	370	if(!fESDTrackCuts) {
	371	Printf("ERROR: No esd track cut");
	372	continue;
	373	}
	374	else {
	375	if(!fESDTrackCuts->AcceptTrack(track))
	376	continue;
	377	}
befb009c	378	trackPt = track->Pt();
	379	fPtSpec->Fill(trackPt);
	380	tracks[nrTracks] = trackPt;
	381
	382	trackEta = track->Eta();
	383	fEta->Fill(trackEta);
	384
	385	trackPhi = track->Phi();
	386
	387	if (trackEta > 0.) {
	388	fEtaPhiPlus->Fill(trackPhi);
	389	}
	390	else if (trackEta < 0.) {
	391	fEtaPhiMinus->Fill(trackPhi);
	392	}
	393
	394	sumPt = sumPt + trackPt;
	395	nrTracks++;
	396
	397	// printf("Track Pt = %.3f Track Eta = %.3f Track Phi = %3f\n",trackPt,trackEta,trackPhi);
	398
	399	} // --- End track loop ---
	400
	401
	402	// --- Calculation of various values and filling of histograms (for remaining events with N_acc > 0) ---
	403
	404	if(nrTracks != 0) {
	405
	406	// VtxZ after all track cuts
	407	fVtxZTrackCuts->Fill(vtxZ);
	408
	409	// Multiplicity distributions
	410	fMult->Fill(nrTracks);
	411	fMultSum->Fill(nrTracks,nrTracks);
	412
	413	// Number of pairs in event
	414	nrPairs = 0.5 * nrTracks * (nrTracks-1);
	415	fMultNrPairs->Fill(nrTracks,nrPairs);
	416
	417	// Calculation of mean Pt and mean Pt Squared
	418	eventMeanPt = sumPt / nrTracks;
	419	eventMeanPtSq = eventMeanPt * eventMeanPt;
	420
	421	// Mean-Pt and Mean-Pt squared
	422	fEventMeanPt->Fill(eventMeanPt);
	423	fEventMeanPtSq->Fill(eventMeanPtSq);
	424
	425	// Mean-Pt and Mean-Pt squared depending on multiplicity
	426	fMultEventMeanPt->Fill(nrTracks,eventMeanPt);
	427	fMultEventMeanPtSq->Fill(nrTracks,eventMeanPtSq);
	428	fMultSumPt->Fill(nrTracks,sumPt);
	429
	430	// printf("nrTracks: %i sumPt: %.8f meanPt: %.8f meanPtSq: %.8f\n",nrTracks,sumPt,eventMeanPt,eventMeanPtSq);
	431
	432
	433	// Centrality V0
	434	if (fUseCentrality == 1) {
	435
	436	// Centrality distributions
	437	fCent->Fill(centralityVZEROBin);
	438	fCentSum->Fill(centralityVZEROBin,nrTracks);
	439
	440	// Number of pairs in event
	441	fCentNrPairs->Fill(centralityVZEROBin,nrPairs);
442
443	// Mean-Pt and Mean-Pt squared depending on centrality
444	fCentEventMeanPt->Fill(centralityVZEROBin,eventMeanPt);
445	fCentEventMeanPtSq->Fill(centralityVZEROBin,eventMeanPtSq);
446	fCentSumPt->Fill(centralityVZEROBin,sumPt);
447
448	// printf("CentV0bin: %i NrTracks: %i NrPairs: %.0f SumPt: %.1f\n",centralityVZEROBin,nrTracks,nrPairs,sumPt);
449	}
450
451
452	// --- Two-particle correlator ---
453
454	// -- Analysis in multiplicity bins --
455
456	for (int k=1; k<32; k++) {
457	if (nrTracks > 100 * (k-1) && nrTracks <= 100 * k) {
458	evMptMult = nbins[k];
459	break;
460	}
461	}
462	// printf("nrTracks = %3i evMptMult = %.3f \n",nrTracks,evMptMult);
463
464	for (int i=0; i<nrTracks; i++) {
465	for (int j=i+1; j<nrTracks; j++) {
466	twoPartCorrPair = (tracks[i] - evMptMult) * (tracks[j] - evMptMult);
467	twoPartCorrEv = twoPartCorrEv + twoPartCorrPair;
468	}
469	}
470
471	twoPartCorrEvSq = twoPartCorrEv * twoPartCorrEv;
472	fTwoPartCorrEv->Fill(nrTracks,twoPartCorrEv);
473	fTwoPartCorrEvSq->Fill(nrTracks,twoPartCorrEvSq);
474
475	// printf("twoPartCorrEv = %.3f twoPartCorrEvSq = %.3f \n\n",twoPartCorrEv,twoPartCorrEvSq);
476
477	// -- End analysis in multiplicity bins --
478
479
480	// -- Analysis in centrality bins --
481
482	// Centrality V0
483	if (fUseCentrality == 1) {
484
485	// printf("CentV0bin: %i ",centralityVZEROBin);
486
487	if (centralityVZEROBin < 91) {
488
489	// printf("CentV0bin: %i ",centralityVZEROBin);
490
491	for (int i=0; i<nrTracks; i++) {
492	for (int j=i+1; j<nrTracks; j++) {
493	twoPartCorrPairCent = (tracks[i] - evMptCent) * (tracks[j] - evMptCent);
494	twoPartCorrEvCent = twoPartCorrEvCent + twoPartCorrPairCent;
495	}
496	}
497
498	twoPartCorrEvCentSq = twoPartCorrEvCent * twoPartCorrEvCent;
499	fTwoPartCorrEvCent->Fill(centralityVZEROBin,twoPartCorrEvCent);
500	fTwoPartCorrEvCentSq->Fill(centralityVZEROBin,twoPartCorrEvCentSq);
501
502	// printf("CentV0bin: %i evMptCent: %.3f CorrEvCent: %.2f CorrEvCentSq: %.2f\n",centralityVZEROBin,evMptCent,twoPartCorrEvCent,twoPartCorrEvCentSq);
503	}
504	}
505
506	// -- End analysis in centrality bins --
507
508	// --- End two-particle correlator ---
509
510
511	} // --- End calculation of various values and filling of histograms ---
512
513
514	// Post output data
515	PostData(1, fOutputList);
516	}
517
518	void AliAnalysisTaskPtFlucPbPb::Terminate(Option_t *)
519	{
520	// Draw result to the screen
521	// Called once at the end of the query
522
523	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
524	if (!fOutputList) {
525	Printf("ERROR: fOutputList not available\n");
526	return;
527	}
528
529	}