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