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   //
 392   fTree->Branch("fTPCnCluster",fTPCnCluster,"fTPCnCluster[fItrk]/s");
 393   fTree->Branch("fTPCNsignal",fTPCNsignal,"fTPCNsignal[fItrk]/s");
 394   fTree->Branch("fChi2PerClusterTPC",fChi2PerClusterTPC,"fChi2PerClusterTPC[fItrk]/D");
 395   fTree->Branch("fTPCRefit",fTPCRefit,"fTPCRefit[fItrk]/O");
 396   fTree->Branch("fTPCsignal0",fTPCsignal0,"fTPCsignal0[fItrk]/D");
 397   fTree->Branch("fTPCsignal1",fTPCsignal1,"fTPCsignal1[fItrk]/D");
 398   fTree->Branch("fTPCsignal2",fTPCsignal2,"fTPCsignal2[fItrk]/D");
 399   fTree->Branch("fTPCsignal3",fTPCsignal3,"fTPCsignal3[fItrk]/D");
 400   fTree->Branch("fTPCSharedClusters",fTPCSharedClusters,"fTPCSharedClusters[fItrk]/I");
 401   fTree->Branch("fTPCNclsIter1",fTPCNclsIter1,"fTPCNclsIter1[fItrk]/s");
 402   //
 403   fTree->Branch("fITSsignal",fITSsignal,"fITSsignal[fItrk]/D");
 404   fTree->Branch("fITSnCluster",fITSnCluster,"fITSnCluster[fItrk]/I");
 405   fTree->Branch("fChi2PerClusterITS",fChi2PerClusterITS,"fChi2PerClusterITS[fItrk]/D");
 406   fTree->Branch("fITSRefit",fITSRefit,"fITSRefit[fItrk]/O");
 407   //
 408   fTree->Branch("fTOFtime",fTOFtime,"fTOFtime[fItrk]/O");
 409   fTree->Branch("fTOFRefit",fTOFRefit,"fTOFRefit[fItrk]/O");
 410   fTree->Branch("fTOFout",fTOFout,"fTOFout[fItrk]/O");
 411   fTree->Branch("fTOFsignalDz",fTOFsignalDz,"fTOFsignalDz[fItrk]/D");
 412   fTree->Branch("fTOFsignalDx",fTOFsignalDx,"fTOFsignalDx[fItrk]/D");
 413   //
 414   fTree->Branch("fDCAXY",fDCAXY,"fDCAXY[fItrk]/F");
 415   fTree->Branch("fDCAZ",fDCAZ,"fDCAZ[fItrk]/F");
 416   //
 417   fTree->Branch("fTrkPtot",fTrkPtot,"fTrkPtot[fItrk]/D");
 418   fTree->Branch("fTPCPtot",fTPCPtot,"fTPCPtot[fItrk]/D");
 419   fTree->Branch("fTrackPt",fTrackPt,"fTrackPt[fItrk]/D");
 420   fTree->Branch("fDeDx",fDeDx,"fDeDx[fItrk]/D");
 421   fTree->Branch("fSign",fSign,"fSign[fItrk]/D");
 422   fTree->Branch("fMass",fMass,"Mass[fItrk]/F");
 423   //
 424   fTree->Branch("fAssociated",fAssociated,"fAssociated[fItrk]/O");
 425
 426 }
 427
 428 //________________________________________________________________________
 429 void AliAnalysisTaskAntiHe4::UserExec(Option_t *)
 430 {
 431   // Main loop
 432   // Called for each event
 433
 434   //get Event-Handler for the trigger information
 435   fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 436   if (!fEventHandler) {
 437     AliError("Could not get InputHandler");
 438     //return -1;
 439     return;
 440   }
 441
 442   // Monte Carlo
 443   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 444   if (!eventHandler){
 445     //Printf("ERROR: Could not retrieve MC event handler");
 446     fMCtrue = kFALSE;
 447   }
 448
 449   AliMCEvent* mcEvent = 0x0;
 450   AliStack* stack = 0x0;
 451   if (eventHandler) mcEvent = eventHandler->MCEvent();
 452   if (!mcEvent){
 453     //Printf("ERROR: Could not retrieve MC event");
 454     if (fMCtrue) return;
 455   }
 456
 457   if (fMCtrue){
 458     stack = mcEvent->Stack();
 459     if (!stack) return;
 460   }
 461
 462   //look for the generated particles
 463   MCGenerated(stack);
 464
 465   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 466   if (!fESD) {
 467     //Printf("ERROR: fESD not available");
 468     return;
 469   }
 470
 471   if (SetupEvent() < 0) {
 472     PostData(1, fOutputContainer);
 473     return;
 474   }
 475
 476   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 477   if(vertex->GetNContributors()<1) {
 478     // SPD vertex
 479     vertex = fESD->GetPrimaryVertexSPD();
 480     if(vertex->GetNContributors()<1) {
 481       PostData(1, fOutputContainer);
 482       return;
 483     }
 484
 485   }
 486   // check if event is selected by physics selection class
 487   //
 488   Bool_t isSelected = kFALSE;
 489   isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 490   if (!isSelected || TMath::Abs(vertex->GetZv()) > 10) {
 491     PostData(1, fOutputContainer);
 492     return;
 493   }
 494
 495   //
 496   //centrality selection in PbPb
 497   //
 498   Int_t centralityClass10 = -1;
 499   Int_t centralityPercentile = -1;
 500   //
 501   if (fESD->GetEventSpecie() == 4) { //species == 4 == PbPb
 502     // PbPb
 503     AliCentrality *esdCentrality = fESD->GetCentrality();
 504     centralityClass10 = esdCentrality->GetCentralityClass10("V0M"); // centrality percentile determined with V0
 505     centralityPercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0
 506     if(!fMCtrue){
 507       if (centralityPercentile < 0. || centralityPercentile > 80. ) return; //select only events with centralities between 0 and 80 %
 508     }
 509   }
 510   //
 511   fHistCentralityClass10->Fill(centralityClass10);
 512   fHistCentralityPercentile->Fill(centralityPercentile);
 513   //
 514   if(fTriggerFired[0] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(0);
 515   if(fTriggerFired[1] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(1);
 516   if(fTriggerFired[2] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(2);
 517   if(fTriggerFired[3] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(3);
 518   if(fTriggerFired[4] == kTRUE) fHistTriggerStatAfterEventSelection->Fill(4);
 519   //
 520   //
 521   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 522   if (!fESDpid) {
 523     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
 524     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
 525   }
 526   //
 527   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for th // for Anti-Alpha
 528   fEvnt =  fESD->GetEventNumberInFile();
 529   sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fName);
 530   fItrk = 0;
 531   //
 532   Int_t runNumber = 0;
 533   runNumber = fESD->GetRunNumber();
 534   //
 535   Bool_t fillTree = kFALSE;
 536   // Track loop to fill the spectram
 537   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 538
 539     fEventNumber[fItrk] = -1;
 540
 541     fEta[fItrk] = -2;
 542     fTPCNsignal[fItrk] = -1;
 543     fTPCnCluster[fItrk] = -1;
 544     fChi2PerClusterTPC[fItrk] = -1;
 545     fTPCRefit[fItrk] = kFALSE;
 546     fTPCsignal0[fItrk] = -1;
 547     fTPCsignal1[fItrk] = -1;
 548     fTPCsignal2[fItrk] = -1;
 549     fTPCsignal3[fItrk] = -1;
 550     fTPCSharedClusters[fItrk] = -1;
 551     fTPCNclsIter1[fItrk] = -1;
 552
 553     fITSsignal[fItrk] = -1;
 554     fITSnCluster[fItrk] = -1;
 555     fChi2PerClusterITS[fItrk] = -1;
 556     fITSRefit[fItrk] = kFALSE;
 557
 558     fTOFRefit[fItrk] = kFALSE;
 559     fTOFtime[fItrk] = kFALSE;
 560     fTOFout[fItrk] = kFALSE;
 561     fTOFsignalDz[fItrk] = -1;
 562     fTOFsignalDx[fItrk] = -1;
 563
 564     fDCAZ[fItrk] = -1;
 565     fDCAXY[fItrk] = -1;
 566
 567     fTrkPtot[fItrk] = -1;
 568     fTPCPtot[fItrk] = -1;
 569     fTrackPt[fItrk] = -1;
 570     fDeDx[fItrk] = -1;
 571     fSign[fItrk] = -2;
 572     fMass[fItrk] = -1;
 573
 574     fAssociated[fItrk] = kFALSE;
 575
 576     AliESDtrack* track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTracks));
 577     if (!fESDtrackCuts->AcceptTrack(track)) continue;
 578     //
 579     Double_t nClustersTPCPID=track->GetTPCsignalN();
 580     if(nClustersTPCPID < 60) continue;
 581     //
 582     if (!track->GetInnerParam()) continue;
 583     //
 584     Double32_t signal[4] = {0,0,0,0};
 585     Char_t ncl[3];
 586     Char_t nrows[3];
 587     if (track->GetTPCdEdxInfo()) {
 588       track->GetTPCdEdxInfo()->GetTPCSignalRegionInfo(signal,ncl,nrows);
 589     }
 590     //
 591     Double_t ptot = track->GetInnerParam()->GetP();
 592     Double_t ptotInc = track->GetP(); // total momentum of the incoming particle
 593     Double_t sign = track->GetSign();
 594     //
 595     Double_t eta = TMath::Abs(track->Eta());
 596     TBits shared = track->GetTPCSharedMap();
 597
 598     track->GetImpactParameters(dca, cov);
 599     Float_t dcaXY = TMath::Sqrt(dca[0]*dca[0]);
 600     Float_t dcaXYsign = dca[0];
 601     Float_t dcaZ  = TMath::Sqrt(dca[1]*dca[1]);
 602     //
 603     //
 604     Double_t tpcSignal = track->GetTPCsignal();
 605     //
 606     //PID via specific energy loss in the TPC
 607     //define the arrays for the Bethe-Bloch-Parameters
 608     //
 609
 610     SetBBParameters(runNumber);
 611
 612     //define expected signals for the various species
 613     Double_t expSignalProton = AliExternalTrackParam::BetheBlochAleph(ptot/0.93827,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 614     Double_t expSignalDeuteron = AliExternalTrackParam::BetheBlochAleph(ptot/1.875612,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 615     Double_t expSignalTriton = AliExternalTrackParam::BetheBlochAleph(ptot/2.808921,fBBParametersLightParticles[0],fBBParametersLightParticles[1],fBBParametersLightParticles[2],fBBParametersLightParticles[3],fBBParametersLightParticles[4]);
 616
 617     Double_t expSignalHelium3 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/2.80941,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
 618     Double_t expSignalHelium4 = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/3.728401,fBBParametersNuclei[0],fBBParametersNuclei[1],fBBParametersNuclei[2],fBBParametersNuclei[3],fBBParametersNuclei[4]);
 619
 620     //
 621     Float_t mass = 0;
 622     Float_t time = -1;
 623     Float_t beta = 0;
 624     Float_t gamma = -1;
 625     Bool_t  hasTOFout  = kFALSE;
 626     Bool_t  hasTOFtime = kFALSE;
 627     Float_t length = track->GetIntegratedLength();
 628     UInt_t  status = track->GetStatus();
 629     Bool_t  isAssociated = kFALSE;
 630
 631     if (length > 350) {
 632       time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(track->P());
 633       if (time > 0) {
 634         beta = length / (2.99792457999999984e-02 * time);
 635         gamma = 1/TMath::Sqrt(1 - beta*beta);
 636         mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
 637       }
 638     }
 639     //
 640     // fill tree and print candidates (for short list)
 641     //
 642     Float_t cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.1931,
 643                                                            31.9806,
 644                                                            5.04114e-11,
 645                                                            2.13096,
 646                                                            2.38541);
 647     if (fMCtrue) cut = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),0.9931,
 648                                                                 31.9806,
 649                                                                 5.04114e-11,
 650                                                                 2.13096,
 651                                                                 2.38541);
 652     if (eta < 1.0 && tpcSignal > 120 && tpcSignal > cut && tpcSignal < 1000 && track->GetTPCsignalN() > 60 && dcaZ < 15 && dcaXY < 15 && ptot > 1.0 && ptot < 20) {
 653       //
 654     cout << "AntiAlphaEvent" << " "
 655            << AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetTree()->GetCurrentFile()->GetName() << " "
 656            << "event number in file: " << fESD->GetEventNumberInFile()
 657            << " Index " << iTracks
 658            << " ptot: " << ptot
 659            << " sig: " << tpcSignal <<endl;
 660       //
 661       fillTree = kTRUE;
 662       //
 663
 664       sscanf(fInputHandler->GetTree()->GetCurrentFile()->GetName(),"%s", fFileName);
 665       fEventNumber[fItrk] = fESD->GetEventNumberInFile();
 666
 667       fEta[fItrk] = eta;
 668       fKinkIndex[fItrk] = track->GetKinkIndex(0);
 669
 670       fTPCNsignal[fItrk] = track->GetTPCsignalN();
 671       fTPCnCluster[fItrk] = track->GetTPCNcls();
 672       fChi2PerClusterTPC[fItrk] = track->GetTPCchi2()/fTPCnCluster[fItrk];
 673       if(status&AliESDtrack::kTPCrefit)
 674         fTPCRefit[fItrk] = kTRUE;
 675       else fTPCRefit[fItrk] = kFALSE;
 676       fTPCsignal0[fItrk] = signal[0];
 677       fTPCsignal1[fItrk] = signal[1];
 678       fTPCsignal2[fItrk] = signal[2];
 679       fTPCsignal3[fItrk] = signal[3];
 680       fTPCSharedClusters[fItrk] = shared.CountBits();
 681       fTPCNclsIter1[fItrk] = track->GetTPCNclsIter1();
 682
 683       fITSsignal[fItrk] = track->GetITSsignal();
 684       fITSnCluster[fItrk] = track->GetNcls(0);
 685       fChi2PerClusterITS[fItrk] = track->GetITSchi2()/fITSnCluster[fItrk];
 686       if(status&AliESDtrack::kITSrefit)
 687         fITSRefit[fItrk] = kTRUE;
 688       else fITSRefit[fItrk] = kFALSE;
 689
 690
 691       if(status&AliESDtrack::kITSrefit)
 692         fITSRefit[fItrk] = kTRUE;
 693       else fITSRefit[fItrk] = kFALSE;
 694       hasTOFout = status&AliESDtrack::kTOFout;
 695       hasTOFtime  = status&AliESDtrack::kTIME;
 696       fTOFtime[fItrk] = hasTOFtime;
 697       fTOFout[fItrk]  = hasTOFout;
 698       fTOFsignalDz[fItrk] = track->GetTOFsignalDz();
 699       fTOFsignalDx[fItrk] = track->GetTOFsignalDx();
 700
 701       fDCAZ[fItrk] = dcaXY;
 702       fDCAXY[fItrk] = dcaZ;
 703
 704       fTrkPtot[fItrk] = track->P();
 705       fTPCPtot[fItrk] = ptot;
 706       fTrackPt[fItrk] = track->Pt();
 707       fDeDx[fItrk] = tpcSignal;
 708       fSign[fItrk] = sign;
 709       fMass[fItrk] = mass;
 710
 711       if (fMCtrue){ //associated
 712
 713         Int_t label  = track->GetLabel();
 714         TParticle *tparticle = stack->Particle(TMath::Abs(label));
 715
 716         Bool_t isPrimary = stack->IsPhysicalPrimary(TMath::Abs(label));
 717         Bool_t isSecondary = stack->IsSecondaryFromMaterial(TMath::Abs(label));
 718
 719         Long_t pdgCode = tparticle->GetPdgCode();
 720         Double_t pT =(track->Pt())*2;
 721
 722         if(pdgCode == 1000020040)
 723           {
 724             fHistHelium4PtAso->Fill(pT);
 725             if(isPrimary) fHistHelium4PtAsoPrim->Fill(pT);
 726             if(isSecondary)  fHistHelium4PtAsoSec->Fill(pT);
 727             isAssociated = kTRUE;
 728           }
 729
 730         if(pdgCode == -1000020040)
 731           {
 732             fHistAntiHelium4PtAso->Fill(pT);
 733             isAssociated = kTRUE;
 734           }
 735
 736       }
 737
 738       fAssociated[fItrk] = isAssociated;
 739
 740       fItrk++;
 741     }
 742     //
 743     // do pid fill histogram for raw ratios
 744     //
 745     //                       (0.) dca, (1.) sign, (2.) particle Type, (3.) p_tot
 746     Int_t id = -1;
 747     if (ptot > 0.2 && TMath::Abs(tpcSignal - expSignalProton)/expSignalProton < 0.2) id = 0;
 748     if (ptot > 0.3 && TMath::Abs(tpcSignal - expSignalDeuteron)/expSignalDeuteron < 0.2) id = 1;
 749     if (ptot > 0.7 && TMath::Abs(tpcSignal - expSignalTriton)/expSignalTriton < 0.2) id = 2;
 750     if (ptot > 0.5 && (tpcSignal - expSignalHelium3)/expSignalHelium3 > -0.1 &&  (tpcSignal - expSignalHelium3)/expSignalHelium3 < 0.2) id = 3;
 751     //
 752     Double_t vecAntiAlpha[4] = {dcaXYsign, sign, id, ptotInc};
 753     if (id != -1 && tpcSignal > 120) fAntiAlpha->Fill(vecAntiAlpha);
 754     //
 755     // fill final histograms
 756     //
 757     if (!fESDtrackCutsSharp->AcceptTrack(track) || shared.CountBits() > 1 || track->GetTPCsignalN() < 80 || track->GetKinkIndex(0) != 0 || track->GetTPCNclsIter1() < 80) continue;
 758     //
 759     fHistDEDx->Fill(ptot,track->GetTPCsignal(), sign);
 760     if (TMath::Abs(tpcSignal - expSignalHelium4)/expSignalHelium4 < 0.2 && time < 99998) fHistTOF3D->Fill(mass,1,sign);
 761     //
 762     // remove overlapping tracks with special dE/dx cut
 763     //
 764     //if (signal[1] < 6) continue;
 765     //if (signal[0]/signal[1] > 1.6 || signal[2]/signal[1] > 1.6 || signal[3]/signal[1] > 1.3) continue;
 766     //
 767     if(sign<0) {
 768       fHistDeDx->Fill(ptot, track->GetTPCsignal());
 769       if (track->GetTPCsignalN() > 100 &&
 770           TMath::Abs(track->Eta()) < 1.0 &&
 771           signal[3]/signal[1] > 0.6 &&
 772           signal[0]/signal[1] > 0.5 &&
 773           signal[3]/signal[1] < 1.2 &&
 774           track->GetTPCNclsIter1() > 70 &&
 775           track->GetTPCnclsS() < 10) {
 776         fHistDeDxSharp->Fill(ptot, track->GetTPCsignal());
 777       }
 778     }
 779     //
 780     // dE/dx for specific regions
 781     //
 782     for(Int_t iSig = 0; iSig < 4; iSig++) {
 783       if (signal[1] > 6) fHistDeDxRegion->Fill(ptot,signal[iSig]/signal[1],iSig);
 784     }
 785     //
 786     // alpha TOF plot
 787     //
 788     if((track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 > -0.15 && (track->GetTPCsignal()-expSignalHelium4)/expSignalHelium4 < 0.15) {
 789       //TOF
 790       hasTOFout  = status&AliESDtrack::kTOFout;
 791       hasTOFtime = status&AliESDtrack::kTIME;
 792       Bool_t hasTOF     = kFALSE;
 793
 794       if (hasTOFout && hasTOFtime) hasTOF = kTRUE;
 795       //
 796       if (length < 350.) hasTOF = kFALSE;
 797
 798       Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
 799       //
 800       if (length > 350. &&  hasTOF == kTRUE && ptot < 4) {
 801         time = track->GetTOFsignal() - time0;
 802         if (time > 0) {
 803           beta = length / (2.99792457999999984e-02 * time);
 804           if (beta < 0.975) {
 805             gamma = 1/TMath::Sqrt(1 - beta*beta);
 806             mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
 807             if (TMath::Sqrt(track->GetTOFsignalDz()*track->GetTOFsignalDz() + track->GetTOFsignalDx()*track->GetTOFsignalDx()) < 5.) {
 808               fHistAlpha->Fill(mass*mass);
 809               if (mass*mass > 3. && mass*mass < 4.) {
 810                 fHistAlphaSignal->Fill(mass*mass);
 811                 fGraphAlphaSignal->SetPoint(fNCounter, ptot, track->GetTPCsignal());
 812                 fNCounter++;
 813               }
 814               fHistTOFnuclei->Fill(ptot,beta);
 815             }
 816           }
 817         }
 818       }
 819     }
 820
 821   }//end loop over tracks
 822
 823   if (fillTree) fTree->Fill();
 824
 825   // Post output data.
 826   PostData(1, fOutputContainer);
 827   PostData(2, fTree);
 828 }
 829
 830 //________________________________________________________________________
 831 void AliAnalysisTaskAntiHe4::Terminate(const Option_t *)
 832 {
 833   // Draw result to the screen
 834   // Called once at the end of the query
 835
 836   //get output data and darw 'fHistPt'
 837   if (!GetOutputData(0)) return;
 838
 839 }
 840
 841
 842 //________________________________________________________________________
 843 void AliAnalysisTaskAntiHe4::BinLogAxis(const THn *h, Int_t axisNumber) {
 844   //
 845   // Method for the correct logarithmic binning of histograms
 846   //
 847   TAxis *axis = h->GetAxis(axisNumber);
 848   int bins = axis->GetNbins();
 849
 850   Double_t from = axis->GetXmin();
 851   Double_t to = axis->GetXmax();
 852   Double_t *newBins = new Double_t[bins + 1];
 853
 854   newBins[0] = from;
 855   Double_t factor = pow(to/from, 1./bins);
 856
 857   for (int i = 1; i <= bins; i++) {
 858    newBins[i] = factor * newBins[i-1];
 859   }
 860   axis->Set(bins, newBins);
 861   delete [] newBins;
 862
 863 }
 864
 865 //________________________________________________________________________
 866 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH3 *h, Int_t axisNumber) {
 867   //
 868   // Method for the correct logarithmic binning of histograms
 869   //
 870   TAxis *axis = h->GetXaxis();
 871   if (axisNumber == 1) axis = h->GetYaxis();
 872   if (axisNumber == 2) axis = h->GetZaxis();
 873   int bins = axis->GetNbins();
 874
 875   Double_t from = axis->GetXmin();
 876   Double_t to = axis->GetXmax();
 877   Double_t *newBins = new Double_t[bins + 1];
 878
 879   newBins[0] = from;
 880   Double_t factor = pow(to/from, 1./bins);
 881
 882   for (int i = 1; i <= bins; i++) {
 883    newBins[i] = factor * newBins[i-1];
 884   }
 885   axis->Set(bins, newBins);
 886   delete [] newBins;
 887
 888 }
 889 //________________________________________________________________________
 890 void AliAnalysisTaskAntiHe4::BinLogAxis(const TH1 *h) {
 891   //
 892   // Method for the correct logarithmic binning of histograms
 893   //
 894   TAxis *axis = h->GetXaxis();
 895   int bins = axis->GetNbins();
 896
 897   Double_t from = axis->GetXmin();
 898   Double_t to = axis->GetXmax();
 899   Double_t *newBins = new Double_t[bins + 1];
 900
 901   newBins[0] = from;
 902   Double_t factor = pow(to/from, 1./bins);
 903
 904   for (int i = 1; i <= bins; i++) {
 905     newBins[i] = factor * newBins[i-1];
 906   }
 907   axis->Set(bins, newBins);
 908   delete [] newBins;
 909
 910 }
 911 //_____________________________________________________________________________
 912 Int_t AliAnalysisTaskAntiHe4::Initialize() {
 913
 914
 915   // -- Reset Event
 916   // ----------------
 917   ResetEvent();
 918
 919   return 0;
 920 }
 921 //________________________________________________________________________
 922 Int_t AliAnalysisTaskAntiHe4::SetupEvent() {
 923   // Setup Reading of event
 924
 925   ResetEvent();
 926   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 927   // -- Get Trigger
 928   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 929
 930   //Bool_t isTriggered = IsTriggered();
 931   IsTriggered();
 932   return 0;
 933 }
 934 //_____________________________________________________________________________
 935 void AliAnalysisTaskAntiHe4::ResetEvent() {
 936   // Reset event
 937   // -- Reset QA
 938   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 939     fTriggerFired[ii] = kFALSE;
 940
 941   return;
 942 }
 943 //________________________________________________________________________
 944 Bool_t AliAnalysisTaskAntiHe4::IsTriggered() {
 945   // Check if Event is triggered and fill Trigger Histogram
 946
 947   if ((fEventHandler->IsEventSelected() & AliVEvent::kMB))          fTriggerFired[0] = kTRUE;
 948   if ((fEventHandler->IsEventSelected() & AliVEvent::kCentral))     fTriggerFired[1] = kTRUE;
 949   if ((fEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) fTriggerFired[2] = kTRUE;
 950   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      fTriggerFired[3] = kTRUE;
 951   if ((fEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      fTriggerFired[4] = kTRUE;
 952
 953   Bool_t isTriggered = kFALSE;
 954
 955   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 956     if(fTriggerFired[ii]) {
 957       isTriggered = kTRUE;
 958       fHistTriggerStat->Fill(ii);
 959     }
 960   }
 961
 962   return isTriggered;
 963 }
 964 //________________________________________________________________________
 965 void AliAnalysisTaskAntiHe4::SetBBParameters(Int_t runNumber){
 966
 967   if(runNumber < 166500){ //LHC10h
 968
 969     fBBParametersLightParticles[0]   = 1.45802;
 970     fBBParametersLightParticles[1]   = 27.4992;
 971     fBBParametersLightParticles[2]   = 4.00313e-15;
 972     fBBParametersLightParticles[3]   = 2.48485;
 973     fBBParametersLightParticles[4]   = 8.31768;
 974
 975     fBBParametersNuclei[0]  = 1.69155;
 976     fBBParametersNuclei[1]  = 27.4992;
 977     fBBParametersNuclei[2]  = 4.00313e-15;
 978     fBBParametersNuclei[3]  = 2.48485;
 979     fBBParametersNuclei[4]  = 8.31768;
 980
 981   }
 982
 983   if(runNumber > 166500){ //LHC11h
 984
 985     fBBParametersLightParticles[0]   = 4.69637;//1.11243;
 986     fBBParametersLightParticles[1]   = 7.51827;//26.1144;
 987     fBBParametersLightParticles[2]   = 0.0183746;//4.00313e-15;
 988     fBBParametersLightParticles[3]   = 2.60;//2.72969;
 989     fBBParametersLightParticles[4]   = 2.7;//9.15038;
 990
 991     fBBParametersNuclei[0]  = 1.4906;
 992     fBBParametersNuclei[1]  = 27.9758;
 993     fBBParametersNuclei[2]  = 4.00313e-15;
 994     fBBParametersNuclei[3]  = 2.50804;
 995     fBBParametersNuclei[4]  = 8.31768;
 996
 997   }
 998 }
 999 //______________________________________________________________________________
1000 void AliAnalysisTaskAntiHe4::MCGenerated(AliStack* stack)
1001 {
1002
1003   // Monte Carlo for genenerated particles
1004   if (fMCtrue) //MC loop
1005     {
1006
1007       Int_t stackN = 0;
1008
1009       for(stackN = 0; stackN < stack->GetNtrack(); stackN++) //loop over stack
1010         {
1011
1012           Bool_t isPrimary = stack->IsPhysicalPrimary(stackN);
1013           Bool_t isSecondary = stack->IsSecondaryFromMaterial(stackN);
1014
1015           const TParticle *tparticle = stack->Particle(stackN);
1016           Long_t pdgCode = tparticle->GetPdgCode();
1017
1018           Double_t pTGen = tparticle->Pt();
1019           Double_t eta = tparticle->Eta();
1020           //check which particle it is
1021
1022           //Alpha
1023           if(pdgCode == 1000020040)
1024             {
1025               fHistHelium4PtGen->Fill(pTGen);
1026               if(isPrimary) fHistHelium4PtGenPrim->Fill(pTGen);
1027               if(isSecondary) fHistHelium4PtGenSec->Fill(pTGen);
1028               if(TMath::Abs(eta) < 1.0)fHistHelium4PtGenEta->Fill(pTGen);
1029               if(isPrimary && TMath::Abs(eta) < 1.0)fHistHelium4PtGenPrimEta->Fill(pTGen);
1030             }
1031
1032           //Anti-Alpha
1033           if(pdgCode == -1000020040)
1034             {
1035               fHistAntiHelium4PtGen->Fill(pTGen);
1036               if(isPrimary) fHistAntiHelium4PtGenPrim->Fill(pTGen);
1037               if(isSecondary) fHistAntiHelium4PtGenSec->Fill(pTGen);
1038               if(TMath::Abs(eta) < 1.0)fHistAntiHelium4PtGenEta->Fill(pTGen);
1039             }
1040
1041
1042         }//end loop over stack
1043
1044
1045     }//end MC
1046 }