PWGCF/EBYE/MeanPtFluctuations/AliAnalysisTaskPtFluc.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TList.h"
   6 #include "TRandom3.h"
   7 #include "iostream"
   8
   9 #include "AliAnalysisTaskSE.h"
  10 #include "AliAnalysisManager.h"
  11
  12 #include "AliESDEvent.h"
  13 #include "AliESDInputHandler.h"
  14 #include "AliESDtrackCuts.h"
  15
  16 #include "AliMCEvent.h"
  17 #include "AliStack.h"
  18 #include "AliMCEventHandler.h"
  19
  20 #include "AliAnalysisTaskPtFluc.h"
  21
  22
  23 using namespace std;
  24
  25 // Analysis of Pt Fluctuations (pp)
  26 // Author: Stefan Heckel
  27 // Version of pp task:   11.4, 11.06.2012
  28
  29
  30 ClassImp(AliAnalysisTaskPtFluc)
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskPtFluc::AliAnalysisTaskPtFluc(const char *name)
  34   :AliAnalysisTaskSE(name),
  35   fESD(0),
  36   fMCev(0),
  37   fRandom3(0),
  38   fOutputList(0),
  39   fPtSpec(0),
  40   fPtSpec2(0),
  41   fMult(0),
  42   fEta(0),
  43   fEtaPhiPlus(0),
  44   fEtaPhiMinus(0),
  45   fVtxZ(0),
  46   fVtxZCut(0),
  47   fVtxZCont(0),
  48   fVtxZCutDiff(0),
  49   fVtxZPileup(0),
  50   fVtxZTrackCuts(0),
  51   fVtxZDiff1(0),
  52   fVtxZDiff2(0),
  53   fVtxZDiff3(0),
  54   fVtxZDiff1b(0),
  55   fVtxZDiff2b(0),
  56   fVtxZDiff3b(0),
  57   fEventMeanPt(0),
  58   fEventMeanPtSq(0),
  59   fEventMeanPtMult(0),
  60   fMultEventMeanPt(0),
  61   fMultEventMeanPtSq(0),
  62   fTwoPartCorrEv(0),
  63   fTwoPartCorrEvSq(0),
  64   fTwoPartCorrEvSample(0),
  65   fTwoPartCorrEvSampleSq(0),
  66   fESDTrackCuts(0),
  67   fMaxVertexZ(0),
  68   fMaxVertexZDiff1(0),
  69   fNContributors(0),
  70   fMC(0),
  71   fMCType(0)
  72 {
  73   DefineOutput(1, TList::Class());
  74 }
  75
  76 //________________________________________________________________________
  77 AliAnalysisTaskPtFluc::~AliAnalysisTaskPtFluc()
  78 {
  79   if (fRandom3) delete fRandom3; fRandom3 = 0;
  80
  81   if (fOutputList) {
  82     fOutputList->Clear();
  83     delete fOutputList;
  84   }
  85   fOutputList = 0;
  86
  87   if (fESDTrackCuts) delete fESDTrackCuts; fESDTrackCuts = 0;
  88
  89 }
  90
  91 //________________________________________________________________________
  92 void AliAnalysisTaskPtFluc::UserCreateOutputObjects()
  93 {
  94   // Create histograms
  95   // Called once
  96
  97   OpenFile(1, "RECREATE");
  98   fOutputList = new TList();
  99   fOutputList->SetOwner();
 100
 101
 102   fPtSpec = new TH1F("fPtSpec","Pt spectrum",100,0,2.5);
 103   fPtSpec2 = new TH1F("fPtSpec2","Pt spectrum 2 - MC ESD",100,0,2.5);
 104   fMult = new TH1F("fMult","Multiplicity distribution",120,-0.5,119.5);
 105
 106   fEta = new TH1F("fEta","Eta distribution",80,-2,2);
 107   fEtaPhiPlus = new TH1F("fEtaPhiPlus","Phi distribution for positive eta",62,0,6.2);
 108   fEtaPhiMinus = new TH1F("fEtaPhiMinus","Phi distribution for negative eta",62,0,6.2);
 109
 110   fVtxZ = new TH1F("fVtxZ","Vertex Z distribution before cuts",100,-20,20);
 111   fVtxZCut = new TH1F("fVtxZCut","Vertex Z distribution after vtxZ cut",110,-11,11);
 112   fVtxZCont = new TH1F("fVtxZCont","Vertex Z distribution after nCont cut",110,-11,11);
 113   fVtxZCutDiff = new TH1F("fVtxZCutDiff","Vertex Z distribution after cut on vtx Z Diff",110,-11,11);
 114   fVtxZPileup = new TH1F("fVtxZPileup","Vertex Z distribution after pileup rejection",110,-11,11);
 115   fVtxZTrackCuts = new TH1F("fVtxZTrackCuts","Vertex Z distribution after track cuts",110,-11,11);
 116
 117   fVtxZDiff1 = new TH1F("fVtxZDiff1","Difference 1 between vertex Z distributions",100,-5,5);
 118   fVtxZDiff2 = new TH1F("fVtxZDiff2","Difference 2 between vertex Z distributions",100,-5,5);
 119   fVtxZDiff3 = new TH1F("fVtxZDiff3","Difference 3 between vertex Z distributions",100,-5,5);
 120   fVtxZDiff1b = new TH1F("fVtxZDiff1b","Difference 1 between vertex Z distributions after all cuts",100,-5,5);
 121   fVtxZDiff2b = new TH1F("fVtxZDiff2b","Difference 2 between vertex Z distributions after all cuts",100,-5,5);
 122   fVtxZDiff3b = new TH1F("fVtxZDiff3b","Difference 3 between vertex Z distributions after all cuts",100,-5,5);
 123
 124
 125   fEventMeanPt = new TH1F("fEventMeanPt","Mean-Pt distribution",50,0,2.5);
 126   fEventMeanPtSq = new TH1F("fEventMeanPtSq","Mean-Pt squared distribution",100,0,5);
 127   fEventMeanPtMult = new TH2F("fEventMeanPtMult","Mean-Pt for single events vs. multiplicity",20,0,100,200,0.,2.);
 128
 129   fMultEventMeanPt = new TH1F("fMultEventMeanPt","Mean-Pt event by event vs. Multiplicity",100,-0.5,99.5);
 130   fMultEventMeanPtSq = new TH1F("fMultEventMeanPtSq","Mean-Pt event by event squared vs. Multiplicity",100,-0.5,99.5);
 131
 132
 133   fTwoPartCorrEv = new TH1F("fTwoPartCorrEv","Two-particle correlator vs. Multiplicity",100,-0.5,99.5);
 134   fTwoPartCorrEvSq = new TH1F("fTwoPartCorrEvSq","Two-particle correlator squared vs. Multiplicity",100,-0.5,99.5);
 135   fTwoPartCorrEvSample = new TH1F("fTwoPartCorrEvSample","Two-particle correlator vs. Multiplicity",100,-0.5,99.5);
 136   fTwoPartCorrEvSampleSq = new TH1F("fTwoPartCorrEvSampleSq","Two-particle correlator squared vs. Multiplicity",100,-0.5,99.5);
 137
 138
 139   // Add histograms to the output list
 140   fOutputList->Add(fPtSpec);
 141   fOutputList->Add(fPtSpec2);
 142   fOutputList->Add(fMult);
 143   fOutputList->Add(fEta);
 144   fOutputList->Add(fEtaPhiPlus);
 145   fOutputList->Add(fEtaPhiMinus);
 146   fOutputList->Add(fVtxZ);
 147   fOutputList->Add(fVtxZCut);
 148   fOutputList->Add(fVtxZCont);
 149   fOutputList->Add(fVtxZCutDiff);
 150   fOutputList->Add(fVtxZPileup);
 151   fOutputList->Add(fVtxZTrackCuts);
 152   fOutputList->Add(fVtxZDiff1);
 153   fOutputList->Add(fVtxZDiff2);
 154   fOutputList->Add(fVtxZDiff3);
 155   fOutputList->Add(fVtxZDiff1b);
 156   fOutputList->Add(fVtxZDiff2b);
 157   fOutputList->Add(fVtxZDiff3b);
 158   fOutputList->Add(fEventMeanPt);
 159   fOutputList->Add(fEventMeanPtSq);
 160   fOutputList->Add(fEventMeanPtMult);
 161   fOutputList->Add(fMultEventMeanPt);
 162   fOutputList->Add(fMultEventMeanPtSq);
 163   fOutputList->Add(fTwoPartCorrEv);
 164   fOutputList->Add(fTwoPartCorrEvSq);
 165   fOutputList->Add(fTwoPartCorrEvSample);
 166   fOutputList->Add(fTwoPartCorrEvSampleSq);
 167
 168   // Post output data (if histograms are not used later, PostData is at least called here)
 169   PostData(1, fOutputList);
 170
 171 }
 172
 173 //________________________________________________________________________
 174 void AliAnalysisTaskPtFluc::UserExec(Option_t *)
 175 {
 176   // Main loop
 177   // Called for each event
 178
 179
 180   // --- ESD event handler ---
 181   if (!fMC || fMCType < 3) {
 182
 183 //     Printf("ESD 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
 199     if(!fESDTrackCuts) Printf("ERROR: No esd track cut");
 200
 201   } // --- End ESD event handler ---
 202
 203
 204   // --- MC event handler ---
 205
 206   AliStack *stack = NULL;
 207
 208   if (fMC) {
 209
 210 //     Printf("MC handler");
 211
 212     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 213
 214     if (!mcH) {
 215       Printf("ERROR: Could not get MCInputHandler");
 216       return;
 217     }
 218     else {
 219      fMCev = mcH->MCEvent();
 220     }
 221
 222     if (!fMCev) {
 223       Printf("ERROR: fMCev not available");
 224       return;
 225     }
 226
 227    stack = fMCev->Stack();
 228
 229   }
 230
 231   // --- End MC event handler ---
 232
 233
 234   // --- Vertex cuts ---
 235
 236   Double_t vtxZGlobal=0., vtxZSPD=0., vtxZTPC=0.;
 237   Double_t vtxNContGlobal=0.;
 238 //   Double_t  vtxNContSPD=0., vtxNContTPC=0.;
 239   Double_t vtxZ=0., vtxNCont=0.;
 240   Double_t vtxZdiff1=0., vtxZdiff2=0., vtxZdiff3=0.;
 241
 242  if (!fMC || fMCType < 3) {
 243
 244   // Global vertex
 245   const AliESDVertex* vtxESD = fESD->GetPrimaryVertexTracks();
 246   vtxZGlobal = vtxESD->GetZ();
 247   vtxNContGlobal = vtxESD->GetNContributors();
 248
 249   // SPD vertex
 250   const AliESDVertex* vtxESDSPD = fESD->GetPrimaryVertexSPD();
 251   vtxZSPD = vtxESDSPD->GetZ();
 252 //   vtxNContSPD = vtxESDSPD->GetNContributors();
 253
 254   // TPC vertex
 255   const AliESDVertex* vtxESDTPC = fESD->GetPrimaryVertexTPC();
 256   vtxZTPC = vtxESDTPC->GetZ();
 257 //   vtxNContTPC = vtxESDTPC->GetNContributors();
 258
 259   vtxZ = vtxZGlobal;
 260   vtxNCont = vtxNContGlobal;
 261
 262   fVtxZ->Fill(vtxZ); // VtxZ before cuts
 263
 264   // Differences between vertices
 265   vtxZdiff1 = vtxZTPC - vtxZGlobal;
 266   vtxZdiff2 = vtxZTPC - vtxZSPD;
 267   vtxZdiff3 = vtxZGlobal - vtxZSPD;
 268   fVtxZDiff1->Fill(vtxZdiff1);
 269   fVtxZDiff2->Fill(vtxZdiff2);
 270   fVtxZDiff3->Fill(vtxZdiff3);
 271
 272   // Event cut on the z-position of the vertex
 273   if(vtxZ > fMaxVertexZ || vtxZ < (-1.*fMaxVertexZ)) {
 274 //     Printf("VertexZ out of range, Zv = %f",vtxZ);
 275     return;
 276   }
 277
 278   fVtxZCut->Fill(vtxZ); // VtxZ after cut on vtxZ
 279
 280   // Event cut on the number of contributors
 281   if(vtxNCont < fNContributors) {
 282 //     Printf("Vertex has no contributors");
 283     return;
 284   }
 285
 286   fVtxZCont->Fill(vtxZ); // VtxZ after cut on nContributors
 287
 288   // Event cut on the difference of the z-positions of the vertices (TPC - global)
 289   if (fMaxVertexZDiff1 > 0.) {
 290     if(vtxZdiff1 > fMaxVertexZDiff1 || vtxZdiff1 < (-1.*fMaxVertexZDiff1)) {
 291 //       Printf("VertexZ Diff (TPC - global) out of range, vtxZdiff1 = %f",vtxZ);
 292       return;
 293     }
 294   }
 295
 296   fVtxZCutDiff->Fill(vtxZ); // VtxZ after cut on vtxZDiff
 297
 298   } // --- End vertex cuts ---
 299
 300
 301   // --- Pileup rejection ---
 302   if (!fMC || fMCType < 3) {
 303
 304     if (fESD->IsPileupFromSPD(3,0.8,3.,2.,5.)) {
 305 //       Printf("IsPileupFromSPD");
 306       return;
 307     }
 308
 309     fVtxZPileup->Fill(vtxZ); // VtxZ after pileup rejection
 310
 311   } // --- End pileup rejection ---
 312
 313
 314   Int_t nrTracks = 0;           // count for number of tracks which pass the track cuts
 315 //   Int_t nrTracks2 = 0;               // count for number of tracks which pass the track cuts - for MC ESD vs. MC truth
 316   Int_t nrESDTracks = 0;        // number of tracks in the ESD file
 317   Int_t nrMCTracks = 0;         // number of tracks in the MC file
 318   Double_t sumPt = 0.;          // sum of track Pts
 319   Double_t twoPartCorrEv = 0.;  // Two-particle correlator of one event for multiplicity bin analysis
 320   Double_t twoPartCorrEvSample = 0.; // Two-part. corr. of one event for whole sample analysis
 321
 322   Double_t tracks[120] = {0.};  // array of track Pts, needed for the two-particle correlator
 323
 324   Double_t trackPt=0., trackEta=0., trackPhi=0.;
 325   Double_t trackPt2=0.; // for MC ESD
 326   Double_t eventMeanPt=0., eventMeanPtSq=0., evMptMult=0.;
 327   Double_t twoPartCorrPair=0., twoPartCorrEvSq=0.;
 328   Double_t twoPartCorrPairSample=0., twoPartCorrEvSampleSq=0.;
 329 //   Double_t nrPairs=0.;
 330
 331   Double_t *nbins = 0x0; // Mean pT values for multiplicity bin analysis
 332   Double_t evMptSample=0.;  // Mean pT value for whole sample analysis
 333
 334   Double_t random=0., trackPt40=0., ptRatio=0.;
 335   Int_t trackPtBin=0;
 336
 337   fRandom3 = new TRandom3();
 338
 339   if (fMCType == 2 || fMCType == 4) { // MC, Type = mod. MC truth
 340    fRandom3->SetSeed(0);
 341   }
 342
 343
 344   Double_t energy = 0.; // beam energy ---> has to be set by the user (in GeV) for MCType >= 3 (w/o ESD) <---
 345
 346   if (!fMC || fMCType < 3) {
 347     energy = fESD->GetBeamEnergy();
 348   }
 349
 350
 351 // --- Mean pT values ---
 352
 353 // !!!!! Have to be calculated in a first run - for each sample, which should be analysed, separately! !!!!!
 354
 355
 356 // -- Mean pT values for the multiplicity bins (starting with N = 0, therefore the first value has to be 0.000) --
 357
 358 // MC 900 GeV 10e13 (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
 359 Double_t nbinsMC900[100] = {0.000, 0.452, 0.468, 0.482, 0.497, 0.511, 0.524, 0.535, 0.544, 0.550, 0.557, 0.563, 0.567, 0.571, 0.576, 0.580, 0.583, 0.586, 0.590, 0.592, 0.595, 0.599, 0.600, 0.602, 0.604, 0.611, 0.604, 0.618, 0.614, 0.612, 0.632, 0.623, 0.617, 0.612, 0.640, 0.632, 0.649, 0.609, 0.658, 0.615, 0.677, 0.614, 0.646, 0.000, 0.000, 0.000, 0.000, 0.000, 0.672, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 360
 361 // MC 2.76 TeV 11b10a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
 362 Double_t nbinsMC276[100] = {0.000, 0.451, 0.467, 0.485, 0.505, 0.523, 0.540, 0.551, 0.563, 0.569, 0.577, 0.583, 0.588, 0.594, 0.598, 0.603, 0.608, 0.611, 0.615, 0.617, 0.622, 0.622, 0.626, 0.630, 0.633, 0.632, 0.637, 0.633, 0.641, 0.643, 0.650, 0.649, 0.655, 0.650, 0.669, 0.664, 0.653, 0.671, 0.668, 0.683, 0.693, 0.662, 0.765, 0.643, 0.677, 0.660, 0.000, 0.689, 0.708, 0.653, 0.635, 0.000, 0.719, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 363
 364 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
 365 Double_t nbinsMC7[100] = {0.000, 0.447, 0.465, 0.485, 0.508, 0.530, 0.548, 0.562, 0.573, 0.581, 0.589, 0.595, 0.601, 0.607, 0.612, 0.617, 0.621, 0.625, 0.629, 0.633, 0.637, 0.640, 0.643, 0.646, 0.649, 0.652, 0.654, 0.657, 0.659, 0.662, 0.664, 0.666, 0.668, 0.670, 0.672, 0.674, 0.677, 0.677, 0.680, 0.682, 0.684, 0.685, 0.687, 0.688, 0.693, 0.692, 0.692, 0.693, 0.696, 0.698, 0.696, 0.700, 0.702, 0.697, 0.701, 0.708, 0.712, 0.714, 0.705, 0.717, 0.702, 0.701, 0.716, 0.716, 0.761, 0.706, 0.750, 0.715, 0.731, 0.728, 0.711, 0.733, 0.000, 0.714, 0.704, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.652, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 366 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> values for MC truth <---
 367 // Double_t nbinsMC7[100] = {0.000, 0.441, 0.455, 0.473, 0.493, 0.515, 0.534, 0.550, 0.561, 0.570, 0.578, 0.585, 0.591, 0.597, 0.602, 0.607, 0.612, 0.616, 0.621, 0.624, 0.628, 0.631, 0.635, 0.638, 0.641, 0.644, 0.646, 0.649, 0.652, 0.654, 0.656, 0.659, 0.661, 0.663, 0.665, 0.666, 0.669, 0.670, 0.672, 0.674, 0.676, 0.678, 0.679, 0.683, 0.683, 0.685, 0.686, 0.686, 0.692, 0.693, 0.690, 0.691, 0.690, 0.692, 0.701, 0.701, 0.700, 0.704, 0.701, 0.704, 0.709, 0.698, 0.707, 0.711, 0.692, 0.739, 0.689, 0.715, 0.691, 0.720, 0.790, 0.670, 0.000, 0.717, 0.833, 0.716, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.619, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 368
 369
 370 // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> TPC standalone <------> values for MC ESD <---
 371 // Double_t nbinsMC7[100] = {0.000, 0.483, 0.490, 0.502, 0.515, 0.527, 0.535, 0.541, 0.545, 0.548, 0.552, 0.555, 0.558, 0.561, 0.563, 0.566, 0.568, 0.571, 0.573, 0.574, 0.576, 0.578, 0.579, 0.581, 0.582, 0.583, 0.584, 0.585, 0.585, 0.586, 0.586, 0.587, 0.588, 0.589, 0.588, 0.588, 0.588, 0.588, 0.587, 0.589, 0.588, 0.588, 0.587, 0.587, 0.586, 0.583, 0.588, 0.587, 0.589, 0.583, 0.586, 0.582, 0.584, 0.583, 0.585, 0.582, 0.587, 0.576, 0.580, 0.588, 0.566, 0.581, 0.594, 0.585, 0.531, 0.000, 0.580, 0.499, 0.534, 0.534, 0.580, 0.432, 0.674, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.655, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 372 // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> values for MC truth <---
 373 // Double_t nbinsMC7[100] = {0.000, 0.479, 0.483, 0.493, 0.506, 0.518, 0.527, 0.533, 0.537, 0.540, 0.543, 0.546, 0.549, 0.551, 0.554, 0.556, 0.558, 0.560, 0.562, 0.564, 0.566, 0.568, 0.569, 0.570, 0.571, 0.572, 0.573, 0.574, 0.575, 0.576, 0.576, 0.576, 0.577, 0.577, 0.577, 0.579, 0.577, 0.577, 0.577, 0.576, 0.577, 0.577, 0.577, 0.575, 0.575, 0.575, 0.573, 0.575, 0.571, 0.575, 0.568, 0.569, 0.570, 0.576, 0.566, 0.570, 0.568, 0.575, 0.565, 0.562, 0.564, 0.545, 0.563, 0.527, 0.551, 0.540, 0.571, 0.561, 0.562, 0.553, 0.000, 0.586, 0.000, 0.603, 0.480, 0.489, 0.566, 0.513, 0.715, 0.000, 0.000, 0.000, 0.441, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 374
 375
 376 // Data 900 GeV 10c.pass3 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 377 Double_t nbinsData900[100] = {0.000, 0.488, 0.486, 0.489, 0.495, 0.504, 0.513, 0.521, 0.528, 0.535, 0.541, 0.545, 0.549, 0.554, 0.558, 0.561, 0.565, 0.568, 0.570, 0.573, 0.576, 0.577, 0.580, 0.583, 0.582, 0.587, 0.587, 0.589, 0.591, 0.594, 0.601, 0.599, 0.598, 0.601, 0.608, 0.609, 0.600, 0.606, 0.611, 0.605, 0.612, 0.603, 0.576, 0.606, 0.617, 0.609, 0.572, 0.607, 0.546, 0.000, 0.604, 0.000, 0.689, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 378
 379 // Data 2.76 TeV 11a.pass2 w/oSDD (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 380 Double_t nbinsData276[100] = {0.000, 0.492, 0.491, 0.496, 0.503, 0.513, 0.523, 0.533, 0.541, 0.549, 0.555, 0.560, 0.565, 0.570, 0.574, 0.578, 0.582, 0.585, 0.589, 0.591, 0.594, 0.597, 0.599, 0.601, 0.604, 0.606, 0.608, 0.610, 0.612, 0.613, 0.615, 0.618, 0.619, 0.621, 0.622, 0.624, 0.625, 0.626, 0.626, 0.628, 0.632, 0.633, 0.633, 0.635, 0.635, 0.638, 0.634, 0.641, 0.642, 0.644, 0.638, 0.646, 0.644, 0.647, 0.644, 0.646, 0.640, 0.650, 0.650, 0.655, 0.666, 0.650, 0.659, 0.644, 0.646, 0.657, 0.629, 0.626, 0.602, 0.652, 0.629, 0.000, 0.659, 0.580, 0.000, 0.000, 0.000, 0.000, 0.000, 0.706, 0.000, 0.596, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000};
 381
 382 // Data 7 TeV 10d.pass2 + 10e.pass2 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 383 Double_t nbinsData7[100] = {0.000, 0.492, 0.492, 0.497, 0.505, 0.516, 0.527, 0.538, 0.547, 0.555, 0.562, 0.568, 0.573, 0.578, 0.583, 0.587, 0.591, 0.595, 0.598, 0.601, 0.604, 0.607, 0.610, 0.612, 0.615, 0.617, 0.619, 0.621, 0.623, 0.625, 0.627, 0.629, 0.631, 0.632, 0.634, 0.636, 0.637, 0.639, 0.640, 0.642, 0.643, 0.645, 0.646, 0.647, 0.648, 0.650, 0.650, 0.652, 0.653, 0.654, 0.655, 0.657, 0.657, 0.659, 0.659, 0.659, 0.660, 0.663, 0.663, 0.662, 0.663, 0.665, 0.668, 0.665, 0.667, 0.668, 0.672, 0.671, 0.669, 0.672, 0.672, 0.671, 0.675, 0.666, 0.683, 0.670, 0.673, 0.676, 0.684, 0.672, 0.673, 0.675, 0.689, 0.678, 0.678, 0.720, 0.675, 0.690, 0.664, 0.622, 0.000, 0.698, 0.678, 0.686, 0.751, 0.657, 0.000, 0.663, 0.000, 0.677};
 384
 385 // -- End mean pT values for the multiplicity bins --
 386
 387
 388 // -- Selection of MC/Data and energy; whole sample values --
 389
 390 if (fMC) { // - MC -
 391
 392   if (energy < 1000.) {
 393
 394 //     Printf(" -- MC, 900 GeV -- ");
 395
 396     evMptSample = 0.5281; // MC 900 GeV 10e13 (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 397
 398     nbins = nbinsMC900;
 399   }
 400   else if (energy > 1000. && energy < 3000.) {
 401
 402 //     Printf(" -- MC, 2.76 TeV -- ");
 403
 404     evMptSample = 0.5574; // MC 2.76 TeV 11b10a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 405
 406     nbins = nbinsMC276;
 407   }
 408   else {
 409
 410 //     Printf(" -- MC, 7 TeV -- ");
 411
 412     evMptSample = 0.5809; // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 413 //     evMptSample = 0.5732; // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2) ---> value for MC truth <---
 414
 415 //     evMptSample = 0.5481; // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 416 //     evMptSample = 0.5409; // MC 7 TeV 10f6 (Phojet) (Pt-Range: 0.15 - 2) ---> value for MC truth <---
 417
 418     nbins = nbinsMC7;
 419   }
 420
 421 } // - End MC -
 422 else { // - Data -
 423
 424   if (energy < 1000.) {
 425
 426 //     Printf(" -- Data, 900 GeV -- ");
 427
 428     evMptSample = 0.5263; // Data 900 GeV 10c.pass3 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 429
 430     nbins = nbinsData900;
 431   }
 432   else if (energy > 1000. && energy < 3000.) {
 433
 434 //     Printf(" -- Data, 2.76 TeV -- ");
 435
 436     evMptSample = 0.5548; // Data 2.76 TeV 11a.pass2 w/oSDD (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 437
 438     nbins = nbinsData276;
 439   }
 440   else {
 441
 442 //     Printf(" -- Data, 7 TeV -- ");
 443
 444     evMptSample = 0.5738; // Data 7 TeV 10d.pass2 + 10e.pass2 (Pt-Range: 0.15 - 2) ---> TPC standalone <---
 445
 446     nbins = nbinsData7;
 447   }
 448
 449 } // - End data -
 450
 451 // -- End selection of MC/Data and energy; whole sample values --
 452
 453 // --- End mean pT values ---
 454
 455
 456 // --- Ratio of pT distributions MC ESD / MC truth ---
 457
 458 // MC 7 TeV 10f6a (Pythia6 - Perugia0) (Pt-Range: 0.15 - 2, array-range: 0 - 2.1 GeV/c)
 459 Double_t nbinsPtRatio[84] = {0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.823, 0.908, 0.908, 0.901, 0.903, 0.910, 0.918, 0.931, 0.941, 0.947, 0.951, 0.953, 0.956, 0.959, 0.959, 0.960, 0.961, 0.962, 0.963, 0.963, 0.963, 0.964, 0.964, 0.964, 0.963, 0.963, 0.963, 0.964, 0.964, 0.962, 0.963, 0.964, 0.963, 0.963, 0.962, 0.962, 0.961, 0.963, 0.964, 0.964, 0.965, 0.965, 0.965, 0.964, 0.964, 0.964, 0.964, 0.964, 0.963, 0.964, 0.961, 0.957, 0.958, 0.959, 0.959, 0.958, 0.958, 0.957, 0.959, 0.956, 0.956, 0.954, 0.954, 0.956, 0.955, 0.954, 0.954, 0.954, 0.955, 0.954, 0.951, 0.952, 0.951, 0.949, 0.000, 0.000, 0.000, 0.000};
 460
 461 // --- End ratio of pT distributions MC ESD / MC truth ---
 462
 463
 464   if (!fMC || fMCType < 3) {
 465     nrESDTracks = fESD->GetNumberOfTracks();
 466   }
 467
 468   if (fMC) {
 469     nrMCTracks = stack->GetNtrack();
 470   }
 471
 472 //   Printf("\n\n nrESDTracks: %i   nrMCTracks : %i \n",nrESDTracks,nrMCTracks);
 473
 474
 475 if (fMC) { // - MC -
 476
 477  if (fMCType == 0) { // MC, Type = ESD
 478
 479
 480   // --- Loop over all tracks of one ESD event ---
 481   for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
 482
 483     AliESDtrack* track = fESD->GetTrack(iTracks);
 484
 485     if (!track) {
 486       Printf("ERROR: Could not receive track %d\n", iTracks);
 487       continue;
 488     }
 489
 490       if(!fESDTrackCuts->AcceptTrack(track))continue;
 491
 492         trackPt = track->Pt();
 493         fPtSpec->Fill(trackPt);
 494         tracks[nrTracks] = trackPt;
 495
 496         trackEta = track->Eta();
 497         fEta->Fill(trackEta);
 498
 499         trackPhi = track->Phi();
 500
 501         if (trackEta > 0.) {
 502           fEtaPhiPlus->Fill(trackPhi);
 503         }
 504         else if (trackEta < 0.) {
 505           fEtaPhiMinus->Fill(trackPhi);
 506         }
 507
 508         sumPt = sumPt + trackPt;
 509         nrTracks++;
 510
 511 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
 512
 513   } // --- End track loop ---
 514
 515
 516  } // End MC, Type = ESD
 517  else if (fMCType > 0) { // MC, Type = MC truth
 518
 519
 520   // --- Loop over all tracks of one MC event - MC truth / stack ---
 521   for (Int_t iTracks = 0; iTracks < nrMCTracks; iTracks++) {
 522
 523     TParticle *track = stack->Particle(iTracks);
 524
 525     if (!track) {
 526       Printf("ERROR: Could not receive track %d\n", iTracks);
 527       continue;
 528     }
 529
 530         // only charged particles
 531         TParticlePDG* pdg = track->GetPDG();
 532         if(!pdg) continue;
 533
 534         Double_t charge = pdg->Charge()/3.;
 535         if ( TMath::Abs(charge) < 0.001 ) continue;
 536
 537         // only physical primary
 538         Bool_t prim = stack->IsPhysicalPrimary(iTracks);
 539         if(!prim) continue;
 540
 541         // get pt, eta and phi and cut in eta and pt
 542         trackPt = track->Pt();
 543         trackEta = track->Eta();
 544         trackPhi = track->Phi();
 545
 546 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f  ",trackPt,trackEta,trackPhi);
 547
 548         if (trackEta > 0.8) continue;
 549         if (trackEta < -0.8) continue;
 550
 551         if (trackPt < 0.15) continue;
 552         if (trackPt > 2.) continue;
 553
 554 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f  ",trackPt,trackEta,trackPhi);
 555 //      Printf("-> Track accepted!");
 556
 557
 558         // -- Reject tracks according to ratio of pT distribution MC ESD / MC truth --
 559         if (fMCType == 2 || fMCType == 4) { // MC, Type = mod. MC truth
 560
 561           trackPt40 = trackPt*40.;
 562           trackPtBin = floor(trackPt40);
 563
 564           ptRatio = nbinsPtRatio[trackPtBin];
 565
 566 //        Printf("trackPt: %f   trackPt40: %f   trackPtBin: %i   ptRatio: %.3f   ",trackPt,trackPt40,trackPtBin,ptRatio);
 567
 568           random = fRandom3->Rndm();
 569 //        Printf("Random: %f  \n",random);
 570
 571           if (random > ptRatio) continue;
 572 //        Printf("accepted random: %f \n",random);
 573
 574         } // -- End reject tracks according to ratio of pT distribution MC ESD / MC truth --
 575
 576
 577         fPtSpec->Fill(trackPt);
 578         tracks[nrTracks] = trackPt;
 579
 580         fEta->Fill(trackEta);
 581
 582         if (trackEta > 0.) {
 583           fEtaPhiPlus->Fill(trackPhi);
 584         }
 585         else if (trackEta < 0.) {
 586           fEtaPhiMinus->Fill(trackPhi);
 587         }
 588
 589         sumPt = sumPt + trackPt;
 590         nrTracks++;
 591
 592 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
 593
 594   } // --- End stack track loop ---
 595
 596 //   Printf("nrTracks MC truth: %i \n",nrTracks);
 597
 598
 599   // MC, Type = MC truth with corresponding ESD
 600   if (fMCType < 3) {
 601
 602     // --- Loop over all tracks of one MC event - MC ESD ---
 603     for (Int_t jTracks = 0; jTracks < nrESDTracks; jTracks++) {
 604       AliESDtrack* track = fESD->GetTrack(jTracks);
 605       if (!track) {
 606         Printf("ERROR: Could not receive track %d\n", jTracks);
 607         continue;
 608       }
 609
 610         if(!fESDTrackCuts->AcceptTrack(track))continue;
 611
 612           trackPt2 = track->Pt();
 613           fPtSpec2->Fill(trackPt2);
 614
 615 //        nrTracks2++;
 616
 617 //        Printf("Track Pt ESD = %.3f\n\n",trackPt2);
 618
 619     } // --- End ESD track loop ---
 620
 621 //     Printf("nrTracks ESD: %i \n\n",nrTracks2);
 622
 623   } // End MC, Tpye = MC truth with corresponding ESD
 624
 625  } // End MC, Type = MC truth
 626
 627
 628 } // - End MC -
 629 else { // - Data -
 630
 631
 632   // --- Loop over all tracks of one ESD event ---
 633   for (Int_t iTracks = 0; iTracks < nrESDTracks; iTracks++) {
 634
 635     AliESDtrack* track = fESD->GetTrack(iTracks);
 636
 637     if (!track) {
 638       Printf("ERROR: Could not receive track %d\n", iTracks);
 639       continue;
 640     }
 641
 642       if(!fESDTrackCuts->AcceptTrack(track))continue;
 643
 644         trackPt = track->Pt();
 645         fPtSpec->Fill(trackPt);
 646         tracks[nrTracks] = trackPt;
 647
 648         trackEta = track->Eta();
 649         fEta->Fill(trackEta);
 650
 651         trackPhi = track->Phi();
 652
 653         if (trackEta > 0.) {
 654           fEtaPhiPlus->Fill(trackPhi);
 655         }
 656         else if (trackEta < 0.) {
 657           fEtaPhiMinus->Fill(trackPhi);
 658         }
 659
 660         sumPt = sumPt + trackPt;
 661         nrTracks++;
 662
 663 //      Printf("Track Pt = %.3f   Track Eta = %.3f   Track Phi = %3f\n",trackPt,trackEta,trackPhi);
 664
 665   } // --- End track loop ---
 666
 667 } // - End data -
 668
 669
 670   // --- Calculation of various values and filling of histograms (for remaining events with N_acc > 0) ---
 671
 672   if(nrTracks != 0) {
 673
 674     if (!fMC || fMCType < 3) {
 675
 676       // VtxZ after all track cuts
 677       fVtxZTrackCuts->Fill(vtxZ);
 678
 679       // Differences between vertices after all cuts
 680       fVtxZDiff1b->Fill(vtxZdiff1);
 681       fVtxZDiff2b->Fill(vtxZdiff2);
 682       fVtxZDiff3b->Fill(vtxZdiff3);
 683
 684     }
 685
 686     // Multiplicity distribution
 687     fMult->Fill(nrTracks);
 688
 689 //     // Number of pairs in event
 690 //     nrPairs = 0.5 * nrTracks * (nrTracks-1);
 691
 692     // Calculation of mean Pt and mean Pt Squared
 693     eventMeanPt = sumPt / nrTracks;
 694     eventMeanPtSq = eventMeanPt * eventMeanPt;
 695
 696     // Mean-Pt and Mean-Pt squared
 697     fEventMeanPt->Fill(eventMeanPt);
 698     fEventMeanPtSq->Fill(eventMeanPtSq);
 699     fEventMeanPtMult->Fill(nrTracks,eventMeanPt);
 700
 701     // Mean-Pt and Mean-Pt squared depending on multiplicity
 702     fMultEventMeanPt->Fill(nrTracks,eventMeanPt);
 703     fMultEventMeanPtSq->Fill(nrTracks,eventMeanPtSq);
 704
 705 //     Printf("nrTracks: %i   sumPt: %.8f   meanPt: %.8f   meanPtSq: %.8f\n",nrTracks,sumPt,eventMeanPt,eventMeanPtSq);
 706
 707
 708     // --- Two-particle correlator ---
 709
 710     evMptMult = nbins[nrTracks];
 711 //     Printf("nrTracks = %3i   evMptMult = %.3f   evMptSample = %.3f \n\n",nrTracks,evMptMult,evMptSample);
 712
 713     for (int i=0; i<nrTracks; i++) {
 714
 715       for (int j=i+1; j<nrTracks; j++) {
 716
 717         twoPartCorrPair = (tracks[i] - evMptMult) * (tracks[j] - evMptMult); // Multiplicity bins
 718         twoPartCorrEv = twoPartCorrEv + twoPartCorrPair;
 719
 720         twoPartCorrPairSample = (tracks[i] - evMptSample) * (tracks[j] - evMptSample); // Whole sample
 721         twoPartCorrEvSample = twoPartCorrEvSample + twoPartCorrPairSample;
 722
 723       }
 724     }
 725
 726     // Multiplicity bins
 727     twoPartCorrEvSq = twoPartCorrEv * twoPartCorrEv;
 728     fTwoPartCorrEv->Fill(nrTracks,twoPartCorrEv);
 729     fTwoPartCorrEvSq->Fill(nrTracks,twoPartCorrEvSq);
 730
 731     // Whole sample
 732     twoPartCorrEvSampleSq = twoPartCorrEvSample * twoPartCorrEvSample;
 733     fTwoPartCorrEvSample->Fill(nrTracks,twoPartCorrEvSample);
 734     fTwoPartCorrEvSampleSq->Fill(nrTracks,twoPartCorrEvSampleSq);
 735
 736     // --- End two-particle correlator ---
 737
 738
 739   } // --- End calculation of various values and filling of histograms ---
 740
 741
 742   // Post output data
 743   PostData(1, fOutputList);
 744
 745   // Delet pointer
 746   if (fRandom3) delete fRandom3; fRandom3 = 0;
 747
 748 }
 749
 750 void AliAnalysisTaskPtFluc::Terminate(Option_t *)
 751 {
 752   // Called once at the end of the query
 753
 754   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 755   if (!fOutputList) {
 756     Printf("ERROR: fOutputList not available\n");
 757     return;
 758   }
 759
 760 }