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 "AliESDpid.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),fESDpid(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),fESDpid(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   fESDpid = new AliESDpid();
 148   fESDpid->GetTPCResponse().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    Double_t tpcMom = track->GetP();
 513    const AliExternalTrackParam *in = track->GetTPCInnerParam();
 514    if (in)
 515       tpcMom = in->GetP();
 516    // PION
 517    if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kPion)) < 2) {
 518      fHistPtEtaPion->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 519      //
 520      if (trackCounter < 300 && fMCtrue) {
 521        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 522        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 523        if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,eta,sign*pT,k2sigmaCorr);
 524        if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 525          fHistEffPion->Fill(1,eta,sign*pT,k2sigmaCorr);
 526          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 527          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 528          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);
 529        }
 530        if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,eta,sign*pT,k2sigmaCorr);
 531        if (TMath::Abs(trackMC->GetPdgCode()) == 13)  fHistEffPion->Fill(4,eta,sign*pT,k2sigmaCorr);
 532        }
 533    }
 534    // KAON
 535    if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2) {
 536      fHistPtEtaKaon->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 537      //
 538      if (trackCounter < 300 && fMCtrue) {
 539        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 540        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 541        if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,eta,sign*pT,k2sigmaCorr);
 542        if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 543          fHistEffKaon->Fill(1,eta,sign*pT,k2sigmaCorr);
 544          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 545          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 546          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);
 547        }
 548          if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,eta,sign*pT,k2sigmaCorr);
 549      }
 550    }
 551    // KAON NO KINK
 552    if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2 && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 553    // PROTON
 554    if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)) < 2) {
 555      fHistPtEtaProton->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 556      //
 557      track->GetImpactParameters(dca,cov);
 558      Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
 559      fHistPtEtaProtonDCA->Fill(primVtxDCA, eta, sign*pT, k2sigmaCorr);
 560      //
 561      if (trackCounter < 300 && fMCtrue) {
 562        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 563        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 564        if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 565          fHistEffProton->Fill(0.,eta,sign*pT,k2sigmaCorr);
 566          if (eta < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
 567        }
 568        if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 569          fHistEffProton->Fill(1,eta,sign*pT,k2sigmaCorr);
 570          if (eta < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
 571          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 572          Int_t code = TMath::Abs(trackMother->GetPdgCode());
 573          if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
 574            fHistEffProton->Fill(3,eta,sign*pT,k2sigmaCorr);
 575            if (eta < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
 576          }
 577          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());
 578        }
 579        if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,eta,sign*pT,k2sigmaCorr);
 580      }
 581    }
 582    // ELECTRON
 583    if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron))) fHistPtEtaElectron->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
 584
 585    delete track;
 586
 587  } // end of track loop
 588
 589  fHistTrackPerEvent->Fill(trackCounter);
 590
 591  // Post output data
 592
 593
 594  PostData(1, fListHist);
 595
 596 }
 597
 598
 599 //________________________________________________________________________
 600 void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
 601 {
 602   // Draw result to the screen
 603   // Called once at the end of the query
 604   cout << "*** TERMINATE ***" << endl;
 605
 606 }
 607
 608
 609
 610
 611 //________________________________________________________________________
 612 void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
 613
 614   // Method for the correct logarithmic binning of histograms
 615
 616   TAxis *axis = h->GetXaxis();
 617   int bins = axis->GetNbins();
 618
 619   Double_t from = axis->GetXmin();
 620   Double_t to = axis->GetXmax();
 621   Double_t *newBins = new Double_t[bins + 1];
 622
 623   newBins[0] = from;
 624   Double_t factor = pow(to/from, 1./bins);
 625
 626   for (int i = 1; i <= bins; i++) {
 627    newBins[i] = factor * newBins[i-1];
 628   }
 629   axis->Set(bins, newBins);
 630   delete newBins;
 631
 632 }
 633
 634
 635 //________________________________________________________________________
 636 Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
 637   //
 638   // get the process type of the event.
 639   //
 640
 641   // can only read pythia headers, either directly or from cocktalil header
 642   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 643
 644   if (!pythiaGenHeader) {
 645
 646     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 647     if (!genCocktailHeader) {
 648       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 649       return -1;
 650     }
 651
 652     TList* headerList = genCocktailHeader->GetHeaders();
 653     if (!headerList) {
 654       return -1;
 655     }
 656
 657     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 658       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 659       if (pythiaGenHeader)
 660         break;
 661     }
 662
 663     if (!pythiaGenHeader) {
 664       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
 665       return -1;
 666     }
 667   }
 668
 669   if (adebug) {
 670     //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
 671   }
 672
 673   return pythiaGenHeader->ProcessType();
 674 }
 675
 676
 677 //________________________________________________________________________
 678 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBin, const Double_t * AlephParams) {
 679   //
 680   // The multiple Gauss fitting is happening here...
 681   // input histogram is of the form (Pt, eta, difference in dEdx)
 682   //
 683
 684   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 685   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 686   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 687   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 688   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 689   //
 690   foProton->SetParameters(AlephParams);
 691   foPion->SetParameters(AlephParams);
 692   foElec->SetParameters(AlephParams);
 693   foKaon->SetParameters(AlephParams);
 694   foMuon->SetParameters(AlephParams);
 695
 696
 697   const Double_t kSigma = 0.06; // expected resolution (roughly)
 698
 699   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 700   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 701   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 702
 703   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 704   // sort pT-bins ..
 705   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};
 706   Int_t indexes[kPtBins+1];
 707   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 708   Double_t binsPt[kPtBins+1];
 709   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 710   histPt->GetXaxis()->Set(kPtBins, binsPt);
 711   //
 712
 713   TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
 714   canvMany->Print("canvManyProton.ps[");
 715
 716   for(Int_t x=1; x < kPtBins+1; x++) {
 717     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 718     for(Int_t y=1; y < kEtaBins+1; y++) {
 719       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 720       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 721       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 722       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 723       histDeDx->SetTitle(Form("pt_%f",pT));
 724       Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
 725       //
 726       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
 727       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};
 728       gausFit->SetParameters(parameters);
 729       gausFit->SetParLimits(0,0.,maxYield);
 730       gausFit->SetParLimits(1,-0.05,0.05);
 731       gausFit->SetParLimits(2,0.04,0.08);
 732       //
 733       Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 734       gausFit->SetParLimits(0,0.,maxYield);
 735       gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
 736       gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
 737       //
 738       Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 739       gausFit->SetParLimits(0,0.,maxYield);
 740       gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
 741       gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
 742       //
 743       Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
 744       gausFit->SetParLimits(0,0.,maxYield);
 745       gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
 746       gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
 747       histDeDx->Fit("gausFit","QNR");
 748       gausFit->GetParameters(parameters);
 749       // visualization
 750       if (y == EtaBin) {
 751         canvMany->cd(x);
 752         gPad->SetLogy();
 753         histDeDx->SetMarkerStyle(21);
 754         histDeDx->SetMarkerSize(0.7);
 755         histDeDx->Draw("E");
 756         gausFit->SetLineWidth(2);
 757         gausFit->Draw("same");
 758         //
 759         TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
 760         TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
 761         TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
 762         TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
 763         gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
 764         gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
 765         gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
 766         gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
 767         //
 768         gausFit0.SetLineColor(4);
 769         gausFit1.SetLineColor(2);
 770         gausFit2.SetLineColor(6);
 771         gausFit3.SetLineColor(8);
 772         //
 773         gausFit0.SetLineWidth(1);
 774         gausFit1.SetLineWidth(1);
 775         gausFit2.SetLineWidth(1);
 776         gausFit3.SetLineWidth(1);
 777         //
 778         gausFit0.Draw("same");
 779         gausFit1.Draw("same");
 780         gausFit2.Draw("same");
 781         gausFit3.Draw("same");
 782         canvMany->Print("canvManyProton.ps");
 783       }
 784
 785       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
 786       delete gausFit;
 787       //
 788       // stupid solution --> remove as soon as possible
 789       /*yield = 0;
 790       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 791       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 792       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 793         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 794         yield += gausYield->Eval(delta);
 795       }*/
 796       //
 797       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 798       //delete gausYield;
 799     }
 800   }
 801   canvMany->Print("canvManyProton.ps]");
 802
 803   TCanvas * canvPt = new TCanvas();
 804   canvPt->cd();
 805   histPt->Draw();
 806
 807   return histPt;
 808 }
 809
 810
 811
 812 //________________________________________________________________________
 813 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
 814   //
 815   // The multiple Gauss fitting is happening here...
 816   //
 817
 818   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 819   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 820   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 821   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 822   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 823   //
 824   foProton->SetParameters(AlephParams);
 825   foPion->SetParameters(AlephParams);
 826   foElec->SetParameters(AlephParams);
 827   foKaon->SetParameters(AlephParams);
 828   foMuon->SetParameters(AlephParams);
 829
 830   // input histogram is of the form (Pt, eta, difference in dEdx)
 831
 832   const Double_t kSigma = 0.06; // expected resolution (roughly)
 833
 834   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 835   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 836   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 837
 838   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 839   // sort pT-bins ..
 840   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};
 841   Int_t indexes[kPtBins+1];
 842   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 843   Double_t binsPt[kPtBins+1];
 844   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 845   histPt->GetXaxis()->Set(kPtBins, binsPt);
 846   //
 847
 848   TCanvas * canvMany = new TCanvas();
 849   canvMany->Divide(8,5);
 850
 851   for(Int_t x=1; x < kPtBins+1; x++) {
 852     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 853     for(Int_t y=1; y < kEtaBins+1; y++) {
 854       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 855       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 856       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 857       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 858       //
 859       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
 860       //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};
 861       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};
 862       gausFit->SetParameters(parameters);
 863       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
 864       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
 865       //
 866       gausFit->FixParameter(4,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 867       gausFit->FixParameter(5,kSigma);
 868       gausFit->FixParameter(7,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 869       gausFit->FixParameter(8,kSigma);
 870       gausFit->FixParameter(10,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
 871       gausFit->FixParameter(11,kSigma);
 872       histDeDx->Fit("gausFit","Q0");
 873       if (y == EtaBin && pT < 0) {
 874         canvMany->cd(x);
 875         gPad->SetLogy();
 876         histDeDx->Draw();
 877         gausFit->Draw("same");
 878       }
 879       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
 880       // stupid solution --> remove as soon as possible
 881       yield = 0;
 882       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 883       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 884       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 885         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 886         yield += gausYield->Eval(delta);
 887       }
 888       //
 889       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 890       //delete gausFit;
 891       delete gausYield;
 892     }
 893   }
 894
 895   TCanvas * canvPt = new TCanvas();
 896   canvPt->cd();
 897   histPt->Draw();
 898
 899   return histPt;
 900
 901 }
 902
 903
 904 //________________________________________________________________________
 905 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
 906   //
 907   // The multiple Gauss fitting is happening here...
 908   //
 909
 910   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 911   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 912   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 913   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 914   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 915   //
 916   foProton->SetParameters(AlephParams);
 917   foPion->SetParameters(AlephParams);
 918   foElec->SetParameters(AlephParams);
 919   foKaon->SetParameters(AlephParams);
 920   foMuon->SetParameters(AlephParams);
 921
 922   const Double_t kSigma = 0.06; // expected resolution (roughly)
 923
 924   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
 925   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 926   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 927
 928   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 929   // sort pT-bins ..
 930   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};
 931   Int_t indexes[kPtBins+1];
 932   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 933   Double_t binsPt[kPtBins+1];
 934   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 935   histPt->GetXaxis()->Set(kPtBins, binsPt);
 936   //
 937
 938   TCanvas * canvMany = new TCanvas();
 939   canvMany->Divide(8,5);
 940
 941   for(Int_t x=1; x < kPtBins+1; x++) {
 942     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 943     for(Int_t y=1; y < kEtaBins+1; y++) {
 944       Double_t eta = input->GetYaxis()->GetBinCenter(y);
 945       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 946       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 947       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
 948       //
 949       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
 950       //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};
 951       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};
 952       gausFit->SetParameters(parameters);
 953       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
 954       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
 955       //
 956       gausFit->FixParameter(4,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
 957       gausFit->FixParameter(5,kSigma);
 958       gausFit->FixParameter(7,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
 959       gausFit->FixParameter(8,kSigma);
 960       gausFit->FixParameter(10,(foPion->Eval(ptot)- foKaon->Eval(ptot))/foKaon->Eval(ptot));
 961       gausFit->FixParameter(11,kSigma);
 962       histDeDx->Fit("gausFit","Q0");
 963       if (y == EtaBin && pT < 0) {
 964         canvMany->cd(x);
 965         gPad->SetLogy();
 966         histDeDx->Draw();
 967         gausFit->Draw("same");
 968       }
 969       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
 970       // stupid solution --> remove as soon as possible
 971       yield = 0;
 972       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
 973       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
 974       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
 975         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
 976         yield += gausYield->Eval(delta);
 977       }
 978       //
 979       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
 980       //delete gausFit;
 981       delete gausYield;
 982     }
 983   }
 984
 985   TCanvas * canvPt = new TCanvas();
 986   canvPt->cd();
 987   histPt->Draw();
 988
 989   return histPt;
 990
 991 }
 992
 993
 994 //________________________________________________________________________
 995 void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const  TList * ListOfHistogramsData, const Char_t *filename) {
 996
 997   // Define ranges
 998   const Int_t kMultLow  = 1;
 999   const Int_t kMultHigh = 10;
1000
1001   const Int_t kEtaHigh  = 4;
1002
1003   // Extract histograms for the MC efficiency study
1004   TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1005   TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1006   TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
1007
1008   TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
1009   TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1010   TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
1011
1012   TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1013   TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1014   TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1015
1016   // Extract data for the final spectra
1017   TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
1018   TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1019   TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
1020
1021   // "MC" of the real data for consistency checks
1022   TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1023   TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1024   TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
1025
1026   TFile spectraFile(filename,"RECREATE");
1027
1028   // 1. Protons
1029   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1030     TH1D *protonRaw     = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1031     protonRaw->Sumw2();
1032     TH1D *protonMC      = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1033     protonMC->Sumw2();
1034     TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1035     protonRawPrim->Sumw2();
1036     TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1037     protonSecReac->Sumw2();
1038     TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1039     protonSecWeak->Sumw2();
1040     TH1D *protonMis     = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1041     protonMis->Sumw2();
1042     protonSecReac->Add(protonSecWeak,-1);
1043     //
1044     // 1. total (physical) efficiency
1045     //
1046     TH1D *protonEffTot = new TH1D(*protonRaw);
1047     protonEffTot->Divide(protonRaw,protonMC);
1048     protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1049     protonEffTot->Write();
1050     //
1051     // 2. contributions from secondaries created by interactions
1052     //
1053     TH1D *protonEffSecReac = new TH1D(*protonRaw);
1054     TH1D *dummy = new TH1D(*protonRaw);
1055     dummy->Add(protonSecReac,-1);
1056     protonEffSecReac->Divide(protonRaw,dummy);
1057     protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1058     protonEffSecReac->Write();
1059     delete dummy;
1060     //
1061     // 3. contributions from secondaries from weak decays
1062     //
1063     TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1064     TH1D *dummy2 = new TH1D(*protonRaw);
1065     dummy2->Add(protonSecWeak,-1);
1066     protonEffSecWeak->Divide(protonRaw,dummy2);
1067     protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1068     protonEffSecWeak->Write();
1069     delete dummy2;
1070     //
1071     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1072     //
1073     TH1D *protonEffAbs = new TH1D(*protonRaw);
1074     protonEffAbs->Divide(protonRawPrim,protonMC);
1075     protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1076     protonEffAbs->Write();
1077     //
1078     // 5. contamination
1079     //
1080     TH1D *protonEffMis = new TH1D(*protonMis);
1081     protonEffMis->Divide(protonMis,protonRaw);
1082     protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1083     protonEffMis->Write();
1084
1085     //
1086     // PROCESS THE REAL DATA FROM HERE ON
1087     //
1088     TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1089     histPtProtonEta->Write();
1090     //
1091     // Apply corrections as you would like them to have...
1092     //
1093     TH1D *protonFinal = new TH1D(*histPtProtonEta);
1094     protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1095     protonFinal->Sumw2();
1096     protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1097     protonFinal->Write();
1098     //
1099     // if available save also the MC truth for consistency checks
1100     //
1101     TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1102     histPtProtonEtaMcData->Write();
1103   }
1104
1105   // 2. Kaons
1106   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1107     TH1D *kaonRaw     = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1108     kaonRaw->Sumw2();
1109     TH1D *kaonMC      = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1110     kaonMC->Sumw2();
1111     TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1112     kaonRawPrim->Sumw2();
1113     TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1114     kaonSecReac->Sumw2();
1115     TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1116     kaonSecWeak->Sumw2();
1117     TH1D *kaonMis     = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1118     kaonMis->Sumw2();
1119     kaonSecReac->Add(kaonSecWeak,-1);
1120     //
1121     // 1. total (physical) efficiency
1122     //
1123     TH1D *kaonEffTot = new TH1D(*kaonRaw);
1124     kaonEffTot->Divide(kaonRaw,kaonMC);
1125     kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1126     kaonEffTot->Write();
1127     //
1128     // 2. contributions from secondaries created by interactions
1129     //
1130     TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1131     TH1D *dummy = new TH1D(*kaonRaw);
1132     dummy->Add(kaonSecReac,-1);
1133     kaonEffSecReac->Divide(kaonRaw,dummy);
1134     kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1135     kaonEffSecReac->Write();
1136     delete dummy;
1137     //
1138     // 3. contributions from secondaries from weak decays
1139     //
1140     TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1141     TH1D *dummy2 = new TH1D(*kaonRaw);
1142     dummy2->Add(kaonSecWeak,-1);
1143     kaonEffSecWeak->Divide(kaonRaw,dummy2);
1144     kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1145     kaonEffSecWeak->Write();
1146     delete dummy2;
1147     //
1148     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1149     //
1150     TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1151     kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1152     kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1153     kaonEffAbs->Write();
1154     //
1155     // 5. contamination
1156     //
1157     TH1D *kaonEffMis = new TH1D(*kaonMis);
1158     kaonEffMis->Divide(kaonMis,kaonRaw);
1159     kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1160     kaonEffMis->Write();
1161
1162     //
1163     // PROCESS THE REAL DATA FROM HERE ON
1164     //
1165     TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1166     histPtKaonEta->Write();
1167     //
1168     // Apply corrections as you would like them to have...
1169     //
1170     TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1171     kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1172     kaonFinal->Sumw2();
1173     kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1174     kaonFinal->Write();
1175     //
1176     // if available save also the MC truth for consistency checks
1177     //
1178     TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1179     histPtKaonEtaMcData->Write();
1180   }
1181
1182
1183   // 3. Pions
1184   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1185     TH1D *pionRaw     = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1186     pionRaw->Sumw2();
1187     TH1D *pionMC      = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1188     pionMC->Sumw2();
1189     TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1190     pionRawPrim->Sumw2();
1191     TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1192     pionSecReac->Sumw2();
1193     TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1194     pionSecWeak->Sumw2();
1195     TH1D *pionMis     = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1196     pionMis->Sumw2();
1197     pionSecReac->Add(pionSecWeak,-1);
1198     //
1199     // 1. total (physical) efficiency
1200     //
1201     TH1D *pionEffTot = new TH1D(*pionRaw);
1202     pionEffTot->Divide(pionRaw,pionMC);
1203     pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1204     pionEffTot->Write();
1205     //
1206     // 2. contributions from secondaries created by interactions
1207     //
1208     TH1D *pionEffSecReac = new TH1D(*pionRaw);
1209     TH1D *dummy = new TH1D(*pionRaw);
1210     dummy->Add(pionSecReac,-1);
1211     pionEffSecReac->Divide(pionRaw,dummy);
1212     pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1213     pionEffSecReac->Write();
1214     delete dummy;
1215     //
1216     // 3. contributions from secondaries from weak decays
1217     //
1218     TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1219     TH1D *dummy2 = new TH1D(*pionRaw);
1220     dummy2->Add(pionSecWeak,-1);
1221     pionEffSecWeak->Divide(pionRaw,dummy2);
1222     pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1223     pionEffSecWeak->Write();
1224     delete dummy2;
1225     //
1226     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1227     //
1228     TH1D *pionEffAbs = new TH1D(*pionRaw);
1229     pionEffAbs->Divide(pionRawPrim,pionMC);
1230     pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1231     pionEffAbs->Write();
1232     //
1233     // 5. contamination
1234     //
1235     TH1D *pionEffMis = new TH1D(*pionMis);
1236     pionEffMis->Divide(pionMis,pionRaw);
1237     pionEffMis->SetName(Form("PionContam_%i",iEta));
1238     pionEffMis->Write();
1239
1240     //
1241     // PROCESS THE REAL DATA FROM HERE ON
1242     //
1243     TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1244     histPtPionEta->Write();
1245     //
1246     // Apply corrections as you would like them to have...
1247     //
1248     TH1D *pionFinal = new TH1D(*histPtPionEta);
1249     pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1250     pionFinal->Sumw2();
1251     pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1252     pionFinal->Write();
1253     //
1254     // if available save also the MC truth for consistency checks
1255     //
1256     TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1257     histPtPionEtaMcData->Write();
1258
1259   }
1260
1261   //
1262   // Close the file after everthing is done
1263   //
1264   spectraFile.Close();
1265
1266 }
1267
1268
1269
1270