PWG2/SPECTRA/AliAnalysisTaskChargedHadronSpectra.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 //                                                                       //
  19 //                                                                       //
  20 // Analysis for identified particle spectra measured with TPC dE/dx.     //
  21 //                                                                       //
  22 //                                                                       //
  23 ///////////////////////////////////////////////////////////////////////////
  24
  25 #include "Riostream.h"
  26 #include "TChain.h"
  27 #include "TTree.h"
  28 #include "TH1F.h"
  29 #include "TH2F.h"
  30 #include "TH3F.h"
  31 #include "TList.h"
  32 #include "TMath.h"
  33 #include "TCanvas.h"
  34 #include "TObjArray.h"
  35 #include "TF1.h"
  36 #include "TFile.h"
  37
  38 #include "AliAnalysisTaskSE.h"
  39 #include "AliAnalysisManager.h"
  40
  41 #include "AliHeader.h"
  42 #include "AliGenPythiaEventHeader.h"
  43 #include "AliGenCocktailEventHeader.h"
  44
  45 #include "AliPID.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliESDVertex.h"
  48 #include "AliESDEvent.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliESDtrack.h"
  51 #include "AliTPCpidESD.h"
  52
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56
  57 #include "AliLog.h"
  58
  59 #include "AliAnalysisTaskChargedHadronSpectra.h"
  60
  61
  62 ClassImp(AliAnalysisTaskChargedHadronSpectra)
  63
  64 //________________________________________________________________________
  65 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra()
  66   : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fPidObject(0),
  67     fMCtrue(0),
  68     fAlephParameters(),
  69     fHistPtMCKaon(0),
  70     fHistPtMCProton(0),
  71     fHistPtMCPion(0),
  72     fHistPtMCElectron(0),
  73     fHistPtMCMuon(0),
  74     fHistPtEtaKaon(0),
  75     fHistPtEtaKaonNoKink(0),
  76     fHistPtEtaProton(0),
  77     fHistPtEtaProtonDCA(0),
  78     fHistPtEtaPion(0),
  79     fHistPtEtaElectron(0),
  80     fHistClassicKaon(0),
  81     fHistClassicProton(0),
  82     fHistClassicPion(0),
  83     fHistClassicElectron(0),
  84     fDeDx(0),
  85     fHistTrackPerEvent(0),
  86     fHistTrackPerEventMC(0),
  87     fSecProtons(0),
  88     fVertexZ(0),
  89     fHistEtaNcls(0),
  90     fHistEtaPhi(0),
  91     fHistEffProton(0),
  92     fHistEffProtonDCA(0),
  93     fHistEffPion(0),
  94     fHistEffKaon(0),
  95     fHighPtElectrons(0),
  96     fHighPtHadrons(0)
  97 {
  98   // default Constructor
  99 }
 100
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra(const char *name)
 104   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fPidObject(0),
 105     fMCtrue(0),
 106     fAlephParameters(),
 107     fHistPtMCKaon(0),
 108     fHistPtMCProton(0),
 109     fHistPtMCPion(0),
 110     fHistPtMCElectron(0),
 111     fHistPtMCMuon(0),
 112     fHistPtEtaKaon(0),
 113     fHistPtEtaKaonNoKink(0),
 114     fHistPtEtaProton(0),
 115     fHistPtEtaProtonDCA(0),
 116     fHistPtEtaPion(0),
 117     fHistPtEtaElectron(0),
 118     fHistClassicKaon(0),
 119     fHistClassicProton(0),
 120     fHistClassicPion(0),
 121     fHistClassicElectron(0),
 122     fDeDx(0),
 123     fHistTrackPerEvent(0),
 124     fHistTrackPerEventMC(0),
 125     fSecProtons(0),
 126     fVertexZ(0),
 127     fHistEtaNcls(0),
 128     fHistEtaPhi(0),
 129     fHistEffProton(0),
 130     fHistEffProtonDCA(0),
 131     fHistEffPion(0),
 132     fHistEffKaon(0),
 133     fHighPtElectrons(0),
 134     fHighPtHadrons(0)
 135 {
 136   //
 137   // standard constructur which should be used - PID objects is initialized
 138   //
 139
 140   fMCtrue = kTRUE; // change two things for processing real data: 1. set this to false! / 2. change ALEPH parameters
 141   fAlephParameters[0] = 4.23232575531564326e+00;//50*0.76176e-1;
 142   fAlephParameters[1] = 8.68482806165147636e+00;//10.632;
 143   fAlephParameters[2] = 1.34000000000000005e-05;//0.13279e-4;
 144   fAlephParameters[3] = 2.30445734159456084e+00;//1.8631;
 145   fAlephParameters[4] = 2.25624744086878559e+00;//1.9479;
 146
 147   fPidObject = new AliTPCpidESD();
 148   fPidObject->SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
 149
 150   // Constructor
 151   Printf("*** CONSTRUCTOR CALLED ****");
 152   // Define input and output slots here
 153   // Input slot #0 works with a TChain
 154   //DefineInput(0, TChain::Class()); <-> not needed in AliAnalysisTaskSE
 155
 156   // Output slot #0 writes into a TList container
 157   DefineOutput(1, TList::Class());
 158
 159
 160
 161 }
 162
 163
 164 //________________________________________________________________________
 165 void AliAnalysisTaskChargedHadronSpectra::UserCreateOutputObjects()
 166 {
 167   // Create histograms
 168   // Called once
 169   fListHist = new TList();
 170   //fListHist->SetOwner(); // Whoever knows how the data handling is ...?
 171
 172   const Int_t kPtBins = 2*56;
 173   const Double_t kPtMax  = 16.0;
 174   const Int_t kEtaBins = 4;
 175   const Double_t kEtaMax = 0.8;
 176   const Int_t kDeDxBins = 200;
 177   const Double_t kDeDxMax = 1;
 178   const Int_t kMultBins = 10;
 179   const Int_t kMultMax = 300;
 180
 181   // sort pT-bins ..
 182   Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
 183   Int_t indexes[kPtBins+1];
 184   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 185   Double_t binsPt[kPtBins+1];
 186   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 187
 188
 189   // MC histograms
 190   fHistPtMCKaon = new TH3F("HistPtMCKaon", "PtEtaKaon; mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 191   fHistPtMCProton = new TH3F("HistPtMCProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 192   fHistPtMCPion = new TH3F("HistPtMCPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 193   fHistPtMCElectron = new TH3F("HistPtMCElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 194   fHistPtMCMuon = new TH3F("HistPtMCMuon", "PtEtaMuon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 195   //
 196   fHistPtMCKaon->GetZaxis()->Set(kPtBins, binsPt);
 197   fHistPtMCProton->GetZaxis()->Set(kPtBins, binsPt);
 198   fHistPtMCPion->GetZaxis()->Set(kPtBins, binsPt);
 199   fHistPtMCElectron->GetZaxis()->Set(kPtBins, binsPt);
 200   fHistPtMCMuon->GetZaxis()->Set(kPtBins, binsPt);
 201   //
 202   fListHist->Add(fHistPtMCKaon);
 203   fListHist->Add(fHistPtMCProton);
 204   fListHist->Add(fHistPtMCPion);
 205   fListHist->Add(fHistPtMCElectron);
 206   fListHist->Add(fHistPtMCMuon);
 207
 208   // reconstructed particle histograms
 209   fHistPtEtaKaon = new TH3F("HistPtEtaKaon", "PtEtaKaon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 210   fHistPtEtaKaonNoKink = new TH3F("HistPtEtaKaonNoKink", "PtEtaKaonNoKink;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 211   fHistPtEtaProton = new TH3F("HistPtEtaProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 212   fHistPtEtaProtonDCA = new TH3F("HistPtEtaProtonDCA", "PtEtaProton;DCA (cm); #eta; p_{T} (GeV)",200,0,15,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 213   fHistPtEtaPion = new TH3F("HistPtEtaPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 214   fHistPtEtaElectron = new TH3F("HistPtEtaElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 215   //
 216   fHistPtEtaKaon->GetZaxis()->Set(kPtBins, binsPt);
 217   fHistPtEtaKaonNoKink->GetZaxis()->Set(kPtBins, binsPt);
 218   fHistPtEtaProton->GetZaxis()->Set(kPtBins, binsPt);
 219   fHistPtEtaProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
 220   fHistPtEtaPion->GetZaxis()->Set(kPtBins, binsPt);
 221   fHistPtEtaElectron->GetZaxis()->Set(kPtBins, binsPt);
 222   //
 223   fListHist->Add(fHistPtEtaKaon);
 224   fListHist->Add(fHistPtEtaKaonNoKink);
 225   fListHist->Add(fHistPtEtaProton);
 226   fListHist->Add(fHistPtEtaProtonDCA);
 227   fListHist->Add(fHistPtEtaPion);
 228   fListHist->Add(fHistPtEtaElectron);
 229
 230   // histograms for the classical analysis
 231   fHistClassicKaon = new TH3F("HistClassicKaon", "PtEtaKaon;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 232   fHistClassicProton = new TH3F("HistClassicProton", "PtEtaProton;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 233   fHistClassicPion = new TH3F("HistClassicPion", "PtEtaPion;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 234   fHistClassicElectron = new TH3F("HistClassicElectron", "PtEtaElectron;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 235   //
 236   fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
 237   fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
 238   fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
 239   fHistClassicElectron->GetXaxis()->Set(kPtBins, binsPt);
 240   //
 241   fListHist->Add(fHistClassicKaon);
 242   fListHist->Add(fHistClassicProton);
 243   fListHist->Add(fHistClassicPion);
 244   fListHist->Add(fHistClassicElectron);
 245
 246   // histograms of general interest
 247   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,20.,1000,0.,500);
 248   BinLogX(fDeDx);
 249   fListHist->Add(fDeDx);
 250   fHistTrackPerEvent = new TH1F("HistTrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",101,-0.5,100);
 251   fListHist->Add(fHistTrackPerEvent);
 252   fHistTrackPerEventMC  = new TH2F("HistTrackPerEventMC", "Tracks per event;Number of Tracks;Number of Events",10,-0.5,9,101,-0.5,100);
 253   fListHist->Add(fHistTrackPerEventMC);
 254   fSecProtons = new TH2F("SecProtons", "xy;x;y",100,-10,10,100,-10,10);
 255   fListHist->Add(fSecProtons);
 256   fVertexZ = new TH1F("VertexZ", "vertex position;z (cm);counts",200,-50,50);
 257   fListHist->Add(fVertexZ);
 258   //
 259   fHistEtaNcls = new TH2F("HistEtaNcls", "EtaNcls;#eta;Ncls",100,-1.5,1.5,80,0,160);
 260   fListHist->Add(fHistEtaNcls);
 261   fHistEtaPhi = new TH2F("HistEtaPhi", "EtaNcls;#eta;#phi",100,-1.5,1.5,100,0,2*TMath::Pi());
 262   fListHist->Add(fHistEtaPhi);
 263
 264   // histograms for a refined efficiency study
 265   fHistEffProton = new TH3F("HistEffProton", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 266   fListHist->Add(fHistEffProton);
 267   fHistEffProtonDCA  = new TH3F("HistEffProtonDCA", "Eff;p_{T} (GeV); code",10,-0.5,9.5,200,0,15,kPtBins,-kPtMax,kPtMax);
 268   fListHist->Add(fHistEffProtonDCA);
 269   fHistEffPion = new TH3F("HistEffPion", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 270   fListHist->Add(fHistEffPion);
 271   fHistEffKaon = new TH3F("HistEffKaon", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 272   fListHist->Add(fHistEffKaon);
 273   //
 274   fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
 275   fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
 276   fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
 277   fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
 278   //
 279
 280
 281   // histograms for dN/dpT
 282   fHighPtElectrons = new TH1F("HistHighPtElectrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
 283   fListHist->Add(fHighPtElectrons);
 284   fHighPtHadrons = new TH1F("HistHighPtHadrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
 285   fListHist->Add(fHighPtHadrons);
 286
 287   //
 288   // Create Objects which are needed for the rest of the analysis
 289   //
 290   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 291   fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2);  // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
 292   fESDtrackCuts->SetMaxNsigmaToVertex(3);
 293   fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
 294   fESDtrackCuts->SetRequireTPCRefit(kFALSE);
 295   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 296   fESDtrackCuts->SetMinNClustersTPC(100);
 297   fESDtrackCuts->SetMaxChi2PerClusterTPC(3.5);
 298   //fESDtrackCuts->SetMaxDCAToVertex(3);
 299   //
 300   fListHist->Add(fESDtrackCuts);
 301   //
 302   //
 303
 304 }
 305
 306 //________________________________________________________________________
 307 void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *)
 308 {
 309   // main event loop
 310
 311   AliLog::SetGlobalLogLevel(AliLog::kError);
 312   //
 313   // Check Monte Carlo information and other access first:
 314   //
 315   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 316   if (!eventHandler) {
 317     //Printf("ERROR: Could not retrieve MC event handler");
 318     if (fMCtrue) return;
 319   }
 320
 321   AliMCEvent* mcEvent = eventHandler->MCEvent();
 322   if (!mcEvent) {
 323     //Printf("ERROR: Could not retrieve MC event");
 324     if (fMCtrue) return;
 325   }
 326
 327   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
 328   if (!fESD) {
 329     //Printf("ERROR: fESD not available");
 330     return;
 331   }
 332
 333   if (!fESDtrackCuts) {
 334     Printf("ERROR: fESDtrackCuts not available");
 335     return;
 336   }
 337
 338   Int_t trackCounter = 0;
 339
 340  // monitor vertex position
 341  const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
 342  if (!vertex) {
 343    fHistTrackPerEvent->Fill(trackCounter);
 344    PostData(1, fListHist);
 345    return;
 346  } else {
 347    if (vertex->GetZv() != 0) fVertexZ->Fill(vertex->GetZv());
 348    if (TMath::Abs(vertex->GetZv()) > 10) {
 349      fHistTrackPerEvent->Fill(trackCounter);
 350      PostData(1, fListHist);
 351      return;
 352    }
 353  }
 354
 355
 356
 357  // 1st track loop to determine multiplicities
 358  for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 359    AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
 360    if (!track) continue;
 361    if (fESDtrackCuts->AcceptTrack(track) &&  TMath::Abs(track->Eta())< 0.9) {
 362      trackCounter++;
 363      //
 364    }
 365    delete track;
 366  }
 367
 368
 369  // MC loop
 370  AliStack* stack = mcEvent->Stack();
 371
 372  if (fMCtrue) {
 373    //
 374    //
 375    //
 376    AliHeader * header = mcEvent->Header();
 377    Int_t processtype = GetPythiaEventProcessType(header);
 378    // non diffractive
 379    if (processtype !=92 && processtype !=93 && processtype != 94) fHistTrackPerEventMC->Fill(1., trackCounter);
 380    // single diffractive
 381    if ((processtype == 92 || processtype == 93)) fHistTrackPerEventMC->Fill(2., trackCounter);
 382    // double diffractive
 383    if (processtype == 94) fHistTrackPerEventMC->Fill(3., trackCounter);
 384    //
 385    Int_t mult = trackCounter;//stack->GetNprimary();
 386    fHistTrackPerEventMC->Fill(0., trackCounter);
 387    //
 388    for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 389      TParticle * trackMC = stack->Particle(i);
 390      Int_t pdg = trackMC->GetPdgCode();
 391      //
 392      Double_t xv = trackMC->Vx();
 393      Double_t yv = trackMC->Vy();
 394      Double_t zv = trackMC->Vz();
 395      Double_t dxy = 0;
 396      dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
 397      Double_t dz = 0;
 398      dz = TMath::Abs(zv); // so stupid to avoid warnings
 399      //
 400      // vertex cut - selection of primaries
 401      //
 402      //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
 403      if (pdg == 11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select e- only; before vertex
 404      if (pdg == -11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select e+ only; before vertex
 405      //
 406      if (!stack->IsPhysicalPrimary(i)) continue;
 407      //
 408      if (pdg == 321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select K+ only
 409      if (pdg == 211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select Pi+ only
 410      if (pdg == 2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select p+ only
 411      if (pdg == -13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select mu+ only
 412      //
 413      if (pdg == -321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select K- only
 414      if (pdg == -211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select Pi- only
 415      if (pdg == -2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select p- only
 416      if (pdg == 13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select mu- only
 417      //
 418      // special Kaon checks - those which decay before entering the TPC fiducial volume
 419      //
 420      if (TMath::Abs(pdg)==321) {
 421        TParticle * trackDaughterMC = stack->Particle(TMath::Abs(trackMC->GetFirstDaughter()));
 422        Int_t pdgDaughter = trackDaughterMC->GetPdgCode();
 423        if (pdgDaughter == TMath::Abs(13) && pdg == 321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // Kaon control hist
 424        if (pdgDaughter == TMath::Abs(13) && pdg == -321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // Kaon control hist
 425      }
 426    }
 427  }
 428
 429  //end MC tracks loop
 430
 431
 432  // ESD loop
 433  // definition of PID functions --> to be put to CreateOutputObjects()
 434
 435  TF1 foProton("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 436  TF1 foPion("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 437  TF1 foElec("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 438  TF1 foKaon("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 439  TF1 foMuon("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 440
 441  foProton.SetParameters(fAlephParameters);
 442  foPion.SetParameters(fAlephParameters);
 443  foElec.SetParameters(fAlephParameters);
 444  foKaon.SetParameters(fAlephParameters);
 445  foMuon.SetParameters(fAlephParameters);
 446
 447  const Float_t k2sigmaCorr = 1/0.9545;
 448
 449  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 450
 451  for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 452    //
 453    // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
 454    //
 455    AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
 456    if (!track) continue;
 457    AliESDtrack *trackDummy = fESD->GetTrack(i);
 458    if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
 459      delete track;
 460      continue;
 461    }
 462    AliExternalTrackParam trackIn(*trackDummy->GetInnerParam());
 463    Double_t ptot = trackIn.GetP(); // momentum for dEdx determination
 464    Double_t pT = track->Pt();
 465    Double_t eta = TMath::Abs(track->Eta());
 466    //
 467    fHistEtaNcls->Fill(track->Eta(),track->GetTPCNcls());
 468    fHistEtaPhi->Fill(track->Eta(),track->Phi());
 469    // cut for dead regions in the detector
 470    //   if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
 471    //
 472    // 2.a) apply some standard track cuts according to general recommendations
 473    //
 474
 475    if (!fESDtrackCuts->AcceptTrack(track)) {
 476      //
 477      if (fMCtrue) {
 478        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 479        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 480        if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(6,eta,track->GetSign()*pT); // Kaon control hist
 481      }
 482      //
 483      delete track;
 484      continue;
 485    }
 486
 487    //
 488    // 3. make the PID
 489    //
 490    Double_t sign = track->GetSign();
 491
 492    Double_t tpcSignal = track->GetTPCsignal();
 493    //
 494    fDeDx->Fill(ptot, tpcSignal);
 495    //
 496    //
 497    //
 498    fHistClassicKaon->Fill(sign*pT,eta,(tpcSignal-foKaon.Eval(ptot))/foKaon.Eval(ptot));
 499    fHistClassicProton->Fill(sign*pT,eta,(tpcSignal-foProton.Eval(ptot))/foProton.Eval(ptot));
 500    fHistClassicPion->Fill(sign*pT,eta,(tpcSignal-foPion.Eval(ptot))/foPion.Eval(ptot));
 501    fHistClassicElectron->Fill(sign*pT,eta,(tpcSignal-foElec.Eval(ptot))/foElec.Eval(ptot));
 502    //
 503    // fill histograms for dNdpT
 504    //
 505    if (eta < 0.15) {
 506      fHighPtHadrons->Fill(pT);
 507      Double_t delta = (tpcSignal - foElec.Eval(ptot))/foElec.Eval(ptot);
 508      if (delta > 0) fHighPtElectrons->Fill(pT);
 509    }
 510    //
 511    /* 2sigma PID with 2sigma eff correction! */
 512    // PION
 513    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kPion)) < 2) {
 514      fHistPtEtaPion->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 515      //
 516      if (trackCounter < 300 && fMCtrue) {
 517        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 518        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 519        if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,eta,sign*pT,k2sigmaCorr);
 520        if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 521          fHistEffPion->Fill(1,eta,sign*pT,k2sigmaCorr);
 522          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 523          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 524          if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffPion->Fill(3,eta,sign*pT,k2sigmaCorr);
 525        }
 526        if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,eta,sign*pT,k2sigmaCorr);
 527        if (TMath::Abs(trackMC->GetPdgCode()) == 13)  fHistEffPion->Fill(4,eta,sign*pT,k2sigmaCorr);
 528        }
 529    }
 530    // KAON
 531    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kKaon)) < 2) {
 532      fHistPtEtaKaon->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 533      //
 534      if (trackCounter < 300 && fMCtrue) {
 535        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 536        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 537        if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,eta,sign*pT,k2sigmaCorr);
 538        if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 539          fHistEffKaon->Fill(1,eta,sign*pT,k2sigmaCorr);
 540          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 541          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 542          if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffKaon->Fill(3,eta,sign*pT,k2sigmaCorr);
 543        }
 544          if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,eta,sign*pT,k2sigmaCorr);
 545      }
 546    }
 547    // KAON NO KINK
 548    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kKaon)) < 2 && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 549    // PROTON
 550    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kProton)) < 2) {
 551      fHistPtEtaProton->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 552      //
 553      track->GetImpactParameters(dca,cov);
 554      Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
 555      fHistPtEtaProtonDCA->Fill(primVtxDCA, eta, sign*pT, k2sigmaCorr);
 556      //
 557      if (trackCounter < 300 && fMCtrue) {
 558        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 559        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 560        if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 561          fHistEffProton->Fill(0.,eta,sign*pT,k2sigmaCorr);
 562          if (eta < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
 563        }
 564        if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 565          fHistEffProton->Fill(1,eta,sign*pT,k2sigmaCorr);
 566          if (eta < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
 567          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 568          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 569          if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
 570            fHistEffProton->Fill(3,eta,sign*pT,k2sigmaCorr);
 571            if (eta < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
 572          }
 573          if (code!=3122&&code!=3222&&code!=3212&&code!=3112&&code!=3322&&code!=3312&&code!=3332&&code!=130&&code!=310) fSecProtons->Fill(trackMC->Vx(),trackMC->Vy());
 574        }
 575        if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,eta,sign*pT,k2sigmaCorr);
 576      }
 577    }
 578    // ELECTRON
 579    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kElectron))) fHistPtEtaElectron->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 580
 581    delete track;
 582
 583  } // end of track loop
 584
 585  fHistTrackPerEvent->Fill(trackCounter);
 586
 587  // Post output data
 588
 589
 590  PostData(1, fListHist);
 591
 592 }
 593
 594
 595 //________________________________________________________________________
 596 void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
 597 {
 598   // Draw result to the screen
 599   // Called once at the end of the query
 600   cout << "*** TERMINATE ***" << endl;
 601
 602 }
 603
 604
 605
 606
 607 //________________________________________________________________________
 608 void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
 609
 610   // Method for the correct logarithmic binning of histograms
 611
 612   TAxis *axis = h->GetXaxis();
 613   int bins = axis->GetNbins();
 614
 615   Double_t from = axis->GetXmin();
 616   Double_t to = axis->GetXmax();
 617   Double_t *newBins = new Double_t[bins + 1];
 618
 619   newBins[0] = from;
 620   Double_t factor = pow(to/from, 1./bins);
 621
 622   for (int i = 1; i <= bins; i++) {
 623    newBins[i] = factor * newBins[i-1];
 624   }
 625   axis->Set(bins, newBins);
 626   delete newBins;
 627
 628 }
 629
 630
 631 //________________________________________________________________________
 632 Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
 633   //
 634   // get the process type of the event.
 635   //
 636
 637   // can only read pythia headers, either directly or from cocktalil header
 638   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 639
 640   if (!pythiaGenHeader) {
 641
 642     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 643     if (!genCocktailHeader) {
 644       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 645       return -1;
 646     }
 647
 648     TList* headerList = genCocktailHeader->GetHeaders();
 649     if (!headerList) {
 650       return -1;
 651     }
 652
 653     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 654       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 655       if (pythiaGenHeader)
 656         break;
 657     }
 658
 659     if (!pythiaGenHeader) {
 660       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
 661       return -1;
 662     }
 663   }
 664
 665   if (adebug) {
 666     //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
 667   }
 668
 669   return pythiaGenHeader->ProcessType();
 670 }
 671
 672
 673 //________________________________________________________________________
 674 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBin, const Double_t * AlephParams) {
 675   //
 676   // The multiple Gauss fitting is happening here...
 677   // input histogram is of the form (Pt, eta, difference in dEdx)
 678   //
 679
 680   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 681   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 682   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 683   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 684   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 685   //
 686   foProton->SetParameters(AlephParams);
 687   foPion->SetParameters(AlephParams);
 688   foElec->SetParameters(AlephParams);
 689   foKaon->SetParameters(AlephParams);
 690   foMuon->SetParameters(AlephParams);
 691
 692
 693   const Double_t kSigma = 0.06; // expected resolution (roughly)
 694
 695   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 696   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 697   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 698
 699   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 700   // sort pT-bins ..
 701   Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
 702   Int_t indexes[kPtBins+1];
 703   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 704   Double_t binsPt[kPtBins+1];
 705   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 706   histPt->GetXaxis()->Set(kPtBins, binsPt);
 707   //
 708
 709   TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
 710   canvMany->Print("canvManyProton.ps[");
 711
 712   for(Int_t x=1; x < kPtBins+1; x++) {
 713     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 714     for(Int_t y=1; y < kEtaBins+1; y++) {
 715       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 716       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 717       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 718       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 719       histDeDx->SetTitle(Form("pt_%f",pT));
 720       Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
 721       //
 722       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
 723       Double_t parameters[12] = {maxYield/2.,0,kSigma,maxYield/2.,(foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma};
 724       gausFit->SetParameters(parameters);
 725       gausFit->SetParLimits(0,0.,maxYield);
 726       gausFit->SetParLimits(1,-0.05,0.05);
 727       gausFit->SetParLimits(2,0.04,0.08);
 728       //
 729       Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 730       gausFit->SetParLimits(0,0.,maxYield);
 731       gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
 732       gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
 733       //
 734       Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 735       gausFit->SetParLimits(0,0.,maxYield);
 736       gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
 737       gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
 738       //
 739       Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 740       gausFit->SetParLimits(0,0.,maxYield);
 741       gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
 742       gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
 743       histDeDx->Fit("gausFit","QNR");
 744       gausFit->GetParameters(parameters);
 745       // visualization
 746       if (y == EtaBin) {
 747         canvMany->cd(x);
 748         gPad->SetLogy();
 749         histDeDx->SetMarkerStyle(21);
 750         histDeDx->SetMarkerSize(0.7);
 751         histDeDx->Draw("E");
 752         gausFit->SetLineWidth(2);
 753         gausFit->Draw("same");
 754         //
 755         TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
 756         TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
 757         TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
 758         TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
 759         gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
 760         gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
 761         gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
 762         gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
 763         //
 764         gausFit0.SetLineColor(4);
 765         gausFit1.SetLineColor(2);
 766         gausFit2.SetLineColor(6);
 767         gausFit3.SetLineColor(8);
 768         //
 769         gausFit0.SetLineWidth(1);
 770         gausFit1.SetLineWidth(1);
 771         gausFit2.SetLineWidth(1);
 772         gausFit3.SetLineWidth(1);
 773         //
 774         gausFit0.Draw("same");
 775         gausFit1.Draw("same");
 776         gausFit2.Draw("same");
 777         gausFit3.Draw("same");
 778         canvMany->Print("canvManyProton.ps");
 779       }
 780
 781       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
 782       delete gausFit;
 783       //
 784       // stupid solution --> remove as soon as possible
 785       /*yield = 0;
 786       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 787       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 788       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 789         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 790         yield += gausYield->Eval(delta);
 791       }*/
 792       //
 793       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 794       //delete gausYield;
 795     }
 796   }
 797   canvMany->Print("canvManyProton.ps]");
 798
 799   TCanvas * canvPt = new TCanvas();
 800   canvPt->cd();
 801   histPt->Draw();
 802
 803   return histPt;
 804 }
 805
 806
 807
 808 //________________________________________________________________________
 809 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
 810   //
 811   // The multiple Gauss fitting is happening here...
 812   //
 813
 814   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 815   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 816   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 817   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 818   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 819   //
 820   foProton->SetParameters(AlephParams);
 821   foPion->SetParameters(AlephParams);
 822   foElec->SetParameters(AlephParams);
 823   foKaon->SetParameters(AlephParams);
 824   foMuon->SetParameters(AlephParams);
 825
 826   // input histogram is of the form (Pt, eta, difference in dEdx)
 827
 828   const Double_t kSigma = 0.06; // expected resolution (roughly)
 829
 830   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 831   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 832   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 833
 834   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 835   // sort pT-bins ..
 836   Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
 837   Int_t indexes[kPtBins+1];
 838   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 839   Double_t binsPt[kPtBins+1];
 840   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 841   histPt->GetXaxis()->Set(kPtBins, binsPt);
 842   //
 843
 844   TCanvas * canvMany = new TCanvas();
 845   canvMany->Divide(8,5);
 846
 847   for(Int_t x=1; x < kPtBins+1; x++) {
 848     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 849     for(Int_t y=1; y < kEtaBins+1; y++) {
 850       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 851       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 852       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 853       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 854       //
 855       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
 856       //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foKaon->Eval(ptot)/foPion->Eval(ptot)),kSigma};
 857       Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma};
 858       gausFit->SetParameters(parameters);
 859       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
 860       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
 861       //
 862       gausFit->FixParameter(4,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 863       gausFit->FixParameter(5,kSigma);
 864       gausFit->FixParameter(7,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 865       gausFit->FixParameter(8,kSigma);
 866       gausFit->FixParameter(10,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 867       gausFit->FixParameter(11,kSigma);
 868       histDeDx->Fit("gausFit","Q0");
 869       if (y == EtaBin && pT < 0) {
 870         canvMany->cd(x);
 871         gPad->SetLogy();
 872         histDeDx->Draw();
 873         gausFit->Draw("same");
 874       }
 875       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
 876       // stupid solution --> remove as soon as possible
 877       yield = 0;
 878       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 879       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 880       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 881         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 882         yield += gausYield->Eval(delta);
 883       }
 884       //
 885       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 886       //delete gausFit;
 887       delete gausYield;
 888     }
 889   }
 890
 891   TCanvas * canvPt = new TCanvas();
 892   canvPt->cd();
 893   histPt->Draw();
 894
 895   return histPt;
 896
 897 }
 898
 899
 900 //________________________________________________________________________
 901 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
 902   //
 903   // The multiple Gauss fitting is happening here...
 904   //
 905
 906   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 907   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 908   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 909   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 910   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 911   //
 912   foProton->SetParameters(AlephParams);
 913   foPion->SetParameters(AlephParams);
 914   foElec->SetParameters(AlephParams);
 915   foKaon->SetParameters(AlephParams);
 916   foMuon->SetParameters(AlephParams);
 917
 918   const Double_t kSigma = 0.06; // expected resolution (roughly)
 919
 920   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 921   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 922   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 923
 924   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 925   // sort pT-bins ..
 926   Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
 927   Int_t indexes[kPtBins+1];
 928   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 929   Double_t binsPt[kPtBins+1];
 930   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 931   histPt->GetXaxis()->Set(kPtBins, binsPt);
 932   //
 933
 934   TCanvas * canvMany = new TCanvas();
 935   canvMany->Divide(8,5);
 936
 937   for(Int_t x=1; x < kPtBins+1; x++) {
 938     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 939     for(Int_t y=1; y < kEtaBins+1; y++) {
 940       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 941       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 942       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 943       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 944       //
 945       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
 946       //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foPion->Eval(ptot)/foPion->Eval(ptot)),kSigma};
 947       Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foPion->Eval(ptot)-foKaon->Eval(ptot))/foPion->Eval(ptot),kSigma};
 948       gausFit->SetParameters(parameters);
 949       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
 950       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
 951       //
 952       gausFit->FixParameter(4,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
 953       gausFit->FixParameter(5,kSigma);
 954       gausFit->FixParameter(7,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
 955       gausFit->FixParameter(8,kSigma);
 956       gausFit->FixParameter(10,(foPion->Eval(ptot)- foKaon->Eval(ptot))/foKaon->Eval(ptot));
 957       gausFit->FixParameter(11,kSigma);
 958       histDeDx->Fit("gausFit","Q0");
 959       if (y == EtaBin && pT < 0) {
 960         canvMany->cd(x);
 961         gPad->SetLogy();
 962         histDeDx->Draw();
 963         gausFit->Draw("same");
 964       }
 965       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
 966       // stupid solution --> remove as soon as possible
 967       yield = 0;
 968       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 969       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 970       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 971         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 972         yield += gausYield->Eval(delta);
 973       }
 974       //
 975       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 976       //delete gausFit;
 977       delete gausYield;
 978     }
 979   }
 980
 981   TCanvas * canvPt = new TCanvas();
 982   canvPt->cd();
 983   histPt->Draw();
 984
 985   return histPt;
 986
 987 }
 988
 989
 990 //________________________________________________________________________
 991 void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const  TList * ListOfHistogramsData, const Char_t *filename) {
 992
 993   // Define ranges
 994   const Int_t kMultLow  = 1;
 995   const Int_t kMultHigh = 10;
 996
 997   const Int_t kEtaHigh  = 4;
 998
 999   // Extract histograms for the MC efficiency study
1000   TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1001   TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1002   TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
1003
1004   TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
1005   TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1006   TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
1007
1008   TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1009   TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1010   TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1011
1012   // Extract data for the final spectra
1013   TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
1014   TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1015   TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
1016
1017   // "MC" of the real data for consistency checks
1018   TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1019   TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1020   TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
1021
1022   TFile spectraFile(filename,"RECREATE");
1023
1024   // 1. Protons
1025   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1026     TH1D *protonRaw     = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1027     protonRaw->Sumw2();
1028     TH1D *protonMC      = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1029     protonMC->Sumw2();
1030     TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1031     protonRawPrim->Sumw2();
1032     TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1033     protonSecReac->Sumw2();
1034     TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1035     protonSecWeak->Sumw2();
1036     TH1D *protonMis     = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1037     protonMis->Sumw2();
1038     protonSecReac->Add(protonSecWeak,-1);
1039     //
1040     // 1. total (physical) efficiency
1041     //
1042     TH1D *protonEffTot = new TH1D(*protonRaw);
1043     protonEffTot->Divide(protonRaw,protonMC);
1044     protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1045     protonEffTot->Write();
1046     //
1047     // 2. contributions from secondaries created by interactions
1048     //
1049     TH1D *protonEffSecReac = new TH1D(*protonRaw);
1050     TH1D *dummy = new TH1D(*protonRaw);
1051     dummy->Add(protonSecReac,-1);
1052     protonEffSecReac->Divide(protonRaw,dummy);
1053     protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1054     protonEffSecReac->Write();
1055     delete dummy;
1056     //
1057     // 3. contributions from secondaries from weak decays
1058     //
1059     TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1060     TH1D *dummy2 = new TH1D(*protonRaw);
1061     dummy2->Add(protonSecWeak,-1);
1062     protonEffSecWeak->Divide(protonRaw,dummy2);
1063     protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1064     protonEffSecWeak->Write();
1065     delete dummy2;
1066     //
1067     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1068     //
1069     TH1D *protonEffAbs = new TH1D(*protonRaw);
1070     protonEffAbs->Divide(protonRawPrim,protonMC);
1071     protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1072     protonEffAbs->Write();
1073     //
1074     // 5. contamination
1075     //
1076     TH1D *protonEffMis = new TH1D(*protonMis);
1077     protonEffMis->Divide(protonMis,protonRaw);
1078     protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1079     protonEffMis->Write();
1080
1081     //
1082     // PROCESS THE REAL DATA FROM HERE ON
1083     //
1084     TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1085     histPtProtonEta->Write();
1086     //
1087     // Apply corrections as you would like them to have...
1088     //
1089     TH1D *protonFinal = new TH1D(*histPtProtonEta);
1090     protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1091     protonFinal->Sumw2();
1092     protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1093     protonFinal->Write();
1094     //
1095     // if available save also the MC truth for consistency checks
1096     //
1097     TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1098     histPtProtonEtaMcData->Write();
1099   }
1100
1101   // 2. Kaons
1102   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1103     TH1D *kaonRaw     = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1104     kaonRaw->Sumw2();
1105     TH1D *kaonMC      = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1106     kaonMC->Sumw2();
1107     TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1108     kaonRawPrim->Sumw2();
1109     TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1110     kaonSecReac->Sumw2();
1111     TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1112     kaonSecWeak->Sumw2();
1113     TH1D *kaonMis     = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1114     kaonMis->Sumw2();
1115     kaonSecReac->Add(kaonSecWeak,-1);
1116     //
1117     // 1. total (physical) efficiency
1118     //
1119     TH1D *kaonEffTot = new TH1D(*kaonRaw);
1120     kaonEffTot->Divide(kaonRaw,kaonMC);
1121     kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1122     kaonEffTot->Write();
1123     //
1124     // 2. contributions from secondaries created by interactions
1125     //
1126     TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1127     TH1D *dummy = new TH1D(*kaonRaw);
1128     dummy->Add(kaonSecReac,-1);
1129     kaonEffSecReac->Divide(kaonRaw,dummy);
1130     kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1131     kaonEffSecReac->Write();
1132     delete dummy;
1133     //
1134     // 3. contributions from secondaries from weak decays
1135     //
1136     TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1137     TH1D *dummy2 = new TH1D(*kaonRaw);
1138     dummy2->Add(kaonSecWeak,-1);
1139     kaonEffSecWeak->Divide(kaonRaw,dummy2);
1140     kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1141     kaonEffSecWeak->Write();
1142     delete dummy2;
1143     //
1144     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1145     //
1146     TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1147     kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1148     kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1149     kaonEffAbs->Write();
1150     //
1151     // 5. contamination
1152     //
1153     TH1D *kaonEffMis = new TH1D(*kaonMis);
1154     kaonEffMis->Divide(kaonMis,kaonRaw);
1155     kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1156     kaonEffMis->Write();
1157
1158     //
1159     // PROCESS THE REAL DATA FROM HERE ON
1160     //
1161     TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1162     histPtKaonEta->Write();
1163     //
1164     // Apply corrections as you would like them to have...
1165     //
1166     TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1167     kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1168     kaonFinal->Sumw2();
1169     kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1170     kaonFinal->Write();
1171     //
1172     // if available save also the MC truth for consistency checks
1173     //
1174     TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1175     histPtKaonEtaMcData->Write();
1176   }
1177
1178
1179   // 3. Pions
1180   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1181     TH1D *pionRaw     = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1182     pionRaw->Sumw2();
1183     TH1D *pionMC      = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1184     pionMC->Sumw2();
1185     TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1186     pionRawPrim->Sumw2();
1187     TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1188     pionSecReac->Sumw2();
1189     TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1190     pionSecWeak->Sumw2();
1191     TH1D *pionMis     = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1192     pionMis->Sumw2();
1193     pionSecReac->Add(pionSecWeak,-1);
1194     //
1195     // 1. total (physical) efficiency
1196     //
1197     TH1D *pionEffTot = new TH1D(*pionRaw);
1198     pionEffTot->Divide(pionRaw,pionMC);
1199     pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1200     pionEffTot->Write();
1201     //
1202     // 2. contributions from secondaries created by interactions
1203     //
1204     TH1D *pionEffSecReac = new TH1D(*pionRaw);
1205     TH1D *dummy = new TH1D(*pionRaw);
1206     dummy->Add(pionSecReac,-1);
1207     pionEffSecReac->Divide(pionRaw,dummy);
1208     pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1209     pionEffSecReac->Write();
1210     delete dummy;
1211     //
1212     // 3. contributions from secondaries from weak decays
1213     //
1214     TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1215     TH1D *dummy2 = new TH1D(*pionRaw);
1216     dummy2->Add(pionSecWeak,-1);
1217     pionEffSecWeak->Divide(pionRaw,dummy2);
1218     pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1219     pionEffSecWeak->Write();
1220     delete dummy2;
1221     //
1222     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1223     //
1224     TH1D *pionEffAbs = new TH1D(*pionRaw);
1225     pionEffAbs->Divide(pionRawPrim,pionMC);
1226     pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1227     pionEffAbs->Write();
1228     //
1229     // 5. contamination
1230     //
1231     TH1D *pionEffMis = new TH1D(*pionMis);
1232     pionEffMis->Divide(pionMis,pionRaw);
1233     pionEffMis->SetName(Form("PionContam_%i",iEta));
1234     pionEffMis->Write();
1235
1236     //
1237     // PROCESS THE REAL DATA FROM HERE ON
1238     //
1239     TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1240     histPtPionEta->Write();
1241     //
1242     // Apply corrections as you would like them to have...
1243     //
1244     TH1D *pionFinal = new TH1D(*histPtPionEta);
1245     pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1246     pionFinal->Sumw2();
1247     pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1248     pionFinal->Write();
1249     //
1250     // if available save also the MC truth for consistency checks
1251     //
1252     TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1253     histPtPionEtaMcData->Write();
1254
1255   }
1256
1257   //
1258   // Close the file after everthing is done
1259   //
1260   spectraFile.Close();
1261
1262 }
1263
1264
1265
1266