PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskAntiHe4.cxx

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