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   // sort pT-bins ..
 212   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, -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};
 213   Int_t indexes[kPtBins+1];
 214   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 215   Double_t binsPt[kPtBins+1];
 216   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 217
 218
 219   // MC histograms
 220   fHistPtMCKaon = new TH3F("HistPtMCKaon", "PtEtaKaon; mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 221   fHistPtMCProton = new TH3F("HistPtMCProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 222   fHistPtMCPion = new TH3F("HistPtMCPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 223   //
 224   fHistPtMCKaon->GetZaxis()->Set(kPtBins, binsPt);
 225   fHistPtMCProton->GetZaxis()->Set(kPtBins, binsPt);
 226   fHistPtMCPion->GetZaxis()->Set(kPtBins, binsPt);
 227   //
 228   fListHist->Add(fHistPtMCKaon);
 229   fListHist->Add(fHistPtMCProton);
 230   fListHist->Add(fHistPtMCPion);
 231   //
 232   // reconstructed particle histograms
 233   fHistPtEtaKaon = new TH3F("HistPtEtaKaon", "PtEtaKaon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 234   fHistPtEtaKaonNoKink = new TH3F("HistPtEtaKaonNoKink", "PtEtaKaonNoKink;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 235   fHistPtEtaProton = new TH3F("HistPtEtaProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 236   fHistPtEtaProtonDCA = new TH3F("HistPtEtaProtonDCA", "PtEtaProton;DCA (cm); #eta; p_{T} (GeV)",1000,0,50,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 237   fHistPtEtaPion = new TH3F("HistPtEtaPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 238   //
 239   fHistPtEtaKaon->GetZaxis()->Set(kPtBins, binsPt);
 240   fHistPtEtaKaonNoKink->GetZaxis()->Set(kPtBins, binsPt);
 241   fHistPtEtaProton->GetZaxis()->Set(kPtBins, binsPt);
 242   fHistPtEtaProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
 243   fHistPtEtaPion->GetZaxis()->Set(kPtBins, binsPt);
 244   //
 245   fListHist->Add(fHistPtEtaKaon);
 246   fListHist->Add(fHistPtEtaKaonNoKink);
 247   fListHist->Add(fHistPtEtaProton);
 248   fListHist->Add(fHistPtEtaProtonDCA);
 249   fListHist->Add(fHistPtEtaPion);
 250
 251   // histograms for the classical analysis
 252   fHistClassicKaon = new TH3F("HistClassicKaon", "PtEtaKaon;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 253   fHistClassicProton = new TH3F("HistClassicProton", "PtEtaProton;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 254   fHistClassicPion = new TH3F("HistClassicPion", "PtEtaPion;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
 255   //
 256   fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
 257   fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
 258   fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
 259   //
 260   fListHist->Add(fHistClassicKaon);
 261   fListHist->Add(fHistClassicProton);
 262   fListHist->Add(fHistClassicPion);
 263   // histograms of general interest
 264   fDeDx = new TH3F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.1,100.,1000,0.,1000,2,-2,2);
 265   BinLogX(fDeDx);
 266   fListHist->Add(fDeDx);
 267   fHistTrackPerEvent = new TH2F("HistTrackPerEvent", "Tracks per event;Number of Tracks;Number of Events; code",101,-0.5,100,10,-0.5,9.5);
 268   fListHist->Add(fHistTrackPerEvent);
 269   fHistTrackPerEventMC  = new TH3F("HistTrackPerEventMC", "Tracks per event;code;TrackPerEvent;IsSelected",10,-0.5,9,101,-0.5,100,2,-0.5,1.5);
 270   fListHist->Add(fHistTrackPerEventMC);
 271   fSecProtons = new TH2F("SecProtons", "xy;x;y",100,-10,10,100,-10,10);
 272   fListHist->Add(fSecProtons);
 273   fVertexZ = new TH3F("VertexZ", "vertex position;x (cm); y (cm); z (cm); counts",100,-5,5,100,-5,5,100,-25,25);
 274   fListHist->Add(fVertexZ);
 275   //
 276   fHistEtaNcls = new TH2F("HistEtaNcls", "EtaNcls;#eta;Ncls",100,-1.5,1.5,80,0,160);
 277   fListHist->Add(fHistEtaNcls);
 278   fHistEtaPhi = new TH2F("HistEtaPhi", "EtaNcls;#eta;#phi",100,-1.5,1.5,100,0,2*TMath::Pi());
 279   fListHist->Add(fHistEtaPhi);
 280
 281   // histograms for a refined efficiency study
 282   fHistEffProton = new TH3F("HistEffProton", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 283   fListHist->Add(fHistEffProton);
 284   fHistEffProtonDCA  = new TH3F("HistEffProtonDCA", "Eff;p_{T} (GeV); code",10,-0.5,9.5,200,0,15,kPtBins,-kPtMax,kPtMax);
 285   fListHist->Add(fHistEffProtonDCA);
 286   fHistEffPion = new TH3F("HistEffPion", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 287   fListHist->Add(fHistEffPion);
 288   fHistEffKaon = new TH3F("HistEffKaon", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
 289   fListHist->Add(fHistEffKaon);
 290   //
 291   fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
 292   fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
 293   fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
 294   fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
 295   //
 296   // create the histograms with all necessary information
 297   //
 298   //  (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rapPion, (5.) rapKaon, (6.) rapProton, (7.) dca, (8.) nclIter1, (9.) pullPion, (10.) pullKaon, (11.) pullProton
 299   //                               0,   1, 2, 3,  4 ,   5,   6, 7,  8,  9, 10, 11
 300   Int_t    binsHistReal[12] = {  301,  16,20, 2,  16,  16,  16,60, 10,100,100,100};
 301   Double_t xminHistReal[12] = { -0.5,-0.8, 0,-2,-0.8,-0.8,-0.8, 0, 60,-10,-10,-10};
 302   Double_t xmaxHistReal[12] = {300.5, 0.8, 2, 2, 0.8, 0.8, 0.8,15,160, 10, 10, 10};
 303   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);
 304   fListHist->Add(fHistRealTracks);
 305   //
 306   // (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rap, (5.) dca, (6.) ID, (7.) code
 307   //                            0,   1, 2, 3,   4, 5,   6,  7
 308   Int_t    binsHistMC[8] = {  301,  16,20, 2,  16,60,   5,  10};
 309   Double_t xminHistMC[8] = { -0.5,-0.8, 0,-2,-0.8, 0,-0.5,-0.5};
 310   Double_t xmaxHistMC[8] = {300.5, 0.8, 2, 2, 0.8,15, 4.5, 9.5};
 311   fHistMCparticles = new THnSparseS("fHistMCparticles"," (0.) multiplicity, (1.) eta, (2.) pT, (3.) sign, (4.) rap, (5.) dca, (6.) ID, (7.) code",8,binsHistMC,xminHistMC,xmaxHistMC);
 312   fListHist->Add(fHistMCparticles);
 313   //
 314   //
 315   //
 316
 317
 318 }
 319
 320 //________________________________________________________________________
 321 void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *)
 322 {
 323   //
 324   // main event loop
 325   //
 326   AliLog::SetGlobalLogLevel(AliLog::kError);
 327   //
 328   // Check Monte Carlo information and other access first:
 329   //
 330   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 331   if (!eventHandler) {
 332     //Printf("ERROR: Could not retrieve MC event handler");
 333     if (fMCtrue) return;
 334   }
 335   //
 336   AliMCEvent* mcEvent = 0x0;
 337   AliStack* stack = 0x0;
 338   if (eventHandler) mcEvent = eventHandler->MCEvent();
 339   if (!mcEvent) {
 340     //Printf("ERROR: Could not retrieve MC event");
 341     if (fMCtrue) return;
 342   }
 343   if (fMCtrue) {
 344     stack = mcEvent->Stack();
 345     if (!stack) return;
 346   }
 347   //
 348   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
 349   if (!fESD) {
 350     //Printf("ERROR: fESD not available");
 351     return;
 352   }
 353
 354   if (!fESDtrackCuts) {
 355     Printf("ERROR: fESDtrackCuts not available");
 356     return;
 357   }
 358
 359   Int_t trackCounter = 0;
 360   //
 361   // check if event is selected by physics selection class
 362   //
 363   fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks(),0);
 364   Bool_t isSelected = kFALSE;
 365   isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 366   if (isSelected) fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks(),1);
 367   //
 368   // monitor vertex position
 369   //
 370   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 371   if(vertex->GetNContributors()<1) {
 372     // SPD vertex
 373     vertex = fESD->GetPrimaryVertexSPD();
 374     if(vertex->GetNContributors()<1) vertex = 0x0;
 375   }
 376
 377   // 1st track loop to determine multiplicities
 378   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 379     AliESDtrack *track = 0x0;
 380     if (fTrackingMode == 0) track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
 381     if (fTrackingMode == 1) track = new AliESDtrack(*fESD->GetTrack(i));  // alternative: global track GLOBAL GLOBAL  GLOBAL GLOBAL  GLOBAL GLOBAL  GLOBAL
 382     if (!track) continue;
 383     if (fESDtrackCuts->AcceptTrack(track) &&  TMath::Abs(track->Eta())< 0.9) {
 384       trackCounter++;
 385       //
 386     }
 387     delete track;
 388   }
 389
 390   fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
 391   fHistTrackPerEvent->Fill(trackCounter,2);
 392
 393   if (fMCtrue) {
 394     //
 395     //
 396     //
 397     AliHeader * header = mcEvent->Header();
 398     Int_t processtype = GetPythiaEventProcessType(header);
 399     // non diffractive
 400     if (processtype !=92 && processtype !=93 && processtype != 94) fHistTrackPerEventMC->Fill(1., trackCounter, isSelected);
 401     // single diffractive
 402     if ((processtype == 92 || processtype == 93)) fHistTrackPerEventMC->Fill(2., trackCounter, isSelected);
 403     // double diffractive
 404     if (processtype == 94) fHistTrackPerEventMC->Fill(3., trackCounter, isSelected);
 405     //
 406     Int_t mult = trackCounter;//stack->GetNprimary();
 407     //
 408     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 409       TParticle * trackMC = stack->Particle(i);
 410       Int_t pdg = trackMC->GetPdgCode();
 411       //
 412       Double_t xv = trackMC->Vx();
 413       Double_t yv = trackMC->Vy();
 414       Double_t zv = trackMC->Vz();
 415       Double_t dxy = 0;
 416       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
 417       Double_t dz = 0;
 418       dz = TMath::Abs(zv); // so stupid to avoid warnings
 419       //
 420       // vertex cut - selection of primaries
 421       //
 422       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
 423       //
 424       if (!stack->IsPhysicalPrimary(i)) continue;
 425       //
 426       if (pdg == 321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select K+ only
 427       if (pdg == 211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select Pi+ only
 428       if (pdg == 2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Y()), trackMC->Pt()); // select p+ only
 429       //
 430       if (pdg == -321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select K- only
 431       if (pdg == -211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select Pi- only
 432       if (pdg == -2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // select p- only
 433       //
 434       // special Kaon checks - those which decay before entering the TPC fiducial volume
 435       // <-> should be revisited by looping from first to last daugther and check if UniqueID is 4
 436       //
 437       if (TMath::Abs(pdg)==321) {
 438         TParticle * trackDaughterMC = stack->Particle(TMath::Abs(trackMC->GetFirstDaughter()));
 439         Int_t pdgDaughter = trackDaughterMC->GetPdgCode();
 440         if (TMath::Abs(pdgDaughter) == 13 && pdg == 321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Y()), trackMC->Pt()); // Kaon control hist
 441         if (TMath::Abs(pdgDaughter) == 13 && pdg == -321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Y()), -trackMC->Pt()); // Kaon control hist
 442       }
 443     }
 444   }
 445   //
 446   //end MC tracks loop <-> now check event selection criteria
 447   //
 448   if (!isSelected) {
 449     PostData(1, fListHist);
 450     return;
 451   }
 452   //
 453   if (!vertex) {
 454     fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
 455     PostData(1, fListHist);
 456     return;
 457   } else {
 458     fVertexZ->Fill(vertex->GetXv(),vertex->GetYv(),vertex->GetZv());
 459     if (TMath::Abs(vertex->GetZv()) > 10) {
 460       fHistTrackPerEventMC->Fill(0., trackCounter, isSelected);
 461       PostData(1, fListHist);
 462       return;
 463     }
 464   }
 465   //
 466   // tarck loop
 467   //
 468   const Float_t kNsigmaCut = 3;
 469   const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
 470
 471   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 472
 473   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 474     //
 475     // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
 476     //
 477     AliESDtrack *track = 0x0;
 478     if (fTrackingMode == 0) track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
 479     if (fTrackingMode == 1) track = new AliESDtrack(*fESD->GetTrack(i)); // alternative: global track GLOBAL GLOBAL  GLOBAL GLOBAL
 480     if (!track) continue;
 481     AliESDtrack *trackDummy = fESD->GetTrack(i);
 482     if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
 483       delete track;
 484       continue;
 485     }
 486     //
 487     if (track->GetTPCNclsIter1() < 80) {
 488       delete track;
 489       continue;
 490     }
 491     //
 492     AliExternalTrackParam trackIn(*trackDummy->GetInnerParam());
 493     Double_t ptot = trackIn.GetP(); // momentum for dEdx determination
 494     Double_t pT = track->Pt();
 495     Double_t eta = TMath::Abs(track->Eta());
 496     //
 497     fHistEtaNcls->Fill(track->Eta(),track->GetTPCNcls());
 498     fHistEtaPhi->Fill(track->Eta(),track->Phi());
 499     //
 500     // cut for dead regions in the detector
 501     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
 502     //
 503     // 2.a) apply some standard track cuts according to general recommendations
 504     //
 505     if (!fESDtrackCuts->AcceptTrack(track)) {
 506       //
 507       if (fMCtrue) {
 508         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 509         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 510         if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(6,eta,track->GetSign()*pT); // Kaon control hist
 511       }
 512       //
 513       delete track;
 514       continue;
 515     }
 516
 517     if (fTrackingMode == 1) {
 518       /*if (!SelectOnImpPar(track)) {
 519         delete track;
 520         continue;
 521         }*/
 522     }
 523
 524     //
 525     // calculate rapidities and kinematics
 526     //
 527     //
 528     Double_t pvec[3];
 529     track->GetPxPyPz(pvec);
 530     Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
 531     Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
 532     Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
 533     //
 534     Double_t rapPion = TMath::Abs(0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2])));
 535     Double_t rapKaon = TMath::Abs(0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2])));
 536     Double_t rapProton = TMath::Abs(0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2])));
 537     //
 538     // 3. make the PID
 539     //
 540     Double_t sign = track->GetSign();
 541     Double_t tpcSignal = track->GetTPCsignal();
 542     //
 543     if (eta < 0.8 && track->GetTPCsignalN() > 70) fDeDx->Fill(ptot, tpcSignal, sign);
 544     //
 545     // 3.a. calculate expected signals
 546     //
 547     Float_t sigKaon     = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kKaon);
 548     Float_t sigProton   = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kProton);
 549     Float_t sigPion     = fESDpid->GetTPCResponse().GetExpectedSignal(ptot, AliPID::kPion);
 550     //
 551     // 3.b. fill histograms
 552     //
 553     fHistClassicKaon->Fill(sign*pT,rapKaon,(tpcSignal-sigKaon)/sigKaon);
 554     fHistClassicProton->Fill(sign*pT,rapProton,(tpcSignal-sigProton)/sigProton);
 555     fHistClassicPion->Fill(sign*pT,rapPion,(tpcSignal-sigPion)/sigPion);
 556     //
 557     /* 2sigma PID with 2sigma eff correction! */
 558     // PION
 559     if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kPion)) < kNsigmaCut) {
 560       fHistPtEtaPion->Fill(trackCounter, rapPion, sign*pT, k2sigmaCorr);
 561       //
 562       if (trackCounter < 300 && fMCtrue) {
 563         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 564         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 565         if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,rapPion,sign*pT,k2sigmaCorr);
 566         if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 567           fHistEffPion->Fill(1,rapPion,sign*pT,k2sigmaCorr);
 568           if (trackMC->GetFirstMother() > 0) {
 569             TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 570             Int_t code = TMath::Abs(trackMother->GetPdgCode());
 571             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);
 572           }
 573         }
 574         if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,rapPion,sign*pT,k2sigmaCorr);
 575         if (TMath::Abs(trackMC->GetPdgCode()) == 13)  fHistEffPion->Fill(4,rapPion,sign*pT,k2sigmaCorr);
 576       }
 577     }
 578     // KAON
 579     if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < kNsigmaCut) {
 580       fHistPtEtaKaon->Fill(trackCounter, rapKaon, sign*pT, k2sigmaCorr);
 581       //
 582       if (trackCounter < 300 && fMCtrue) {
 583         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 584         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 585         if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,rapKaon,sign*pT,k2sigmaCorr);
 586         if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 587           fHistEffKaon->Fill(1,rapKaon,sign*pT,k2sigmaCorr);
 588           if (trackMC->GetFirstMother() > 0) {
 589             TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 590             Int_t code = TMath::Abs(trackMother->GetPdgCode());
 591             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);
 592           }
 593         }
 594         if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,rapKaon,sign*pT,k2sigmaCorr);
 595       }
 596     }
 597     // KAON NO KINK
 598     if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < kNsigmaCut && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, rapKaon, sign*pT, k2sigmaCorr);
 599     // PROTON
 600     if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)) < kNsigmaCut) {
 601       fHistPtEtaProton->Fill(trackCounter, rapProton, sign*pT, k2sigmaCorr);
 602       //
 603       track->GetImpactParameters(dca,cov);
 604       Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
 605       fHistPtEtaProtonDCA->Fill(primVtxDCA, rapProton, sign*pT, k2sigmaCorr);
 606       //
 607       if (trackCounter < 300 && fMCtrue) {
 608         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 609         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 610         if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 611           fHistEffProton->Fill(0.,rapProton,sign*pT,k2sigmaCorr);
 612           if (rapProton < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
 613         }
 614         if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
 615           fHistEffProton->Fill(1,rapProton,sign*pT,k2sigmaCorr);
 616           if (rapProton < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
 617           if (trackMC->GetFirstMother() > 0) {
 618             TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
 619             Int_t code = TMath::Abs(trackMother->GetPdgCode());
 620             if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
 621               fHistEffProton->Fill(3,rapProton,sign*pT,k2sigmaCorr);
 622               if (rapProton < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
 623             }
 624             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());
 625           }
 626         }
 627         if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,rapProton,sign*pT,k2sigmaCorr);
 628       }
 629     }
 630
 631     delete track;
 632
 633   } // end of track loop
 634
 635   // Post output data
 636   PostData(1, fListHist);
 637
 638 }
 639
 640
 641 //________________________________________________________________________
 642 void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
 643 {
 644   // Draw result to the screen
 645   // Called once at the end of the query
 646   Printf("*** CONSTRUCTOR CALLED ****");
 647
 648 }
 649
 650
 651 //________________________________________________________________________
 652 Bool_t AliAnalysisTaskChargedHadronSpectra::SelectOnImpPar(AliESDtrack* t) {
 653   //
 654   // cut on transverse impact parameter
 655   //
 656   Float_t d0z0[2],covd0z0[3];
 657   t->GetImpactParameters(d0z0,covd0z0);
 658   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
 659   Float_t d0max = 7.*sigma;
 660   //
 661   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
 662   if (t->Pt() > 1) sigmaZ = 0.0216;
 663   Float_t d0maxZ = 5.*sigmaZ;
 664   //
 665   if (TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
 666   return kFALSE;
 667 }
 668
 669
 670 //________________________________________________________________________
 671 void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
 672
 673   // Method for the correct logarithmic binning of histograms
 674
 675   TAxis *axis = h->GetXaxis();
 676   int bins = axis->GetNbins();
 677
 678   Double_t from = axis->GetXmin();
 679   Double_t to = axis->GetXmax();
 680   Double_t *newBins = new Double_t[bins + 1];
 681
 682   newBins[0] = from;
 683   Double_t factor = pow(to/from, 1./bins);
 684
 685   for (int i = 1; i <= bins; i++) {
 686    newBins[i] = factor * newBins[i-1];
 687   }
 688   axis->Set(bins, newBins);
 689   delete [] newBins;
 690
 691 }
 692
 693
 694 //________________________________________________________________________
 695 Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
 696   //
 697   // get the process type of the event.
 698   //
 699
 700   // can only read pythia headers, either directly or from cocktalil header
 701   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 702
 703   if (!pythiaGenHeader) {
 704
 705     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 706     if (!genCocktailHeader) {
 707       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 708       return -1;
 709     }
 710
 711     TList* headerList = genCocktailHeader->GetHeaders();
 712     if (!headerList) {
 713       return -1;
 714     }
 715
 716     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 717       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 718       if (pythiaGenHeader)
 719         break;
 720     }
 721
 722     if (!pythiaGenHeader) {
 723       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
 724       return -1;
 725     }
 726   }
 727
 728   if (adebug) {
 729     //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
 730   }
 731
 732   return pythiaGenHeader->ProcessType();
 733 }
 734
 735
 736 //________________________________________________________________________
 737 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBinLow, Int_t EtaBinHigh,const  Double_t * AlephParams) {
 738   //
 739   // The multiple Gauss fitting is happening here...
 740   // input histogram is of the form (Pt, eta, difference in dEdx)
 741   //
 742
 743   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
 744   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
 745   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
 746   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
 747   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
 748   //
 749   foProton->SetParameters(AlephParams);
 750   foPion->SetParameters(AlephParams);
 751   foElec->SetParameters(AlephParams);
 752   foKaon->SetParameters(AlephParams);
 753   foMuon->SetParameters(AlephParams);
 754
 755
 756   const Double_t kSigma = 0.06; // expected resolution (roughly)
 757
 758   const Int_t kPtBins = 2*30/*input->GetXaxis()->GetNbins()*/;
 759   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
 760   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
 761
 762   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
 763   // sort pT-bins ..
 764   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, -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};
 765   Int_t indexes[kPtBins+1];
 766   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 767   Double_t binsPt[kPtBins+1];
 768   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 769   histPt->GetXaxis()->Set(kPtBins, binsPt);
 770   //
 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   // sort pT-bins ..
 902   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, -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};
 903   Int_t indexes[kPtBins+1];
 904   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
 905   Double_t binsPt[kPtBins+1];
 906   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
 907   histPt->GetXaxis()->Set(kPtBins, binsPt);
 908   //
 909
 910   TCanvas * canvMany = new TCanvas("canvManyPion","canvManyPion");
 911   canvMany->Print("canvManyPion.ps[");
 912
 913   for(Int_t x=1; x < kPtBins+1; x++) {
 914     Double_t pT = input->GetXaxis()->GetBinCenter(x);
 915     for(Int_t y=1; y < kEtaBins+1; y++) {
 916       Double_t rapidity =  input->GetYaxis()->GetBinCenter(y);
 917       Double_t eta = TMath::ASinH(TMath::Sqrt(1 + (0.1395*0.1395)/(pT*pT))*TMath::SinH(rapidity));
 918       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
 919       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
 920       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y + EtaBinHigh - EtaBinLow);
 921       histDeDx->SetTitle(Form("pt_%f",pT));
 922       Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
 923       //
 924       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
 925       Double_t parameters[12] = {maxYield/2.,0,kSigma,
 926                                  maxYield/2.,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foPion->Eval(ptot)/foPion->Eval(ptot)),
 927                                  maxYield/2.,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)),
 928                                  maxYield/2.,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot))};
 929       gausFit->SetParameters(parameters);
 930       gausFit->SetParLimits(0,0.,maxYield);
 931       gausFit->SetParLimits(1,-0.05,0.05);
 932       gausFit->SetParLimits(2,0.04,0.08);
 933       //
 934       Double_t protonPosition = (foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
 935       gausFit->SetParLimits(3,0.,maxYield);
 936       gausFit->SetParLimits(4, 0.95*protonPosition, 1.05*protonPosition);
 937       gausFit->SetParLimits(5,0.8*kSigma*(foProton->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foProton->Eval(ptot)/foPion->Eval(ptot)));
 938       //
 939       Double_t elecPosition = (foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
 940       gausFit->SetParLimits(6,0.,maxYield);
 941       gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
 942       gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foPion->Eval(ptot)));
 943       //
 944       Double_t kaonPosition = (foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot);
 945       gausFit->SetParLimits(9,0.,maxYield);
 946       gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
 947       gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foPion->Eval(ptot)));
 948       histDeDx->Fit("gausFit","QNR");
 949       gausFit->GetParameters(parameters);
 950       // visualization
 951       if (y == EtaBinLow) {
 952         canvMany->cd(x);
 953         gPad->SetLogy();
 954         histDeDx->SetMarkerStyle(21);
 955         histDeDx->SetMarkerSize(0.7);
 956         histDeDx->Draw("E");
 957         gausFit->SetLineWidth(2);
 958         gausFit->Draw("same");
 959         //
 960         TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
 961         TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
 962         TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
 963         TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
 964         gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
 965         gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
 966         gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
 967         gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
 968         //
 969         gausFit0.SetLineColor(4);
 970         gausFit1.SetLineColor(2);
 971         gausFit2.SetLineColor(6);
 972         gausFit3.SetLineColor(8);
 973         //
 974         gausFit0.SetLineWidth(1);
 975         gausFit1.SetLineWidth(1);
 976         gausFit2.SetLineWidth(1);
 977         gausFit3.SetLineWidth(1);
 978         //
 979         gausFit0.Draw("same");
 980         gausFit1.Draw("same");
 981         gausFit2.Draw("same");
 982         gausFit3.Draw("same");
 983         if (TMath::Abs(pT) < 2) canvMany->Print("canvManyPion.ps");
 984       }
 985
 986       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
 987       delete gausFit;
 988       //
 989       if (y == EtaBinLow && yield > 0) histPt->Fill(pT,yield);
 990       //delete gausYield;
 991     }
 992   }
 993   canvMany->Print("canvManyPion.ps]");
 994
 995   TCanvas * canvPt = new TCanvas();
 996   canvPt->cd();
 997   histPt->Draw();
 998
 999   return histPt;
1000 }
1001
1002
1003 //________________________________________________________________________
1004 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBinLow, Int_t EtaBinHigh,const  Double_t * AlephParams) {
1005   //
1006   // The multiple Gauss fitting is happening here...
1007   // input histogram is of the form (Pt, eta, difference in dEdx)
1008   //
1009
1010   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
1011   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
1012   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
1013   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
1014   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
1015   //
1016   foProton->SetParameters(AlephParams);
1017   foPion->SetParameters(AlephParams);
1018   foElec->SetParameters(AlephParams);
1019   foKaon->SetParameters(AlephParams);
1020   foMuon->SetParameters(AlephParams);
1021
1022
1023   const Double_t kSigma = 0.055; // expected resolution (roughly)
1024
1025   const Int_t kPtBins = 2*30/*input->GetXaxis()->GetNbins()*/;
1026   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
1027   //const Int_t kEtaBins = input->GetYaxis()->GetNbins();
1028
1029   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
1030   // sort pT-bins ..
1031   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, -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};
1032   Int_t indexes[kPtBins+1];
1033   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
1034   Double_t binsPt[kPtBins+1];
1035   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
1036   histPt->GetXaxis()->Set(kPtBins, binsPt);
1037   //
1038
1039   TCanvas * canvMany = new TCanvas("canvManyKaon","canvManyKaon");
1040   canvMany->Print("canvManyKaon.ps[");
1041
1042   for(Int_t x=1; x < kPtBins+1; x++) {
1043     Double_t pT = input->GetXaxis()->GetBinCenter(x);
1044     for(Int_t y=1; y < 2/*kEtaBins+1*/; y++) {
1045       Double_t rapidity =  input->GetYaxis()->GetBinCenter(y);
1046       Double_t eta = TMath::ASinH(TMath::Sqrt(1 + (0.4936*0.4936)/(pT*pT))*TMath::SinH(rapidity));
1047       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1048       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
1049       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y + EtaBinHigh - EtaBinLow);
1050       histDeDx->SetTitle(Form("pt_%f",pT));
1051       Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
1052       //
1053       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1 /*-5*kSigma*/, 1/*5*kSigma*/);
1054       Double_t parameters[12] = {maxYield*0.1,0,kSigma,
1055                                  maxYield*0.8,(foPion->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)),
1056                                  maxYield*0.01,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)),
1057                                  maxYield*0.1,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot))};
1058       gausFit->SetParameters(parameters);
1059       gausFit->SetParLimits(0,0.,maxYield);
1060       gausFit->SetParLimits(1,-0.15,0.15);
1061       gausFit->SetParLimits(2,0.04,0.08);
1062       //
1063       Double_t pionPosition = (foPion->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
1064       //Double_t elecPosition = (foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
1065       Double_t protonPosition = (foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot);
1066       //
1067       if (TMath::Abs(pT) < 0.4) {
1068         cout << " fixing the parameters " << endl;
1069         gausFit->SetParameter(3,0.); gausFit->SetParameter(6,0.); gausFit->SetParameter(9,0.);
1070         gausFit->FixParameter(3,0); gausFit->FixParameter(6,0); gausFit->FixParameter(9,0);
1071         //gausFit->SetParLimits(0,0.,maxYield);
1072         //gausFit->SetParLimits(1,-0.15,0.15);
1073         //gausFit->SetParLimits(2,0.03,0.18);
1074       } else {
1075         //
1076         gausFit->SetParLimits(0,0.,0.5*maxYield);
1077         gausFit->SetParLimits(1,-0.05,0.05);
1078         gausFit->SetParLimits(2,0.04,0.08);
1079         //
1080         gausFit->SetParLimits(3,0.,maxYield);
1081         gausFit->SetParLimits(4, 0.5*pionPosition, 1.5*pionPosition);
1082         gausFit->SetParLimits(5,0.6*kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foPion->Eval(ptot)/foKaon->Eval(ptot)));
1083         //
1084         gausFit->FixParameter(6,0);
1085         gausFit->FixParameter(7,0);
1086         gausFit->FixParameter(8,0);
1087         /*
1088           gausFit->SetParLimits(6,0.,0.2*maxYield);
1089           gausFit->SetParLimits(7, 0.9*elecPosition, 1.1*elecPosition);
1090           gausFit->SetParLimits(8,0.6*kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foElec->Eval(ptot)/foKaon->Eval(ptot)));
1091         */
1092         //
1093         gausFit->SetParLimits(9,0.,0.2*maxYield);
1094         gausFit->SetParLimits(10, 0.9*protonPosition, 1.1*protonPosition);
1095         gausFit->SetParLimits(11,0.6*kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot)),1.4*kSigma*(foProton->Eval(ptot)/foKaon->Eval(ptot)));
1096       }
1097       histDeDx->Fit("gausFit","QNR");
1098       gausFit->GetParameters(parameters);
1099       // visualization
1100       if (y == EtaBinLow) {
1101         canvMany->cd(x);
1102         //gPad->SetLogy();
1103         histDeDx->SetMarkerStyle(21);
1104         histDeDx->SetMarkerSize(0.7);
1105         histDeDx->Draw("E");
1106         gausFit->SetLineWidth(2);
1107         gausFit->Draw("same");
1108         //
1109         TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
1110         TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
1111         TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
1112         TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
1113         gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
1114         gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
1115         gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
1116         gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
1117         //
1118         gausFit0.SetLineColor(4);
1119         gausFit1.SetLineColor(2);
1120         gausFit2.SetLineColor(6);
1121         gausFit3.SetLineColor(8);
1122         //
1123         gausFit0.SetLineWidth(1);
1124         gausFit1.SetLineWidth(1);
1125         gausFit2.SetLineWidth(1);
1126         gausFit3.SetLineWidth(1);
1127         //
1128         gausFit0.Draw("same");
1129         gausFit1.Draw("same");
1130         gausFit2.Draw("same");
1131         gausFit3.Draw("same");
1132         if (TMath::Abs(pT) < 2) canvMany->Print("canvManyKaon.ps");
1133         //
1134         // propaganda slots for the paper
1135         //
1136         if (pT > 0.374 && pT < 0.376) {
1137           TFile f1("slice1.root","RECREATE");
1138           canvMany->Write();
1139           f1.Close();
1140         }
1141         if (pT > 0.524 && pT < 0.526) {
1142           TFile f1("slice2.root","RECREATE");
1143           canvMany->Write();
1144           f1.Close();
1145         }
1146         if (pT > 0.674 && pT < 0.676) {
1147           TFile f1("slice3.root","RECREATE");
1148           canvMany->Write();
1149           f1.Close();
1150         }
1151         if (pT > 0.624 && pT < 0.626) {
1152           TFile f1("slice4.root","RECREATE");
1153           canvMany->Write();
1154           f1.Close();
1155         }
1156       }
1157
1158       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
1159       cout << pT << " res.: " << parameters[2] << " " << " yield: " << yield << endl;
1160       delete gausFit;
1161       //
1162       if (y == EtaBinLow && yield > 0) histPt->Fill(pT,yield);
1163       //delete gausYield;
1164     }
1165   }
1166   canvMany->Print("canvManyKaon.ps]");
1167
1168   TCanvas * canvPt = new TCanvas();
1169   canvPt->cd();
1170   histPt->Draw();
1171
1172   return histPt;
1173 }
1174
1175
1176 //________________________________________________________________________
1177 void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const  TList * ListOfHistogramsData, const Char_t *filename) {
1178
1179   // Define ranges
1180   const Int_t kMultLow  = 1;
1181   const Int_t kMultHigh = 10;
1182
1183   const Int_t kEtaHigh  = 4;
1184
1185   // Extract histograms for the MC efficiency study
1186   TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1187   TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1188   TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
1189
1190   TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
1191   TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1192   TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
1193
1194   TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1195   TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1196   TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1197
1198   // Extract data for the final spectra
1199   TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
1200   TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1201   TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
1202
1203   // "MC" of the real data for consistency checks
1204   TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1205   TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1206   TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
1207
1208
1209   TFile spectraFile(filename,"RECREATE");
1210
1211   // 1. Protons
1212   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1213     TH1D *protonRaw     = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1214     protonRaw->Sumw2();
1215     TH1D *protonMC      = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1216     protonMC->Sumw2();
1217     TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1218     protonRawPrim->Sumw2();
1219     TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1220     protonSecReac->Sumw2();
1221     TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1222     protonSecWeak->Sumw2();
1223     TH1D *protonMis     = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1224     protonMis->Sumw2();
1225     protonSecReac->Add(protonSecWeak,-1);
1226     //
1227     // 1. total (physical) efficiency
1228     //
1229     TH1D *protonEffTot = new TH1D(*protonRaw);
1230     protonEffTot->Divide(protonRaw,protonMC);
1231     protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1232     protonEffTot->Write();
1233     //
1234     // 2. contributions from secondaries created by interactions
1235     //
1236     TH1D *protonEffSecReac = new TH1D(*protonRaw);
1237     TH1D *dummy = new TH1D(*protonRaw);
1238     dummy->Add(protonSecReac,-1);
1239     protonEffSecReac->Divide(protonRaw,dummy);
1240     protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1241     protonEffSecReac->Write();
1242     delete dummy;
1243     //
1244     // 3. contributions from secondaries from weak decays
1245     //
1246     TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1247     TH1D *dummy2 = new TH1D(*protonRaw);
1248     dummy2->Add(protonSecWeak,-1);
1249     protonEffSecWeak->Divide(protonRaw,dummy2);
1250     protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1251     protonEffSecWeak->Write();
1252     delete dummy2;
1253     //
1254     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1255     //
1256     TH1D *protonEffAbs = new TH1D(*protonRaw);
1257     protonEffAbs->Divide(protonRawPrim,protonMC);
1258     protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1259     protonEffAbs->Write();
1260     //
1261     // 5. contamination
1262     //
1263     TH1D *protonEffMis = new TH1D(*protonMis);
1264     protonEffMis->Divide(protonMis,protonRaw);
1265     protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1266     protonEffMis->Write();
1267
1268     //
1269     // PROCESS THE REAL DATA FROM HERE ON
1270     //
1271     TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1272     histPtProtonEta->Write();
1273     //
1274     // Apply corrections as you would like them to have...
1275     //
1276     TH1D *protonFinal = new TH1D(*histPtProtonEta);
1277     protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1278     protonFinal->Sumw2();
1279     protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1280     protonFinal->Write();
1281     //
1282     // if available save also the MC truth for consistency checks
1283     //
1284     TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1285     histPtProtonEtaMcData->Write();
1286   }
1287
1288   // 2. Kaons
1289   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1290     TH1D *kaonRaw     = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1291     kaonRaw->Sumw2();
1292     TH1D *kaonMC      = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1293     kaonMC->Sumw2();
1294     TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1295     kaonRawPrim->Sumw2();
1296     TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1297     kaonSecReac->Sumw2();
1298     TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1299     kaonSecWeak->Sumw2();
1300     TH1D *kaonMis     = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1301     kaonMis->Sumw2();
1302     kaonSecReac->Add(kaonSecWeak,-1);
1303     //
1304     // 1. total (physical) efficiency
1305     //
1306     TH1D *kaonEffTot = new TH1D(*kaonRaw);
1307     kaonEffTot->Divide(kaonRaw,kaonMC);
1308     kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1309     kaonEffTot->Write();
1310     //
1311     // 2. contributions from secondaries created by interactions
1312     //
1313     TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1314     TH1D *dummy = new TH1D(*kaonRaw);
1315     dummy->Add(kaonSecReac,-1);
1316     kaonEffSecReac->Divide(kaonRaw,dummy);
1317     kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1318     kaonEffSecReac->Write();
1319     delete dummy;
1320     //
1321     // 3. contributions from secondaries from weak decays
1322     //
1323     TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1324     TH1D *dummy2 = new TH1D(*kaonRaw);
1325     dummy2->Add(kaonSecWeak,-1);
1326     kaonEffSecWeak->Divide(kaonRaw,dummy2);
1327     kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1328     kaonEffSecWeak->Write();
1329     delete dummy2;
1330     //
1331     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1332     //
1333     TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1334     kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1335     kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1336     kaonEffAbs->Write();
1337     //
1338     // 5. contamination
1339     //
1340     TH1D *kaonEffMis = new TH1D(*kaonMis);
1341     kaonEffMis->Divide(kaonMis,kaonRaw);
1342     kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1343     kaonEffMis->Write();
1344
1345     //
1346     // PROCESS THE REAL DATA FROM HERE ON
1347     //
1348     TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1349     histPtKaonEta->Write();
1350     //
1351     // Apply corrections as you would like them to have...
1352     //
1353     TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1354     kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1355     kaonFinal->Sumw2();
1356     kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1357     kaonFinal->Write();
1358     //
1359     // if available save also the MC truth for consistency checks
1360     //
1361     TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1362     histPtKaonEtaMcData->Write();
1363   }
1364
1365
1366   // 3. Pions
1367   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1368     TH1D *pionRaw     = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1369     pionRaw->Sumw2();
1370     TH1D *pionMC      = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1371     pionMC->Sumw2();
1372     TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1373     pionRawPrim->Sumw2();
1374     TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1375     pionSecReac->Sumw2();
1376     TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1377     pionSecWeak->Sumw2();
1378     TH1D *pionMis     = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1379     pionMis->Sumw2();
1380     pionSecReac->Add(pionSecWeak,-1);
1381     //
1382     // 1. total (physical) efficiency
1383     //
1384     TH1D *pionEffTot = new TH1D(*pionRaw);
1385     pionEffTot->Divide(pionRaw,pionMC);
1386     pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1387     pionEffTot->Write();
1388     //
1389     // 2. contributions from secondaries created by interactions
1390     //
1391     TH1D *pionEffSecReac = new TH1D(*pionRaw);
1392     TH1D *dummy = new TH1D(*pionRaw);
1393     dummy->Add(pionSecReac,-1);
1394     pionEffSecReac->Divide(pionRaw,dummy);
1395     pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1396     pionEffSecReac->Write();
1397     delete dummy;
1398     //
1399     // 3. contributions from secondaries from weak decays
1400     //
1401     TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1402     TH1D *dummy2 = new TH1D(*pionRaw);
1403     dummy2->Add(pionSecWeak,-1);
1404     pionEffSecWeak->Divide(pionRaw,dummy2);
1405     pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1406     pionEffSecWeak->Write();
1407     delete dummy2;
1408     //
1409     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1410     //
1411     TH1D *pionEffAbs = new TH1D(*pionRaw);
1412     pionEffAbs->Divide(pionRawPrim,pionMC);
1413     pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1414     pionEffAbs->Write();
1415     //
1416     // 5. contamination
1417     //
1418     TH1D *pionEffMis = new TH1D(*pionMis);
1419     pionEffMis->Divide(pionMis,pionRaw);
1420     pionEffMis->SetName(Form("PionContam_%i",iEta));
1421     pionEffMis->Write();
1422
1423     //
1424     // PROCESS THE REAL DATA FROM HERE ON
1425     //
1426     TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1427     histPtPionEta->Write();
1428     //
1429     // Apply corrections as you would like them to have...
1430     //
1431     TH1D *pionFinal = new TH1D(*histPtPionEta);
1432     pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1433     pionFinal->Sumw2();
1434     pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1435     pionFinal->Write();
1436     //
1437     // if available save also the MC truth for consistency checks
1438     //
1439     TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1440     histPtPionEtaMcData->Write();
1441
1442   }
1443
1444   //
1445   // Close the file after everthing is done
1446   //
1447   spectraFile.Close();
1448
1449 }
1450
1451
1452
1453