PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskAntiHe4.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TF1.h"
   4 #include "TH1F.h"
   5 #include "TH2F.h"
   6 #include "TH3F.h"
   7 #include "THn.h"
   8 #include "TCanvas.h"
   9 #include "AliESDtrackCuts.h"
  10 #include "AliESDpid.h"
  11 #include "AliESDVertex.h"
  12 #include "TFile.h"
  13
  14 #include "AliAnalysisTask.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliESDEvent.h"
  18 #include "AliESDInputHandler.h"
  19
  20 #include "AliMCEventHandler.h"
  21 #include "AliMCEvent.h"
  22 #include "AliStack.h"
  23
  24 #include "AliAnalysisTaskAntiHe4.h"
  25
  26
  27 ///////////////////////////////////////////////////////////
  28 //                                                       //
  29 // Analysis for the observation of anti-alpha particles. //
  30 //                                                       //
  31 ///////////////////////////////////////////////////////////
  32
  33 using std::cout;
  34 using std::endl;
  35
  36 ClassImp(AliAnalysisTaskAntiHe4)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4()
  40   : AliAnalysisTaskSE(),
  41   fEventHandler(0),
  42   fESD(0),
  43   fHistCentralityClass10(0),
  44   fHistCentralityPercentile(0),
  45   fHistTriggerStat(0),
  46   fHistTriggerStatAfterEventSelection(0),
  47   fHistDEDx(0),
  48   fHistTOF3D(0),
  49   fHistAlpha(0),
  50   fHistAlphaSignal(0),
  51   fGraphAlphaSignal(0),
  52   fNCounter(0),
  53   fHistDeDx(0),
  54   fHistDeDxRegion(0),
  55   fHistDeDxSharp(0),
  56   fHistTOF2D(0),
  57   fHistTOFnuclei(0x0),
  58   fNTriggers(5),
  59   fBBParametersLightParticles(),
  60   fBBParametersNuclei(),
  61   fMCtrue(0),
  62   fTriggerFired(),
  63   fESDtrackCuts(0),
  64   fESDtrackCutsSharp(0),
  65   fESDpid(0),
  66   fAntiAlpha(0),
  67   fHistHelium4PtGen(0),
  68   fHistHelium4PtGenPrim(0),
  69   fHistHelium4PtGenSec(0),
  70   fHistHelium4PtGenEta(0),
  71   fHistHelium4PtGenPrimEta(0),
  72   fHistAntiHelium4PtGen(0),
  73   fHistAntiHelium4PtGenPrim(0),
  74   fHistAntiHelium4PtGenSec(0),
  75   fHistAntiHelium4PtGenEta(0),
  76   fHistHelium4PtAso(0),
  77   fHistHelium4PtAsoPrim(0),
  78   fHistHelium4PtAsoSec(0),
  79   fHistAntiHelium4PtAso(0),
  80   fTree(0),
  81   fOutputContainer(0)
  82 {
  83   // default Constructor
  84
  85
  86   // Define input and output slots here
  87 }
  88
  89 //________________________________________________________________________
  90 AliAnalysisTaskAntiHe4::AliAnalysisTaskAntiHe4(const char *name)
  91   : AliAnalysisTaskSE(name),
  92     fEventHandler(0),
  93     fESD(0),
  94     fHistCentralityClass10(0),
  95     fHistCentralityPercentile(0),
  96     fHistTriggerStat(0),
  97     fHistTriggerStatAfterEventSelection(0),
  98     fHistDEDx(0),
  99     fHistTOF3D(0),
 100     fHistAlpha(0),
 101     fHistAlphaSignal(0),
 102     fGraphAlphaSignal(0),
 103     fNCounter(0),
 104     fHistDeDx(0),
 105     fHistDeDxRegion(0),
 106     fHistDeDxSharp(0),
 107     fHistTOF2D(0),
 108     fHistTOFnuclei(0x0),
 109     fNTriggers(5),
 110     fBBParametersLightParticles(),
 111     fBBParametersNuclei(),
 112     fMCtrue(0),
 113     fTriggerFired(),
 114     fESDtrackCuts(0),
 115     fESDtrackCutsSharp(0),
 116     fESDpid(0),
 117     fAntiAlpha(0),
 118     fHistHelium4PtGen(0),
 119     fHistHelium4PtGenPrim(0),
 120     fHistHelium4PtGenSec(0),
 121     fHistHelium4PtGenEta(0),
 122     fHistHelium4PtGenPrimEta(0),
 123     fHistAntiHelium4PtGen(0),
 124     fHistAntiHelium4PtGenPrim(0),
 125     fHistAntiHelium4PtGenSec(0),
 126     fHistAntiHelium4PtGenEta(0),
 127     fHistHelium4PtAso(0),
 128     fHistHelium4PtAsoPrim(0),
 129     fHistHelium4PtAsoSec(0),
 130     fHistAntiHelium4PtAso(0),
 131     fTree(0),
 132     fOutputContainer(0)
 133 {
 134   // Constructor
 135
 136   // Define input and output slots here
 137   // Input slot #0 works with a TChain
 138   DefineInput(0, TChain::Class());
 139   // Output slot #0 writes into a TH1 container
 140   DefineOutput(1, TObjArray::Class());
 141   DefineOutput(2, TTree::Class());
 142   //
 143   // cuts for candidates
 144   //
 145   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 146   //
 147   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 148   fESDtrackCuts->SetMinNClustersTPC(70);
 149   fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
 150   fESDtrackCuts->SetMaxDCAToVertexXY(3);
 151   fESDtrackCuts->SetMaxDCAToVertexZ(2);
 152   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 153   //fESDtrackCuts->SetRequireITSRefit(kTRUE);
 154   fESDtrackCuts->SetMinNClustersITS(2);
 155   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 156   //
 157   // cuts for final plots
 158   //
 159   fESDtrackCutsSharp  = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 160   fESDtrackCutsSharp->SetAcceptKinkDaughters(kFALSE);
 161   fESDtrackCutsSharp->SetMinNClustersTPC(80);
 162   fESDtrackCutsSharp->SetMaxChi2PerClusterITS(10);// TO BE INVESTIGATED !!!!!!!!!!!!!!
 163   fESDtrackCutsSharp->SetMaxChi2PerClusterTPC(5);
 164   fESDtrackCutsSharp->SetRequireTPCRefit(kTRUE);
 165   fESDtrackCutsSharp->SetRequireITSRefit(kTRUE);
 166   fESDtrackCutsSharp->SetMinNClustersITS(2);
 167   fESDtrackCutsSharp->SetMaxDCAToVertexXY(0.1);
 168   fESDtrackCutsSharp->SetMaxDCAToVertexZ(0.5);
 169   fESDtrackCutsSharp->SetEtaRange(-0.8,0.8);
 170
 171   //ESD Track cuts  from TestFilterRawTask
 172   /*  fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 173   fESDtrackCuts->SetMinNClustersTPC(80);
 174   fESDtrackCuts->SetMinNClustersITS(2);
 175   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 176   fESDtrackCuts->SetMaxChi2PerClusterTPC(5);
 177   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 178   fESDtrackCuts->SetRequireITSRefit(kTRUE);
 179   fESDtrackCuts->SetEtaRange(-1,1);
 180   fESDtrackCuts->SetMaxDCAToVertexXY(1);
 181   fESDtrackCuts->SetMaxDCAToVertexZ(2);
 182   //test strickter cuts
 183   //fESDtrackCuts->SetMaxDCAToVertexXY(0.1);
 184   //fESDtrackCuts->SetMaxDCAToVertexZ(0.5);
 185   //fESDtrackCuts->SetEtaRange(-0.8,0.8);
 186   */
 187
 188   fMCtrue = kTRUE;
 189
 190 }
 191
 192 //________________________________________________________________________
 193 void AliAnalysisTaskAntiHe4::UserCreateOutputObjects()
 194 {
 195   // Create histograms
 196   // Called once
 197   //
 198   //histogram to count number of events
 199   //
 200   fHistCentralityClass10  = new TH1F("fHistCentralityClass10", "centrality with class10", 11, -0.5, 10.5);
 201   fHistCentralityClass10->GetXaxis()->SetTitle("Centrality");
 202   fHistCentralityClass10->GetYaxis()->SetTitle("Entries");
 203   //
 204   fHistCentralityPercentile  = new TH1F("fHistCentralityPercentile", "centrality with percentile", 101, -0.1, 100.1);
 205   fHistCentralityPercentile->GetXaxis()->SetTitle("Centrality");
 206   fHistCentralityPercentile->GetYaxis()->SetTitle("Entries");
 207   //
 208   //trigger statitics histogram
 209   //
 210   fHistTriggerStat = new TH1F("fHistTriggerStat","Trigger statistics", fNTriggers,-0.5,fNTriggers-0.5);
 211   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 212   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 213     fHistTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 214
 215   fHistTriggerStatAfterEventSelection = new TH1F("fHistTriggerStatAfterEventSelection","Trigger statistics after event selection", fNTriggers,-0.5,fNTriggers-0.5);
 216   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 217     fHistTriggerStatAfterEventSelection->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 218
 219   //dE/dx performance
 220   fHistDEDx= new TH3F("fHistDEDx","DEDx",1000,0.01,100,1000,1,2000,2,-2,2);
 221   BinLogAxis(fHistDEDx, 0);
 222   BinLogAxis(fHistDEDx, 1);
 223   fHistDEDx->GetXaxis()->SetTitle("p_{tot}/sign");
 224   fHistDEDx->GetYaxis()->SetTitle("TPC signal");
 225
 226   fHistDeDx = new TH2F("fHistDeDx", "dE/dx", 1000, 0.1, 6.0, 1000, 0.0, 1000);
 227   fHistDeDx->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
 228   fHistDeDx->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
 229   BinLogAxis(fHistDeDx);
 230
 231   fHistDeDxRegion = new TH3F("fHistDeDxRegion", "dE/dx", 400, 0., 6.0, 300, 0., 3., 4, -0.5, 3.5);
 232   fHistDeDxRegion->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
 233   fHistDeDxRegion->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
 234
 235   fHistDeDxSharp = new TH2F("fHistDeDxSharp", "dE/dx", 1000, 0.1, 6.0, 1000, 0.0, 1000);
 236   fHistDeDxSharp->GetYaxis()->SetTitle("TPC dE/dx signal (a.u.)");
 237   fHistDeDxSharp->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
 238   BinLogAxis(fHistDeDxSharp);
 239
 240   //TOF performance
 241   fHistTOF3D = new TH3F("fHistTOF3D","mass determination from TOF",400,0.25,4.5,2,-0.5,1.5,2,-1,1);
 242
 243   fHistTOF2D = new TH2F("fHistTOF2D", "TOF2D", 500, 0.0, 10., 2250, 0.2, 1.1);
 244   fHistTOF2D->GetYaxis()->SetTitle("#beta");
 245   fHistTOF2D->GetXaxis()->SetTitle("p (GeV/c)");
 246
 247   fHistTOFnuclei = new TH2F("fHistTOFnuclei", "TOF", 500, 0.0, 10., 2250, 0.7, 1.);
 248   fHistTOFnuclei->GetYaxis()->SetTitle("#beta");
 249   fHistTOFnuclei->GetXaxis()->SetTitle("#frac{p}{z} (GeV/c)");
 250
 251   //alphas
 252   fHistAlpha = new TH1F("fHistAlpha", "Anti-Alpha", 22, 1.12, 4.31);
 253   fHistAlpha->GetYaxis()->SetTitle("Counts");
 254   fHistAlpha->GetXaxis()->SetTitle("#frac{m^{2}}{z^{2}} (GeV^{2}/c^{4})");
 255
 256   fHistAlphaSignal  = new TH1F("fHistAlphaSignal", "Anti-Alpha", 22, 1.12, 4.31);
 257   fHistAlphaSignal->GetYaxis()->SetTitle("Counts");
 258   fHistAlphaSignal->GetXaxis()->SetTitle("#frac{m^{2}}{z^{2}} (GeV^{2}/c^{4})");
 259
 260   fGraphAlphaSignal = new TGraph(20);
 261   fNCounter = 0;
 262   //
 263   //  (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
 264   //
 265   TString axisNameMult[4]={"dca","sign","particleType","ptot"};
 266   TString axisTitleMult[4]={"dca","sign","particleType","ptot"};
 267   const Int_t kDcaBins = 76;
 268   Double_t binsDca[kDcaBins+1] = {-3,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3};
 269   //
 270   //                     (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
 271   Int_t    binsAntiAlpha[4] = {77,         2,                  4,        100};
 272   Double_t xminAntiAlpha[4] = {-3,        -2,               -0.5,          0};
 273   Double_t xmaxAntiAlpha[4] = { 3,         2,                3.5,          6};
 274   //
 275   fAntiAlpha = new THnF("fAntiAlpha", "AntiAlpha; (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot", 4, binsAntiAlpha, xminAntiAlpha, xmaxAntiAlpha);
 276   //
 277   fAntiAlpha->GetAxis(0)->Set(kDcaBins, binsDca);
 278   for (Int_t iaxis=0; iaxis<4;iaxis++){
 279     fAntiAlpha->GetAxis(iaxis)->SetName(axisNameMult[iaxis]);
 280     fAntiAlpha->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]);
 281   }
 282   //
 283   //
 284   //Create histograms for MC
 285   //generated histogramms
 286   fHistHelium4PtGen = new TH1F("fHistHelium4PtGen", "Generated  ^{4}He", 200, 0.0, 10.0);
 287   fHistHelium4PtGen->GetYaxis()->SetTitle("Counts");
 288   fHistHelium4PtGen->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 289
 290   fHistHelium4PtGenPrim = new TH1F("fHistHelium4PtGenPrim", "Primary generated  ^{4}He", 200, 0.0, 10.0);
 291   fHistHelium4PtGenPrim->GetYaxis()->SetTitle("Counts");
 292   fHistHelium4PtGenPrim->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 293
 294   fHistHelium4PtGenSec = new TH1F("fHistHelium4PtGenSec", "Secondary generated  ^{4}He", 200, 0.0, 10.0);
 295   fHistHelium4PtGenSec->GetYaxis()->SetTitle("Counts");
 296   fHistHelium4PtGenSec->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 297
 298   fHistHelium4PtGenEta = new TH1F("fHistHelium4PtGenEta", "Generated  ^{4}He in #eta < 0.8", 200, 0.0, 10.0);
 299   fHistHelium4PtGenEta->GetYaxis()->SetTitle("Counts");
 300   fHistHelium4PtGenEta->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 301
 302   fHistHelium4PtGenPrimEta = new TH1F("fHistHelium4PtGenPrimEta", "Primary generated ^{4}He in #eta < 0.8", 200, 0.0, 10.0);
 303   fHistHelium4PtGenPrimEta->GetYaxis()->SetTitle("Counts");
 304   fHistHelium4PtGenPrimEta->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 305
 306   fHistAntiHelium4PtGen = new TH1F("fHistAntiHelium4PtGen", "Generated  #bar{^{4}He}", 200, 0.0, 10.0);
 307   fHistAntiHelium4PtGen->GetYaxis()->SetTitle("Counts");
 308   fHistAntiHelium4PtGen->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 309
 310   fHistAntiHelium4PtGenPrim = new TH1F("fHistAntiHelium4PtGenPrim", "Primary generated  #bar{^{4}He}", 200, 0.0, 10.0);
 311   fHistAntiHelium4PtGenPrim->GetYaxis()->SetTitle("Counts");
 312   fHistAntiHelium4PtGenPrim->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 313
 314   fHistAntiHelium4PtGenSec = new TH1F("fHistAntiHelium4PtGenSec", "Secondary generated  #bar{^{4}He}", 200, 0.0, 10.0);
 315   fHistAntiHelium4PtGenSec->GetYaxis()->SetTitle("Counts");
 316   fHistAntiHelium4PtGenSec->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 317
 318   fHistAntiHelium4PtGenEta = new TH1F("fHistAntiHelium4PtGenEta", "Generated  #bar{^{4}He} in #eta < 0.8", 200, 0.0, 10.0);
 319   fHistAntiHelium4PtGenEta->GetYaxis()->SetTitle("Counts");
 320   fHistAntiHelium4PtGenEta->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 321
 322   //associated histograms
 323   fHistHelium4PtAso = new TH1F("fHistHelium4PtAso", "Associated  ^{4}He", 200, 0.0, 10.0);
 324   fHistHelium4PtAso->GetYaxis()->SetTitle("Counts");
 325   fHistHelium4PtAso->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 326
 327   fHistHelium4PtAsoPrim = new TH1F("fHistHelium4PtAsoPrim", "Associated prim ^{4}He", 200, 0.0, 10.0);
 328   fHistHelium4PtAsoPrim->GetYaxis()->SetTitle("Counts");
 329   fHistHelium4PtAsoPrim->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 330
 331   fHistHelium4PtAsoSec = new TH1F("fHistHelium4PtAsoSec", "Associated sec  ^{4}He", 200, 0.0, 10.0);
 332   fHistHelium4PtAsoSec->GetYaxis()->SetTitle("Counts");
 333   fHistHelium4PtAsoSec->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 334
 335   fHistAntiHelium4PtAso = new TH1F("fHistAntiHelium4PtAso", "Associated  #bar{^{4}He}", 200, 0.0, 10.0);
 336   fHistAntiHelium4PtAso->GetYaxis()->SetTitle("Counts");
 337   fHistAntiHelium4PtAso->GetXaxis()->SetTitle("p_{T}  (GeV/c)");
 338
 339   //
 340   //
 341   fOutputContainer = new TObjArray(1);
 342   fOutputContainer->SetOwner(kTRUE);
 343   fOutputContainer->SetName(GetName());
 344   fOutputContainer->Add(fHistCentralityClass10);
 345   fOutputContainer->Add(fHistCentralityPercentile);
 346   fOutputContainer->Add(fHistTriggerStat);
 347   fOutputContainer->Add(fHistTriggerStatAfterEventSelection);
 348   fOutputContainer->Add(fHistDEDx);
 349   fOutputContainer->Add(fHistDeDx);
 350   fOutputContainer->Add(fHistDeDxRegion);
 351   fOutputContainer->Add(fHistDeDxSharp);
 352   fOutputContainer->Add(fAntiAlpha);
 353   fOutputContainer->Add(fHistAlpha);
 354   fOutputContainer->Add(fHistAlphaSignal);
 355   fOutputContainer->Add(fGraphAlphaSignal);
 356   fOutputContainer->Add(fHistTOF3D);
 357   fOutputContainer->Add(fHistTOF2D);
 358   fOutputContainer->Add(fHistTOFnuclei);
 359   fOutputContainer->Add(fHistHelium4PtGen);
 360   fOutputContainer->Add(fHistHelium4PtGenPrim);
 361   fOutputContainer->Add(fHistHelium4PtGenSec);
 362   fOutputContainer->Add(fHistHelium4PtGenEta);
 363   fOutputContainer->Add(fHistHelium4PtGenPrimEta);
 364   fOutputContainer->Add(fHistAntiHelium4PtGen);
 365   fOutputContainer->Add(fHistAntiHelium4PtGenPrim);
 366   fOutputContainer->Add(fHistAntiHelium4PtGenSec);
 367   fOutputContainer->Add(fHistAntiHelium4PtGenEta);
 368   fOutputContainer->Add(fHistHelium4PtAso);
 369   fOutputContainer->Add(fHistHelium4PtAsoPrim);
 370   fOutputContainer->Add(fHistHelium4PtAsoSec);
 371   fOutputContainer->Add(fHistAntiHelium4PtAso);
 372
 373
 374   //------------ Tree and branch definitions ----------------//
 375   fTree = new TTree("tree", " alpha tree");
 376
 377   //------------ Event variables ------------//
 378   fTree->Branch("Name",Name,"Name/C");
 379   fTree->Branch("Evnt",&evnt, "evnt/I");
 380   fTree->Branch("itrk", &itrk, "itrk/I");
 381
 382   //-------------------------------------------//
 383   //----------- Track variables --------------//
 384
 385   fTree->Branch("TrkPtot",TrkPtot,"TrkPtot[itrk]/F");
 386   fTree->Branch("TPCPtot",TPCPtot,"TPCPtot[itrk]/F");
 387   fTree->Branch("DeDx",DeDx,"DeDx[itrk]/F");
 388   fTree->Branch("Sign",Sign,"Sign[itrk]/F");
 389   fTree->Branch("DCAXY",DCAXY,"DCAXY[itrk]/F");
 390   fTree->Branch("DCAZ",DCAZ,"DCAZ[itrk]/F");
 391   fTree->Branch("ITSnCluster",ITSnCluster,"ITSnCluster[itrk]/F");
 392   fTree->Branch("TPCNsignal",TPCNsignal,"TPCNsignal[itrk]/F");
 393   fTree->Branch("Mass",Mass,"Mass[itrk]/F");
 394   //
 395   fTree->Branch("ITSRefit",ITSRefit,"ITSRefit[itrk]/F");
 396   fTree->Branch("TOFtime",TOFtime,"TOFtime[itrk]/F");
 397   fTree->Branch("TOFRefit",TOFRefit,"TOFRefit[itrk]/F");
 398   fTree->Branch("TOFout",TOFout,"TOFout[itrk]/F");
 399   //
 400   fTree->Branch("ITSsignal",ITSsignal,"ITSsignal[itrk]/F");
 401   fTree->Branch("SharedClusters",SharedClusters,"SharedClusters[itrk]/F");
 402   fTree->Branch("fFileName",fFileName,"fFileName/C");
 403   fTree->Branch("fEventNumber",fEventNumber,"fEventNumber/I");
 404   //
 405   fTree->Branch("fAssociated",fAssociated,"fAssociated[itrk]/O");
 406   fTree->Branch("fTrackPt",fTrackPt,"fTrackPt[itrk]/F");
 407
 408 }
 409
 410 //________________________________________________________________________
 411 void AliAnalysisTaskAntiHe4::UserExec(Option_t *)
 412 {
 413   // Main loop
 414   // Called for each event
 415
 416   //get Event-Handler for the trigger information
 417   fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 418   if (!fEventHandler) {
 419     AliError("Could not get InputHandler");
 420     //return -1;
 421     return;
 422   }
 423
 424   // Monte Carlo
 425   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 426   if (!eventHandler){
 427     //Printf("ERROR: Could not retrieve MC event handler");
 428     fMCtrue = kFALSE;
 429   }
 430
 431   AliMCEvent* mcEvent = 0x0;
 432   AliStack* stack = 0x0;
 433   if (eventHandler) mcEvent = eventHandler->MCEvent();
 434   if (!mcEvent){
 435     //Printf("ERROR: Could not retrieve MC event");
 436     if (fMCtrue) return;
 437   }
 438
 439   if (fMCtrue){
 440     stack = mcEvent->Stack();
 441     if (!stack) return;
 442   }
 443
 444   //look for the generated particles
 445   MCGenerated(stack);
 446
 447   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 448   if (!fESD) {
 449     //Printf("ERROR: fESD not available");
 450     return;
 451   }
 452
 453   if (SetupEvent() < 0) {
 454     PostData(1, fOutputContainer);
 455     return;
 456   }
 457
 458   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 459   if(vertex->GetNContributors()<1) {
 460     // SPD vertex
 461     vertex = fESD->GetPrimaryVertexSPD();
 462     if(vertex->GetNContributors()<1) {
 463       PostData(1, fOutputContainer);
 464       return;
 465     }
 466
 467   }
 468   // check if event is selected by physics selection class
 469   //
 470   Bool_t isSelected = kFALSE;
 471   isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 472   if (!isSelected || TMath::Abs(vertex->GetZv()) > 10) {
 473     PostData(1, fOutputContainer);
 474     return;
 475   }
 476
 477   //
 478   //centrality selection in PbPb
 479   //
 480   Int_t centralityClass10 = -1;
 481   Int_t centralityPercentile = -1;
 482   //
 483   if (fESD->GetEventSpecie() == 4) { //species == 4 == PbPb
 484     // PbPb
 485     AliCentrality *esdCentrality = fESD->GetCentrality();
 486     centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 487     centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
 488     if(!fMCtrue){
 489       if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
 490     }
 491   }
 492   //
 493   if (!fTriggerFired[0] && !fTriggerFired[1] && !fTriggerFired[2]) return; // select only events which pass kMB, kCentral, kSemiCentral
 494   //
 495   fHistCentralityClass10->Fill(centralityClass10);
 496   fHistCentralityPercentile->Fill(centralityPercentile);
 497   //
 498   if(fTriggerFired[0] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(0);
 499   if(fTriggerFired[1] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(1);
 500   if(fTriggerFired[2] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(2);
 501   if(fTriggerFired[3] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(3);
 502   if(fTriggerFired[4] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(4);
 503   //
 504   //
 505   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 506   if (!fESDpid) {
 507     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
 508     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
 509   }
 510   //
 511   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for th // for Anti-Alpha
 512   evnt =  fESD->GetEventNumberInFile();
 513   sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", Name);
 514   itrk = 0;
 515   //
 516   Int_t runNumber = 0;
 517   runNumber = fESD->GetRunNumber();
 518   //
 519   Bool_t fillTree = kFALSE;
 520   // Track loop to fill the spectram
 521   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 522     //
 523     AliESDtrack* track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTracks));
 524     if (!fESDtrackCuts->AcceptTrack(track)) continue;
 525     //
 526     Double_t nClustersTPCPID=track->GetTPCsignalN();
 527     if(nClustersTPCPID < 60) continue;
 528     //
 529     if (!track->GetInnerParam()) continue;
 530     //
 531     Double32_t signal[4] = {0,0,0,0};
 532     Char_t ncl[3];
 533     Char_t nrows[3];
 534     if (track->GetTPCdEdxInfo()) {
 535       track->GetTPCdEdxInfo()->GetTPCSignalRegionInfo(signal,ncl,nrows);
 536     }
 537     //
 538     Double_t ptot = track->GetInnerParam()->GetP();
 539     Double_t ptotInc = track->GetP(); // total momentum of the incoming particle
 540     Double_t sign = track->GetSign();
 541     //
 542     Double_t eta = TMath::Abs(track->Eta());
 543     TBits shared = track->GetTPCSharedMap();
 544
 545     track->GetImpactParameters(dca, cov);
 546     Float_t dcaXY = TMath::Sqrt(dca[0]*dca[0]);
 547     Float_t dcaXYsign = dca[0];
 548     Float_t dcaZ  = TMath::Sqrt(dca[1]*dca[1]);
 549     //
 550     //
 551     Double_t tpcSignal = track->GetTPCsignal();
 552     //
 553     //PID via specific energy loss in the TPC
 554     //define the arrays for the Bethe-Bloch-Parameters
 555     //
 556
 557     SetBBParameters(runNumber);
 558
 559     //define expected signals for the various species
 560     Double_t expSignalProton = AliExternalTrackParam::BetheBlochAleph(ptot/0.93827,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 561     Double_t expSignalDeuteron = AliExternalTrackParam::BetheBlochAleph(ptot/1.875612,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 562     Double_t expSignalTriton = AliExternalTrackParam::BetheBlochAleph(ptot/2.808921,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 563
 564     Double_t expSignalHelium3 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/2.80941,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
 565     Double_t expSignalHelium4 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/3.728401,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
 566
 567     //
 568     Float_t mass = 0;
 569     Float_t time = -1;
 570     Float_t beta = 0;
 571     Float_t gamma = -1;
 572     Bool_t  hasTOFout  = kFALSE;
 573     Bool_t  hasTOFtime = kFALSE;
 574     Float_t length = track->GetIntegratedLength();
 575     UInt_t  status = track->GetStatus();
 576     Bool_t  isAssociated = kFALSE;
 577
 578     if (length > 350) {
 579       time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(track->P());
 580       if (time > 0) {
 581         beta = length / (2.99792457999999984e-02 * time);
 582         gamma = 1/TMath::Sqrt(1 - beta*beta);
 583         mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
 584       }
 585     }
 586     //
 587     // fill tree and print candidates (for short list)
 588     //
 589     Float_t cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.1931,
 590                                                            31.9806,
 591                                                            5.04114e-11,
 592                                                            2.13096,
 593                                                            2.38541);
 594     if (fMCtrue) cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),0.9931,
 595                                                                 31.9806,
 596                                                                 5.04114e-11,
 597                                                                 2.13096,
 598                                                                 2.38541);
 599     if (eta < 0.8 && tpcSignal > 120 && tpcSignal > cut && tpcSignal < 1000 && track->GetTPCsignalN() > 60 && dcaZ < 15 && dcaXY < 15 && ptot > 1.0 && ptot < 20) {
 600       //
 601     cout << "AntiAlphaEvent" << " "
 602            << AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree()->GetCurrentFile()->GetName() << " "
 603            << "event number in file: " << fESD->GetEventNumberInFile()
 604            << " Index " << iTracks
 605            << " ptot: " << ptot
 606            << " sig: " << tpcSignal <<endl;
 607       //
 608       fillTree = kTRUE;
 609       //
 610       TrkPtot[itrk] = track->P();
 611       TPCPtot[itrk] = ptot;
 612       DeDx[itrk] = tpcSignal;
 613       DCAZ[itrk] = dcaZ;
 614       TPCNsignal[itrk] = track->GetTPCsignalN();
 615       ITSnCluster[itrk] = track->GetNcls(0);
 616       Sign[itrk] = sign;
 617       DCAXY[itrk] = dcaXY;
 618       Mass[itrk] = mass;
 619       //
 620       ITSsignal[itrk] = track->GetITSsignal();
 621       SharedClusters[itrk] = shared.CountBits();
 622       sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fFileName);
 623       //fFileName[itrk] = fInputHandler->GetTree()->GetCurrentFile()->GetName();
 624       fEventNumber[itrk] = fESD->GetEventNumberInFile();
 625       //
 626       if(status&AliESDtrack::kITSrefit)
 627         ITSRefit[itrk] = 1;
 628       else ITSRefit[itrk] = 0;
 629       //
 630       if (time < 99998) {
 631         TOFRefit[itrk] = 1;
 632       } else {
 633         TOFRefit[itrk] = 0;
 634       }
 635       //itrk++;
 636       //
 637       hasTOFout = status&AliESDtrack::kTOFout;
 638       hasTOFtime  = status&AliESDtrack::kTIME;
 639       //
 640       TOFtime[itrk] = hasTOFtime;
 641       TOFout[itrk]  = hasTOFout;
 642
 643       if (fMCtrue){ //associated
 644
 645         Int_t label  = track->GetLabel();
 646         TParticle *tparticle = stack->Particle(TMath::Abs(label));
 647
 648         Bool_t isPrimary = stack->IsPhysicalPrimary(TMath::Abs(label));
 649         Bool_t isSecondary = stack->IsSecondaryFromMaterial(TMath::Abs(label));
 650
 651         Long_t pdgCode = tparticle->GetPdgCode();
 652         Double_t pT =(track->Pt())*2;
 653
 654         if(pdgCode == 1000020040)
 655           {
 656             fHistHelium4PtAso->Fill(pT);
 657             if(isPrimary) fHistHelium4PtAsoPrim->Fill(pT);
 658             if(isSecondary)  fHistHelium4PtAsoSec->Fill(pT);
 659             isAssociated = kTRUE;
 660           }
 661
 662         if(pdgCode == -1000020040)
 663           {
 664             fHistAntiHelium4PtAso->Fill(pT);
 665             isAssociated = kTRUE;
 666           }
 667
 668       }
 669
 670       fAssociated[itrk] = isAssociated;
 671       fTrackPt[itrk] = track->Pt();
 672
 673       itrk++;
 674     }
 675     //
 676     // do pid fill histogram for raw ratios
 677     //
 678     //                       (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
 679     Int_t id = -1;
 680     if (ptot > 0.2 && TMath::Abs(tpcSignal - expSignalProton)/expSignalProton < 0.2) id = 0;
 681     if (ptot > 0.3 && TMath::Abs(tpcSignal - expSignalDeuteron)/expSignalDeuteron < 0.2) id = 1;
 682     if (ptot > 0.7 && TMath::Abs(tpcSignal - expSignalTriton)/expSignalTriton < 0.2) id = 2;
 683     if (ptot > 0.5 && (tpcSignal - expSignalHelium3)/expSignalHelium3 > -0.1 &&  (tpcSignal - expSignalHelium3)/expSignalHelium3 < 0.2) id = 3;
 684     //
 685     Double_t vecAntiAlpha[4] = {dcaXYsign, sign, id, ptotInc};
 686     if (id != -1 && tpcSignal > 120) fAntiAlpha->Fill(vecAntiAlpha);
 687     //
 688     // fill final histograms
 689     //
 690     if (!fESDtrackCutsSharp->AcceptTrack(track) || shared.CountBits() > 1 || track->GetTPCsignalN() < 80 || track->GetKinkIndex(0) != 0 || track->GetTPCNclsIter1() < 80) continue;
 691     //
 692     fHistDEDx->Fill(ptot,track->GetTPCsignal(), sign);
 693     if (TMath::Abs(tpcSignal - expSignalHelium4)/expSignalHelium4 < 0.2 && time < 99998) fHistTOF3D->Fill(mass,1,sign);
 694     //
 695     // remove overlapping tracks with special dE/dx cut
 696     //
 697     //if (signal[1] < 6) continue;
 698     //if (signal[0]/signal[1] > 1.6 || signal[2]/signal[1] > 1.6 || signal[3]/signal[1] > 1.3) continue;
 699     //
 700     if(sign<0) {
 701       fHistDeDx->Fill(ptot, track->GetTPCsignal());
 702       if (track->GetTPCsignalN() > 100 &&
 703           TMath::Abs(track->Eta()) < 0.8 &&
 704           signal[3]/signal[1] > 0.6 &&
 705           signal[0]/signal[1] > 0.5 &&
 706           signal[3]/signal[1] < 1.2 &&
 707           track->GetTPCNclsIter1() > 70 &&
 708           track->GetTPCnclsS() < 10) {
 709         fHistDeDxSharp->Fill(ptot, track->GetTPCsignal());
 710       }
 711     }
 712     //
 713     // dE/dx for specific regions
 714     //
 715     for(Int_t iSig = 0; iSig < 4; iSig++) {
 716       if (signal[1] > 6) fHistDeDxRegion->Fill(ptot,signal[iSig]/signal[1],iSig);
 717     }
 718     //
 719     // alpha TOF plot
 720     //
 721     if((track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 > -0.15 && (track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 < 0.15) {
 722       //TOF
 723       hasTOFout  = status&AliESDtrack::kTOFout;
 724       hasTOFtime = status&AliESDtrack::kTIME;
 725       Bool_t hasTOF     = kFALSE;
 726
 727       if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
 728       //
 729       if (length < 350.) hasTOF = kFALSE;
 730
 731       Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
 732       //
 733       if (length > 350. &&  hasTOF == kTRUE && ptot < 4) {
 734         time = track->GetTOFsignal() - time0;
 735         if (time > 0) {
 736           beta = length / (2.99792457999999984e-02 * time);
 737           if (beta < 0.975) {
 738             gamma = 1/TMath::Sqrt(1 - beta*beta);
 739             mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
 740             if (TMath::Sqrt(track->GetTOFsignalDz()*track->GetTOFsignalDz() + track->GetTOFsignalDx()*track->GetTOFsignalDx()) < 5.) {
 741               fHistAlpha->Fill(mass*mass);
 742               if (mass*mass > 3. && mass*mass < 4.) {
 743                 fHistAlphaSignal->Fill(mass*mass);
 744                 fGraphAlphaSignal->SetPoint(fNCounter, ptot, track->GetTPCsignal());
 745                 fNCounter++;
 746               }
 747               fHistTOFnuclei->Fill(ptot,beta);
 748             }
 749           }
 750         }
 751       }
 752     }
 753
 754   }//end loop over tracks
 755
 756   if (fillTree) fTree->Fill();
 757
 758   // Post output data.
 759   PostData(1, fOutputContainer);
 760   PostData(2, fTree);
 761 }
 762
 763 //________________________________________________________________________
 764 void AliAnalysisTaskAntiHe4::Terminate(const Option_t *)
 765 {
 766   // Draw result to the screen
 767   // Called once at the end of the query
 768
 769   //get output data and darw 'fHistPt'
 770   if (!GetOutputData(0)) return;
 771
 772 }
 773
 774
 775 //________________________________________________________________________
 776 void AliAnalysisTaskAntiHe4::BinLogAxis(const THn *h, Int_t axisNumber) {
 777   //
 778   // Method for the correct logarithmic binning of histograms
 779   //
 780   TAxis *axis = h->GetAxis(axisNumber);
 781   int bins = axis->GetNbins();
 782
 783   Double_t from = axis->GetXmin();
 784   Double_t to = axis->GetXmax();
 785   Double_t *newBins = new Double_t[bins + 1];
 786
 787   newBins[0] = from;
 788   Double_t factor = pow(to/from, 1./bins);
 789
 790   for (int i = 1; i <= bins; i++) {
 791    newBins[i] = factor * newBins[i-1];
 792   }
 793   axis->Set(bins, newBins);
 794   delete [] newBins;
 795
 796 }
 797
 798 //________________________________________________________________________
 799 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH3 *h, Int_t axisNumber) {
 800   //
 801   // Method for the correct logarithmic binning of histograms
 802   //
 803   TAxis *axis = h->GetXaxis();
 804   if (axisNumber == 1) axis = h->GetYaxis();
 805   if (axisNumber == 2) axis = h->GetZaxis();
 806   int bins = axis->GetNbins();
 807
 808   Double_t from = axis->GetXmin();
 809   Double_t to = axis->GetXmax();
 810   Double_t *newBins = new Double_t[bins + 1];
 811
 812   newBins[0] = from;
 813   Double_t factor = pow(to/from, 1./bins);
 814
 815   for (int i = 1; i <= bins; i++) {
 816    newBins[i] = factor * newBins[i-1];
 817   }
 818   axis->Set(bins, newBins);
 819   delete [] newBins;
 820
 821 }
 822 //________________________________________________________________________
 823 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH1 *h) {
 824   //
 825   // Method for the correct logarithmic binning of histograms
 826   //
 827   TAxis *axis = h->GetXaxis();
 828   int bins = axis->GetNbins();
 829
 830   Double_t from = axis->GetXmin();
 831   Double_t to = axis->GetXmax();
 832   Double_t *newBins = new Double_t[bins + 1];
 833
 834   newBins[0] = from;
 835   Double_t factor = pow(to/from, 1./bins);
 836
 837   for (int i = 1; i <= bins; i++) {
 838     newBins[i] = factor * newBins[i-1];
 839   }
 840   axis->Set(bins, newBins);
 841   delete [] newBins;
 842
 843 }
 844 //_____________________________________________________________________________
 845 Int_t AliAnalysisTaskAntiHe4::Initialize() {
 846
 847
 848   // -- Reset Event
 849   // ----------------
 850   ResetEvent();
 851
 852   return 0;
 853 }
 854 //________________________________________________________________________
 855 Int_t AliAnalysisTaskAntiHe4::SetupEvent() {
 856   // Setup Reading of event
 857
 858   ResetEvent();
 859   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 860   // -- Get Trigger
 861   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 862
 863   //Bool_t isTriggered = IsTriggered();
 864   IsTriggered();
 865   return 0;
 866 }
 867 //_____________________________________________________________________________
 868 void AliAnalysisTaskAntiHe4::ResetEvent() {
 869   // Reset event
 870   // -- Reset QA
 871   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 872     fTriggerFired[ii] = kFALSE;
 873
 874   return;
 875 }
 876 //________________________________________________________________________
 877 Bool_t AliAnalysisTaskAntiHe4::IsTriggered() {
 878   // Check if Event is triggered and fill Trigger Histogram
 879
 880   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
 881   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
 882   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
 883   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
 884   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
 885
 886   Bool_t isTriggered = kFALSE;
 887
 888   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 889     if(fTriggerFired[ii]) {
 890       isTriggered = kTRUE;
 891       fHistTriggerStat->Fill(ii);
 892     }
 893   }
 894
 895   return isTriggered;
 896 }
 897 //________________________________________________________________________
 898 void AliAnalysisTaskAntiHe4::SetBBParameters(Int_t runNumber){
 899
 900   if(runNumber < 166500){ //LHC10h
 901
 902     fBBParametersLightParticles[0]   = 1.45802;
 903     fBBParametersLightParticles[1]   = 27.4992;
 904     fBBParametersLightParticles[2]   = 4.00313e-15;
 905     fBBParametersLightParticles[3]   = 2.48485;
 906     fBBParametersLightParticles[4]   = 8.31768;
 907
 908     fBBParametersNuclei[0]  = 1.69155;
 909     fBBParametersNuclei[1]  = 27.4992;
 910     fBBParametersNuclei[2]  = 4.00313e-15;
 911     fBBParametersNuclei[3]  = 2.48485;
 912     fBBParametersNuclei[4]  = 8.31768;
 913
 914   }
 915
 916   if(runNumber > 166500){ //LHC11h
 917
 918     fBBParametersLightParticles[0]   = 4.69637;//1.11243;
 919     fBBParametersLightParticles[1]   = 7.51827;//26.1144;
 920     fBBParametersLightParticles[2]   = 0.0183746;//4.00313e-15;
 921     fBBParametersLightParticles[3]   = 2.60;//2.72969;
 922     fBBParametersLightParticles[4]   = 2.7;//9.15038;
 923
 924     fBBParametersNuclei[0]  = 1.4906;
 925     fBBParametersNuclei[1]  = 27.9758;
 926     fBBParametersNuclei[2]  = 4.00313e-15;
 927     fBBParametersNuclei[3]  = 2.50804;
 928     fBBParametersNuclei[4]  = 8.31768;
 929
 930   }
 931 }
 932 //______________________________________________________________________________
 933 void AliAnalysisTaskAntiHe4::MCGenerated(AliStack* stack)
 934 {
 935
 936   // Monte Carlo for genenerated particles
 937   if (fMCtrue) //MC loop
 938     {
 939
 940       Int_t stackN = 0;
 941
 942       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
 943         {
 944
 945           Bool_t isPrimary = stack->IsPhysicalPrimary(stackN);
 946           Bool_t isSecondary = stack->IsSecondaryFromMaterial(stackN);
 947
 948           const TParticle *tparticle = stack->Particle(stackN);
 949           Long_t pdgCode = tparticle->GetPdgCode();
 950
 951           Double_t pTGen = tparticle->Pt();
 952           Double_t eta = tparticle->Eta();
 953           //check which particle it is
 954
 955           //Alpha
 956           if(pdgCode == 1000020040)
 957             {
 958               fHistHelium4PtGen->Fill(pTGen);
 959               if(isPrimary) fHistHelium4PtGenPrim->Fill(pTGen);
 960               if(isSecondary) fHistHelium4PtGenSec->Fill(pTGen);
 961               if(eta < 0.8)fHistHelium4PtGenEta->Fill(pTGen);
 962               if(isPrimary && eta < 0.8)fHistHelium4PtGenPrimEta->Fill(pTGen);
 963             }
 964
 965           //Anti-Alpha
 966           if(pdgCode == -1000020040)
 967             {
 968               fHistAntiHelium4PtGen->Fill(pTGen);
 969               if(isPrimary) fHistAntiHelium4PtGenPrim->Fill(pTGen);
 970               if(isSecondary) fHistAntiHelium4PtGenSec->Fill(pTGen);
 971               if(eta < 0.8)fHistAntiHelium4PtGenEta->Fill(pTGen);
 972             }
 973
 974
 975         }//end loop over stack
 976
 977
 978     }//end MC
 979 }