PWG2/SPECTRA/AliAnalysisTaskChargedHadronSpectra.cxx

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