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