| 88b71b9f |
1 | /************************************************************************** |
| 2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
| 3 | * * |
| 4 | * Author: Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de> * |
| 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 | // The task: |
| 17 | // stores TPC PID quantities in a THnSparse |
| 18 | // output can then be used for e.g. dEdx calibration |
| 19 | // |
| 20 | // Author: |
| 21 | // Yvonne Pachmayer <pachmay@physi.uni-heidelberg.de> |
| 22 | // |
| 23 | |
| 24 | |
| 25 | #include "AliTPCcalibResidualPID.h" |
| 26 | #include "TChain.h" |
| 27 | #include "AliAnalysisManager.h" |
| 28 | #include "AliESDEvent.h" |
| 29 | #include "AliAODEvent.h" |
| 30 | #include "AliMCEvent.h" |
| 31 | #include "AliMCParticle.h" |
| 32 | //#include "AliStack.h" |
| 33 | #include "AliPID.h" |
| 34 | #include "AliTPCcalibResidualPID.h" |
| 35 | #include "AliESDtrack.h" |
| 36 | #include "AliESDtrackCuts.h" |
| 37 | #include "AliESDpid.h" |
| 38 | #include "AliESDInputHandler.h" |
| 39 | #include "AliESDv0KineCuts.h" |
| 40 | #include "AliESDv0.h" |
| 41 | #include "AliCentrality.h" |
| 42 | #include "THnSparse.h" |
| 43 | #include "TH2D.h" |
| 44 | #include "TCanvas.h" |
| 45 | #include "TGraphErrors.h" |
| 46 | #include "TMultiGraph.h" |
| 47 | #include "TAxis.h" |
| 48 | #include "TFile.h" |
| 49 | #include "TSpline.h" |
| 50 | #include "TStyle.h" |
| 51 | #include "AliTPCdEdxInfo.h" |
| 52 | #include "TString.h" |
| 53 | #include "TFitResult.h" |
| 54 | #include "TH1F.h" |
| 55 | #include "TLegend.h" |
| 56 | #include "TVirtualFitter.h" |
| 57 | |
| 58 | using namespace std; |
| 59 | |
| 60 | ClassImp(AliTPCcalibResidualPID) |
| 61 | |
| 62 | Double_t AliTPCcalibResidualPID::fgCutGeo = 1.; |
| 63 | Double_t AliTPCcalibResidualPID::fgCutNcr = 0.85; |
| 64 | Double_t AliTPCcalibResidualPID::fgCutNcl = 0.7; |
| 65 | |
| 66 | //________________________________________________________________________ |
| 67 | AliTPCcalibResidualPID::AliTPCcalibResidualPID() |
| 68 | : AliAnalysisTaskSE(), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0), |
| 69 | fUseTPCCutMIGeo(kFALSE), |
| 70 | fUseMCinfo(kTRUE), |
| 71 | fIsPbpOrpPb(kFALSE), |
| 72 | fZvtxCutEvent(9999.0), |
| 73 | fV0KineCuts(0x0), |
| 74 | fNumTagsStored(0), |
| 75 | fV0tags(0x0), |
| 76 | fV0motherIndex(0x0), |
| 77 | fV0motherPDG(0x0), |
| 78 | fProduceAllPadTypes(0), |
| 79 | fProduceGlobal(0), |
| 80 | fProduceShortPads(0), |
| 81 | fProduceMediumPads(0), |
| 82 | fProduceLongPads(0), |
| 83 | fProduceOroc(0), |
| 84 | fHistPidQA(0), |
| 85 | fHistPidQAshort(0), |
| 86 | fHistPidQAmedium(0), |
| 87 | fHistPidQAlong(0), |
| 88 | fHistPidQAoroc(0), |
| 89 | fProduceTPCSignalSparse(0), |
| 90 | fCorrectdEdxEtaDependence(0), |
| 91 | fCorrectdEdxMultiplicityDependence(0), |
| 92 | fThnspTpc(0), |
| 93 | fQAList(0x0), |
| 94 | fhInvMassGamma(0x0), |
| 95 | fhInvMassK0s(0x0), |
| 96 | fhInvMassLambda(0x0), |
| 97 | fhInvMassAntiLambda(0x0), |
| 98 | fhArmenterosAll(0x0), |
| 99 | fhArmenterosGamma(0x0), |
| 100 | fhArmenterosK0s(0x0), |
| 101 | fhArmenterosLambda(0x0), |
| 102 | fhArmenterosAntiLambda(0x0), |
| 103 | fHistSharedClusQAV0Pi(0x0), |
| 104 | fHistSharedClusQAV0Pr(0x0), |
| 105 | fHistSharedClusQAV0El(0x0) |
| 106 | { |
| 107 | // default Constructor |
| 108 | /* fast compilation test |
| 109 | gSystem->Load("libANALYSIS"); |
| 110 | gSystem->Load("libANALYSISalice"); |
| 111 | .L /lustre/alice/akalweit/train/trunk/util/statsQA/AliTPCcalibResidualPID.cxx++ |
| 112 | */ |
| 113 | |
| 114 | } |
| 115 | |
| 116 | |
| 117 | //________________________________________________________________________ |
| 118 | AliTPCcalibResidualPID::AliTPCcalibResidualPID(const char *name) |
| 119 | : AliAnalysisTaskSE(name), fESD(0), fMC(0), fOutputContainer(0), fESDtrackCuts(0), fESDtrackCutsV0(0), fESDpid(0), |
| 120 | fUseTPCCutMIGeo(kFALSE), |
| 121 | fUseMCinfo(kTRUE), |
| 122 | fIsPbpOrpPb(kFALSE), |
| 123 | fZvtxCutEvent(9999.0), |
| 124 | fV0KineCuts(0x0), |
| 125 | fNumTagsStored(0), |
| 126 | fV0tags(0x0), |
| 127 | fV0motherIndex(0x0), |
| 128 | fV0motherPDG(0x0), |
| 129 | fProduceAllPadTypes(0), |
| 130 | fProduceGlobal(0), |
| 131 | fProduceShortPads(0), |
| 132 | fProduceMediumPads(0), |
| 133 | fProduceLongPads(0), |
| 134 | fProduceOroc(0), |
| 135 | fHistPidQA(0), |
| 136 | fHistPidQAshort(0), |
| 137 | fHistPidQAmedium(0), |
| 138 | fHistPidQAlong(0), |
| 139 | fHistPidQAoroc(0), |
| 140 | fProduceTPCSignalSparse(0), |
| 141 | fCorrectdEdxEtaDependence(0), |
| 142 | fCorrectdEdxMultiplicityDependence(0), |
| 143 | fThnspTpc(0), |
| 144 | fQAList(0x0), |
| 145 | fhInvMassGamma(0x0), |
| 146 | fhInvMassK0s(0x0), |
| 147 | fhInvMassLambda(0x0), |
| 148 | fhInvMassAntiLambda(0x0), |
| 149 | fhArmenterosAll(0x0), |
| 150 | fhArmenterosGamma(0x0), |
| 151 | fhArmenterosK0s(0x0), |
| 152 | fhArmenterosLambda(0x0), |
| 153 | fhArmenterosAntiLambda(0x0), |
| 154 | fHistSharedClusQAV0Pi(0x0), |
| 155 | fHistSharedClusQAV0Pr(0x0), |
| 156 | fHistSharedClusQAV0El(0x0) |
| 157 | { |
| 158 | |
| 159 | //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE); |
| 160 | //fESDtrackCutsV0 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE); |
| 161 | |
| 162 | // Constructor |
| 163 | DefineInput(0, TChain::Class()); |
| 164 | DefineOutput(1, TObjArray::Class()); |
| 165 | DefineOutput(2, TObjArray::Class()); |
| 166 | |
| 167 | } |
| 168 | |
| 169 | |
| 170 | //_________________________________________________ |
| 171 | AliTPCcalibResidualPID::~AliTPCcalibResidualPID() |
| 172 | { |
| 173 | delete fOutputContainer; |
| 174 | fOutputContainer = 0; |
| 175 | |
| 176 | delete fQAList; |
| 177 | fQAList = 0; |
| 178 | |
| 179 | delete fESDtrackCuts; |
| 180 | fESDtrackCuts = 0; |
| 181 | |
| 182 | delete fESDtrackCutsV0; |
| 183 | fESDtrackCutsV0 = 0; |
| 184 | |
| 185 | delete fV0KineCuts; |
| 186 | fV0KineCuts = 0; |
| 187 | |
| 188 | delete fV0tags; |
| 189 | fV0tags = 0; |
| 190 | fNumTagsStored = 0; |
| 191 | |
| 192 | delete fV0motherIndex; |
| 193 | fV0motherIndex = 0; |
| 194 | |
| 195 | delete fV0motherPDG; |
| 196 | fV0motherPDG = 0; |
| 197 | } |
| 198 | |
| 199 | |
| 200 | //________________________________________________________________________ |
| 201 | void AliTPCcalibResidualPID::UserCreateOutputObjects() |
| 202 | { |
| 203 | |
| 204 | // THnSparse binning |
| 205 | const Int_t kNdim = 9; |
| 206 | //v0id, dEdx, TPCsigele, TPCsigpion, TOFbit, eta, TPCclus, centr, p |
| 207 | Int_t bins[kNdim] = { 4, 250, 200, 200, 8, 20, 50, 20, 100}; |
| 208 | Double_t xmin[kNdim] = { 0.5, 30, -10., -10., -1.5, -1., 60., 0., 0.1}; |
| 209 | Double_t xmax[kNdim] = { 4.5, 500., 10., 10., 6.5, 1., 160., 100, 4}; |
| 210 | fThnspTpc= new THnSparseF("tpcsignals", "TPC signal;v0id;tpc signal;tpcsigele;tpcsigpion;tofbit;eta;tpcclus;centr;p (GeV/c)", kNdim, bins, xmin, xmax); |
| 211 | BinLogAxis(fThnspTpc, 8); |
| 212 | |
| 213 | // |
| 214 | // 0.ptot, 1.tpcSig, 2.particle ID, 3. assumed particle, 4. nSigmaTPC (4x), 5. nSigmaTOF (4x), 6. centrality |
| 215 | // Concerning 2 (particle ID): |
| 216 | // - in case of MC: Contains MC ID. Bin 1-4 (<=> Slot 0-3): el, pi, ka, pr |
| 217 | // - in case of data: Contains V0 particles + bin with all others: Bin 1-4 (<=> Slot 0-3): All non-V0s, V0-el, V0-pi, V0-pr |
| 218 | // |
| 219 | |
| 220 | Int_t binsHistQA[7] = {135, 1980, 4, 5, 40, 10, 40}; |
| 221 | Double_t xminHistQA[7] = {0.1, 20, -0.5, -0.5, -10, -5, 0.}; |
| 222 | Double_t xmaxHistQA[7] = {50., 2000, 3.5, 4.5, 10, 5, 20000}; |
| 223 | fHistPidQA = new THnSparseF("fHistPidQA","PID QA",7,binsHistQA,xminHistQA,xmaxHistQA); |
| 224 | BinLogAxis(fHistPidQA, 0); |
| 225 | |
| 226 | // |
| 227 | fHistPidQAshort = new THnSparseF("fHistPidQAshort" ,"PID QA -- short pads",7,binsHistQA,xminHistQA,xmaxHistQA); |
| 228 | fHistPidQAmedium = new THnSparseF("fHistPidQAmedium","PID QA -- med pads",7,binsHistQA,xminHistQA,xmaxHistQA); |
| 229 | fHistPidQAlong = new THnSparseF("fHistPidQAlong" ,"PID QA -- long pads",7,binsHistQA,xminHistQA,xmaxHistQA); |
| 230 | fHistPidQAoroc = new THnSparseF("fHistPidQAoroc" ,"PID QA -- oroc",7,binsHistQA,xminHistQA,xmaxHistQA); |
| 231 | BinLogAxis(fHistPidQAshort, 0); |
| 232 | BinLogAxis(fHistPidQAmedium, 0); |
| 233 | BinLogAxis(fHistPidQAlong, 0); |
| 234 | BinLogAxis(fHistPidQAoroc, 0); |
| 235 | // |
| 236 | fOutputContainer = new TObjArray(2); |
| 237 | fOutputContainer->SetName(GetName()); |
| 238 | fOutputContainer->SetOwner(); |
| 239 | |
| 240 | if(fProduceTPCSignalSparse)fOutputContainer->Add(fThnspTpc); |
| 241 | if(fProduceGlobal)fOutputContainer->Add(fHistPidQA); |
| 242 | if(fProduceAllPadTypes || fProduceShortPads)fOutputContainer->Add(fHistPidQAshort); |
| 243 | if(fProduceAllPadTypes || fProduceMediumPads)fOutputContainer->Add(fHistPidQAmedium); |
| 244 | if(fProduceAllPadTypes || fProduceLongPads)fOutputContainer->Add(fHistPidQAlong); |
| 245 | if(fProduceAllPadTypes || fProduceOroc)fOutputContainer->Add(fHistPidQAoroc); |
| 246 | |
| 247 | // V0 Kine cuts |
| 248 | fV0KineCuts = new AliESDv0KineCuts; |
| 249 | fV0KineCuts->SetGammaCutChi2NDF(5.); |
| 250 | // Only accept V0el with prod. radius within 45 cm -> PID will by systematically biased for larger values! |
| 251 | Float_t gammaProdVertexRadiusCuts[2] = { 3.0, 45. }; |
| 252 | fV0KineCuts->SetGammaCutVertexR(&gammaProdVertexRadiusCuts[0]); |
| 253 | |
| 254 | |
| 255 | fQAList = new TObjArray(4); |
| 256 | fQAList->SetName(GetName()); |
| 257 | fQAList->SetOwner(); |
| 258 | |
| 259 | fhInvMassGamma = new TH1F("fhInvMassGamma", "Invariant Mass of gammas; m_{ee} (GeV/#it{c}^{2}); Entries", 200, 0., 0.2); |
| 260 | fhInvMassK0s = new TH1F("fhInvMassK0s", "Invariant Mass of K0s; m_{#pi#pi} (GeV/#it{c}^{2}); Entries;", 200, 0.45, 0.55); |
| 261 | fhInvMassLambda = new TH1F("fhInvMassLambda", "Invariant Mass of lambdas; m_{p#pi^{-}} (GeV/#it{c}^{2}); Entries", 200, 1.05, 1.15); |
| 262 | fhInvMassAntiLambda = new TH1F("fhInvMassAntiLambda", "Invariant Mass of anti-lambdas; m_{#pi^{+}#bar{p}} (GeV/#it{c}^{2}); Entries", |
| 263 | 200, 1.05, 1.15); |
| 264 | |
| 265 | fQAList->Add(fhInvMassGamma); |
| 266 | fQAList->Add(fhInvMassK0s); |
| 267 | fQAList->Add(fhInvMassLambda); |
| 268 | fQAList->Add(fhInvMassAntiLambda); |
| 269 | |
| 270 | |
| 271 | fhArmenterosAll = new TH2F("fhArmenterosAll", |
| 272 | "Armenteros plot all V0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", |
| 273 | 200, -1., 1., 200, 0., 0.4); |
| 274 | fhArmenterosGamma = new TH2F("fhArmenterosGamma", |
| 275 | "Armenteros plot Gamma;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", |
| 276 | 200, -1., 1., 200, 0., 0.4); |
| 277 | fhArmenterosK0s = new TH2F("fhArmenterosK0s", |
| 278 | "Armenteros plot K0s;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", |
| 279 | 200, -1., 1., 200, 0., 0.4); |
| 280 | fhArmenterosLambda = new TH2F("fhArmenterosLambda", |
| 281 | "Armenteros plot lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", |
| 282 | 200, -1., 1., 200, 0., 0.4); |
| 283 | fhArmenterosAntiLambda = new TH2F("fhArmenterosAntiLambda", |
| 284 | "Armenteros plot anti-lambda;#alpha = (#it{p}^{+}_{L}-#it{p}^{-}_{L})/(#it{p}^{+}_{L}+#it{p}^{-}_{L});#it{q}_{T} (GeV/#it{c})", |
| 285 | 200, -1., 1., 200, 0., 0.4); |
| 286 | |
| 287 | fQAList->Add(fhArmenterosAll); |
| 288 | fQAList->Add(fhArmenterosGamma); |
| 289 | fQAList->Add(fhArmenterosK0s); |
| 290 | fQAList->Add(fhArmenterosLambda); |
| 291 | fQAList->Add(fhArmenterosAntiLambda); |
| 292 | |
| 293 | |
| 294 | const Int_t dimQASharedClusters = 4; |
| 295 | const TString axisTitles[dimQASharedClusters] = { "#it{p}_{TPC} (GeV/#it{c})", "#it{#Delta'}", "#it{N}_{shared cl}", "Pad row"}; |
| 296 | Int_t binsHistQASharedClusters[dimQASharedClusters] = { 100, 100, 160, 160}; |
| 297 | Double_t xminHistQASharedClusters[dimQASharedClusters] = { 0.3, 0.5, 0, -1}; |
| 298 | Double_t xmaxHistQASharedClusters[dimQASharedClusters] = { 20, 1.5, 160, 159}; |
| 299 | |
| 300 | fHistSharedClusQAV0Pi = new THnSparseF("fHistSharedClusQAV0Pi" ,"PID QA shared clusters - V0 pi", dimQASharedClusters, |
| 301 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); |
| 302 | BinLogAxis(fHistSharedClusQAV0Pi, 0); |
| 303 | for (Int_t i = 0; i < dimQASharedClusters; i++) |
| 304 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); |
| 305 | fQAList->Add(fHistSharedClusQAV0Pi); |
| 306 | |
| 307 | fHistSharedClusQAV0Pr = new THnSparseF("fHistSharedClusQAV0Pr" ,"PID QA shared clusters - V0 pr", dimQASharedClusters, |
| 308 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); |
| 309 | BinLogAxis(fHistSharedClusQAV0Pr, 0); |
| 310 | for (Int_t i = 0; i < dimQASharedClusters; i++) |
| 311 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); |
| 312 | fQAList->Add(fHistSharedClusQAV0Pr); |
| 313 | |
| 314 | fHistSharedClusQAV0El = new THnSparseF("fHistSharedClusQAV0El" ,"PID QA shared clusters - V0 el", dimQASharedClusters, |
| 315 | binsHistQASharedClusters, xminHistQASharedClusters, xmaxHistQASharedClusters); |
| 316 | BinLogAxis(fHistSharedClusQAV0El, 0); |
| 317 | for (Int_t i = 0; i < dimQASharedClusters; i++) |
| 318 | fHistSharedClusQAV0Pi->GetAxis(i)->SetTitle(axisTitles[i].Data()); |
| 319 | fQAList->Add(fHistSharedClusQAV0El); |
| 320 | |
| 321 | PostData(1,fOutputContainer); |
| 322 | PostData(2,fQAList); |
| 323 | } |
| 324 | |
| 325 | |
| 326 | //_____________________________________________________________________________ |
| 327 | void AliTPCcalibResidualPID::UserExec(Option_t *) |
| 328 | { |
| 329 | //calls the Process function |
| 330 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
| 331 | |
| 332 | if (!esdH) { |
| 333 | printf("ERROR: Could not get ESDInputHandler \n"); |
| 334 | } |
| 335 | else fESD = esdH->GetEvent(); |
| 336 | |
| 337 | // If MC, forward MCevent |
| 338 | fMC = dynamic_cast<AliMCEvent*>(MCEvent()); |
| 339 | // |
| 340 | Process(fESD, fMC); |
| 341 | // |
| 342 | PostData(1,fOutputContainer); |
| 343 | PostData(2,fQAList); |
| 344 | } |
| 345 | |
| 346 | |
| 347 | |
| 348 | //________________________________________________________________________ |
| 349 | void AliTPCcalibResidualPID::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent) |
| 350 | { |
| 351 | |
| 352 | //called for each event |
| 353 | if (!esdEvent) { |
| 354 | Printf("ERROR: esdEvent not available"); |
| 355 | return; |
| 356 | } |
| 357 | |
| 358 | if(!((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))) printf("ESD inputhandler not available \n"); |
| 359 | |
| 360 | if (!fESDpid) { |
| 361 | fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid(); |
| 362 | } |
| 363 | |
| 364 | if (!fESDpid || !fV0KineCuts) |
| 365 | return; |
| 366 | |
| 367 | |
| 368 | Float_t centralityFper=99; |
| 369 | |
| 370 | AliCentrality *esdCentrality = esdEvent->GetCentrality(); |
| 371 | centralityFper = esdCentrality->GetCentralityPercentile("V0M"); |
| 372 | |
| 373 | if (!GetVertexIsOk(esdEvent)) |
| 374 | return; |
| 375 | |
| 376 | const AliESDVertex* fESDvertex = esdEvent->GetPrimaryVertexTracks(); |
| 377 | if (!fESDvertex) |
| 378 | return; |
| 379 | |
| 380 | Int_t ncontr = fESDvertex->GetNContributors(); |
| 381 | |
| 382 | if (ncontr <= 0) |
| 383 | return; |
| 384 | |
| 385 | // Fill V0 arrays with V0s |
| 386 | FillV0PIDlist(esdEvent); |
| 387 | |
| 388 | // Array with flags wheter QA for this V0 was already done or not |
| 389 | const Int_t numV0s = esdEvent->GetNumberOfV0s(); |
| 390 | Bool_t v0QAadded[numV0s]; |
| 391 | for (Int_t i = 0; i < numV0s; i++) |
| 392 | v0QAadded[i] = kFALSE; |
| 393 | |
| 394 | Int_t nTotTracks = esdEvent->GetNumberOfTracks(); |
| 395 | for (Int_t iTracks = 0; iTracks < nTotTracks; iTracks++){//begin track loop |
| 396 | AliESDtrack *trackESD = esdEvent->GetTrack(iTracks); |
| 397 | if(!trackESD) { |
| 398 | Printf("ERROR: Could not receive track %d (esd loop)", iTracks); |
| 399 | continue; |
| 400 | } |
| 401 | if((TMath::Abs(trackESD->Eta())) > 0.9) |
| 402 | continue; |
| 403 | |
| 404 | // Do not distinguish positively and negatively charged V0's |
| 405 | Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks)); |
| 406 | if (v0tagAllCharges == -99) { |
| 407 | AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, esdEvent->GetNumberOfTracks(), |
| 408 | fV0tags != 0x0)); |
| 409 | continue; |
| 410 | } |
| 411 | |
| 412 | Bool_t isV0el = v0tagAllCharges == 14; |
| 413 | Bool_t isV0pi = v0tagAllCharges == 15; |
| 414 | Bool_t isV0pr = v0tagAllCharges == 16; |
| 415 | Bool_t isV0 = isV0el || isV0pi || isV0pr; |
| 416 | |
| 417 | if (mcEvent && fUseMCinfo) { |
| 418 | // For MC, do not look for V0s, i.e. demand the non-V0 track cuts |
| 419 | if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue; |
| 420 | |
| 421 | if (fUseTPCCutMIGeo) { |
| 422 | // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included. |
| 423 | if (!TPCCutMIGeo(trackESD, esdEvent)) |
| 424 | continue; |
| 425 | } |
| 426 | else { |
| 427 | // If cut on geometry is not active, always cut on num clusters |
| 428 | if (trackESD->GetTPCsignalN() < 60) |
| 429 | continue; |
| 430 | } |
| 431 | } |
| 432 | else { |
| 433 | // For data, take V0 AND non-V0 with separate cuts |
| 434 | //if (!isV0 && fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) continue; |
| 435 | if (isV0) { |
| 436 | if (fESDtrackCutsV0 && !fESDtrackCutsV0->AcceptTrack(trackESD)) |
| 437 | continue; |
| 438 | } |
| 439 | else { |
| 440 | if (fESDtrackCuts && !fESDtrackCuts->AcceptTrack(trackESD)) |
| 441 | continue; |
| 442 | } |
| 443 | |
| 444 | if (fUseTPCCutMIGeo) { |
| 445 | // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included. |
| 446 | if (!isV0) { |
| 447 | if (!TPCCutMIGeo(trackESD, esdEvent)) |
| 448 | continue; |
| 449 | } |
| 450 | else { |
| 451 | // One should not cut on geometry for V0's since they have different topology. (Loosely) cut on num of clusters instead. |
| 452 | if (trackESD->GetTPCsignalN() < 60) |
| 453 | continue; |
| 454 | } |
| 455 | } |
| 456 | else { |
| 457 | // If cut on geometry is not active, always cut on num clusters |
| 458 | if (trackESD->GetTPCsignalN() < 60) |
| 459 | continue; |
| 460 | } |
| 461 | } |
| 462 | |
| 463 | const AliExternalTrackParam *paramIn = trackESD->GetInnerParam(); |
| 464 | Float_t precin=-1; |
| 465 | if(paramIn) { |
| 466 | precin=paramIn->GetP(); |
| 467 | } else { |
| 468 | continue; |
| 469 | } |
| 470 | Int_t precdefault=CompareFloat(precin,-1); |
| 471 | if(precdefault==1) continue; // momentum cut |
| 472 | // |
| 473 | |
| 474 | AliMCParticle* mcTrack = 0x0; |
| 475 | Int_t particleID = -1; |
| 476 | Int_t v0id = 0; |
| 477 | |
| 478 | if (mcEvent && fUseMCinfo) { |
| 479 | // MC - particle ID = MC ID |
| 480 | Int_t label = trackESD->GetLabel(); |
| 481 | |
| 482 | if (label < 0) |
| 483 | continue; // If MC is available, reject tracks with negative ESD label |
| 484 | |
| 485 | mcTrack = dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(TMath::Abs(label))); |
| 486 | if (!mcTrack) { |
| 487 | Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks); |
| 488 | continue; |
| 489 | } |
| 490 | |
| 491 | /*// Only accept MC primaries |
| 492 | if (!mcEvent->Stack()->IsPhysicalPrimary(TMath::Abs(label))) { |
| 493 | continue; |
| 494 | }*/ |
| 495 | |
| 496 | Int_t pdgAbs = TMath::Abs(mcTrack->PdgCode()); |
| 497 | |
| 498 | if (pdgAbs == 11) { // electron |
| 499 | particleID = 0; |
| 500 | } |
| 501 | else if (pdgAbs == 211) { // pion |
| 502 | particleID = 1; |
| 503 | } |
| 504 | else if (pdgAbs == 321) { // kaon |
| 505 | particleID = 2; |
| 506 | } |
| 507 | else if (pdgAbs == 2212) { // proton |
| 508 | particleID = 3; |
| 509 | } |
| 510 | else |
| 511 | continue; // Reject all other particles |
| 512 | } |
| 513 | else { |
| 514 | // Data - particle ID = V0 ID (if V0) |
| 515 | if (isV0pi) { // pion |
| 516 | particleID = 2; |
| 517 | v0id = 2; |
| 518 | } |
| 519 | else if (isV0el) { // electron |
| 520 | particleID = 1; |
| 521 | v0id = 1; |
| 522 | } |
| 523 | else if (isV0pr) { // proton |
| 524 | particleID = 3; |
| 525 | v0id = 3; |
| 526 | } |
| 527 | else { // No V0-ID available (species must be determined via TPC, TOF, ...) |
| 528 | particleID = 0; |
| 529 | } |
| 530 | } |
| 531 | |
| 532 | |
| 533 | Float_t tpcsignal=trackESD->GetTPCsignal(); |
| 534 | Int_t tofbit=-1; |
| 535 | Int_t iTOFpid=0; |
| 536 | Int_t ikTIME=0; |
| 537 | Double_t tpcnsigmaele=-10; |
| 538 | Double_t tpcnsigmapion=-10; |
| 539 | // |
| 540 | if ((trackESD->GetStatus() & AliESDtrack::kTOFpid)) iTOFpid = 1; |
| 541 | if ((trackESD->GetStatus() & AliESDtrack::kTIME)) ikTIME = 1; |
| 542 | Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero(); |
| 543 | // |
| 544 | if((iTOFpid==1) &&(ikTIME==1)){ |
| 545 | Double_t tofnsigmaele= fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kElectron, time0); |
| 546 | Double_t tofnsigmapion=fESDpid->NumberOfSigmasTOF(trackESD,AliPID::kPion, time0); |
| 547 | if (TMath::Abs(tofnsigmapion)<3) tofbit = 1; |
| 548 | if (TMath::Abs(tofnsigmaele)<3) tofbit = 0; |
| 549 | tpcnsigmaele=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kElectron); |
| 550 | tpcnsigmapion=fESDpid->NumberOfSigmasTPC(trackESD,AliPID::kPion); |
| 551 | } |
| 552 | // |
| 553 | Int_t tpcnclusN=trackESD->GetTPCsignalN(); |
| 554 | Double_t eta=trackESD->Eta(); |
| 555 | |
| 556 | // |
| 557 | if(fProduceTPCSignalSparse){ |
| 558 | Double_t contentSignal[9]; |
| 559 | contentSignal[0]=v0id; |
| 560 | contentSignal[1]=tpcsignal; |
| 561 | contentSignal[2]=tpcnsigmaele; |
| 562 | contentSignal[3]=tpcnsigmapion; |
| 563 | contentSignal[4]=tofbit; |
| 564 | contentSignal[5]=eta; |
| 565 | contentSignal[6]=tpcnclusN; |
| 566 | contentSignal[7]=centralityFper; |
| 567 | contentSignal[8]=precin; |
| 568 | // |
| 569 | if (fThnspTpc->GetEntries() < 1e6) fThnspTpc->Fill(contentSignal); |
| 570 | }//should it be created or not |
| 571 | |
| 572 | |
| 573 | |
| 574 | // |
| 575 | // 2nd THnSparse |
| 576 | // |
| 577 | Double_t tpcQA[5] = {fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kElectron), |
| 578 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kPion), |
| 579 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kKaon), |
| 580 | fESDpid->NumberOfSigmasTPC(trackESD, AliPID::kProton), |
| 581 | 0}; |
| 582 | Double_t tofQA[5] = {fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kElectron, time0), |
| 583 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kPion, time0), |
| 584 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kKaon, time0), |
| 585 | fESDpid->NumberOfSigmasTOF(trackESD, AliPID::kProton, time0), |
| 586 | 0 }; |
| 587 | |
| 588 | // id 5 is just again Kaons in restricted eta range |
| 589 | tpcQA[4] = tpcQA[2]; |
| 590 | tofQA[4] = tofQA[2]; |
| 591 | |
| 592 | // |
| 593 | // dE/dx in different pad regions |
| 594 | // |
| 595 | AliTPCdEdxInfo * infoTpcPid = trackESD->GetTPCdEdxInfo(); |
| 596 | Double32_t signal[4]; Char_t ncl[3]; Char_t nrows[3]; |
| 597 | if (infoTpcPid) { |
| 598 | infoTpcPid->GetTPCSignalRegionInfo(signal, ncl, nrows); |
| 599 | } else { |
| 600 | for(Int_t iarr = 0; iarr < 3; iarr++) { |
| 601 | signal[iarr] = 0; |
| 602 | ncl[iarr] = 0; |
| 603 | nrows[iarr] = 0; |
| 604 | } |
| 605 | signal[3] = 0; |
| 606 | } |
| 607 | // |
| 608 | UInt_t status = trackESD->GetStatus(); |
| 609 | Bool_t hasTOFout = status&AliESDtrack::kTOFout; |
| 610 | Bool_t hasTOFtime = status&AliESDtrack::kTIME; |
| 611 | Bool_t hasTOF = kFALSE; |
| 612 | if (hasTOFout && hasTOFtime) hasTOF = kTRUE; |
| 613 | Float_t length = trackESD->GetIntegratedLength(); |
| 614 | // Length check only for primaries! |
| 615 | if (length < 350 && !isV0) hasTOF = kFALSE; |
| 616 | // |
| 617 | |
| 618 | if (!hasTOF) { |
| 619 | // Make sure that number of sigmas is large in this case, so that the track will be rejected if a TOF cut is applied |
| 620 | for (Int_t i = 0; i < 5; i++) { |
| 621 | tofQA[i] = 999; |
| 622 | } |
| 623 | } |
| 624 | |
| 625 | |
| 626 | Double_t processedTPCsignal[5] = { tpcsignal, tpcsignal, tpcsignal, tpcsignal, tpcsignal }; |
| 627 | |
| 628 | |
| 629 | if (fCorrectdEdxEtaDependence && fCorrectdEdxMultiplicityDependence) { |
| 630 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron, |
| 631 | AliTPCPIDResponse::kdEdxDefault); |
| 632 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion, |
| 633 | AliTPCPIDResponse::kdEdxDefault); |
| 634 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon, |
| 635 | AliTPCPIDResponse::kdEdxDefault); |
| 636 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaAndMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton, |
| 637 | AliTPCPIDResponse::kdEdxDefault); |
| 638 | } |
| 639 | else if (fCorrectdEdxEtaDependence) { |
| 640 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault); |
| 641 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault); |
| 642 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault); |
| 643 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetEtaCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault); |
| 644 | } |
| 645 | else if (fCorrectdEdxMultiplicityDependence) { |
| 646 | processedTPCsignal[0] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault); |
| 647 | processedTPCsignal[1] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault); |
| 648 | processedTPCsignal[2] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault); |
| 649 | processedTPCsignal[3] = fESDpid->GetTPCResponse().GetMultiplicityCorrectedTrackdEdx(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault); |
| 650 | } |
| 651 | |
| 652 | // id 5 is just again Kaons in restricted eta range |
| 653 | processedTPCsignal[4] = processedTPCsignal[2]; |
| 654 | |
| 655 | for(Int_t iPart = 0; iPart < 5; iPart++) { |
| 656 | // Only accept "Kaons" within |eta| < 0.2 for index 4 in case of data (no contamination in case of MC, so this index is not used) |
| 657 | if (iPart == 4 && ((mcEvent && fUseMCinfo) || abs(trackESD->Eta()) > 0.2)) |
| 658 | continue; |
| 659 | |
| 660 | Double_t vecHistQA[7] = {precin, processedTPCsignal[iPart], particleID, iPart, tpcQA[iPart], tofQA[iPart], nTotTracks}; |
| 661 | if (fProduceGlobal) fHistPidQA->Fill(vecHistQA); |
| 662 | vecHistQA[1] = signal[0]; vecHistQA[4] = ncl[0]; |
| 663 | if ((fProduceAllPadTypes || fProduceShortPads)) fHistPidQAshort->Fill(vecHistQA); |
| 664 | // |
| 665 | vecHistQA[1] = signal[1]; vecHistQA[4] = ncl[1]; |
| 666 | if ((fProduceAllPadTypes || fProduceMediumPads)) fHistPidQAmedium->Fill(vecHistQA); |
| 667 | // |
| 668 | vecHistQA[1] = signal[2]; vecHistQA[4] = ncl[2]; |
| 669 | if ((fProduceAllPadTypes || fProduceLongPads)) fHistPidQAlong->Fill(vecHistQA); |
| 670 | // |
| 671 | vecHistQA[1] = signal[3]; vecHistQA[4] = nrows[1] + nrows[2]; |
| 672 | if ((fProduceAllPadTypes || fProduceOroc)) fHistPidQAoroc->Fill(vecHistQA); |
| 673 | } |
| 674 | |
| 675 | if (isV0) { |
| 676 | // Fill QA hists for shared clusters dE/dx |
| 677 | Double_t vecHistQA[4] = { precin, -1, trackESD->GetTPCSharedMap().CountBits(), -1 }; |
| 678 | THnSparse* currHist = fHistSharedClusQAV0Pi; |
| 679 | |
| 680 | if (isV0pi) { |
| 681 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, |
| 682 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); |
| 683 | if (expectedDeDx > 0 ) { |
| 684 | vecHistQA[1] = tpcsignal / expectedDeDx; |
| 685 | currHist = fHistSharedClusQAV0Pi; |
| 686 | } |
| 687 | } |
| 688 | else if (isV0pr) { |
| 689 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, |
| 690 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); |
| 691 | if (expectedDeDx > 0 ) { |
| 692 | vecHistQA[1] = tpcsignal / expectedDeDx; |
| 693 | currHist = fHistSharedClusQAV0Pr; |
| 694 | } |
| 695 | } |
| 696 | else if (isV0el) { |
| 697 | Double_t expectedDeDx = fESDpid->GetTPCResponse().GetExpectedSignal(trackESD, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, |
| 698 | fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence); |
| 699 | if (expectedDeDx > 0 ) { |
| 700 | vecHistQA[1] = tpcsignal / expectedDeDx; |
| 701 | currHist = fHistSharedClusQAV0El; |
| 702 | } |
| 703 | } |
| 704 | |
| 705 | Int_t iRowInd = -1; |
| 706 | // iRowInd == -1 for "all rows w/o multiple counting" |
| 707 | currHist->Fill(vecHistQA); |
| 708 | |
| 709 | // Fill hist for every pad row with shared cluster |
| 710 | for (iRowInd = 0; iRowInd < 159; iRowInd++) { |
| 711 | if (trackESD->GetTPCSharedMap().TestBitNumber(iRowInd)) { |
| 712 | vecHistQA[3] = iRowInd; |
| 713 | currHist->Fill(vecHistQA); |
| 714 | } |
| 715 | } |
| 716 | |
| 717 | |
| 718 | |
| 719 | // Check, whether the QA for this V0 mother was already filled into the histograms (is the case if the other daughter has already |
| 720 | // been processed). If not, set flag to kTRUE and fill the QA. |
| 721 | const Int_t iV0 = GetV0motherIndex(iTracks); |
| 722 | if (!v0QAadded[iV0]) { |
| 723 | v0QAadded[iV0] = kTRUE; |
| 724 | |
| 725 | AliESDv0* esdV0 = (AliESDv0*)esdEvent->GetV0(iV0); |
| 726 | |
| 727 | if (!esdV0) { |
| 728 | printf("Error: V0 tagged, but does not exist!\n"); |
| 729 | continue; |
| 730 | } |
| 731 | |
| 732 | Float_t armVar[2] = {0.0,0.0}; |
| 733 | fV0KineCuts->Armenteros(esdV0, armVar); |
| 734 | |
| 735 | const Int_t motherPDG = GetV0motherPDG(iTracks); |
| 736 | |
| 737 | // Mother and daughter match the requirements, otherwise it wouldn't be tagged as a V0. Now just fill the QA histos. |
| 738 | fhArmenterosAll->Fill(armVar[0], armVar[1]); |
| 739 | //if (esdV0->TestBit(BIT(14))) { |
| 740 | if (TMath::Abs(motherPDG) == 22) { |
| 741 | fhInvMassGamma->Fill(esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron)); |
| 742 | fhArmenterosGamma->Fill(armVar[0], armVar[1]); |
| 743 | } |
| 744 | else if (TMath::Abs(motherPDG) == 310) { |
| 745 | //else if (esdV0->TestBit(BIT(15))) { |
| 746 | fhInvMassK0s->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kPion)); |
| 747 | fhArmenterosK0s->Fill(armVar[0], armVar[1]); |
| 748 | } |
| 749 | else if (motherPDG == 3122) { |
| 750 | //else if (esdV0->TestBit(BIT(16))) { |
| 751 | fhInvMassLambda->Fill(esdV0->GetEffMass(AliPID::kProton, AliPID::kPion)); |
| 752 | fhArmenterosLambda->Fill(armVar[0], armVar[1]); |
| 753 | } |
| 754 | else if (motherPDG == -3122) { |
| 755 | //else if (esdV0->TestBit(BIT(17))) { |
| 756 | fhInvMassAntiLambda->Fill(esdV0->GetEffMass(AliPID::kPion, AliPID::kProton)); |
| 757 | fhArmenterosAntiLambda->Fill(armVar[0], armVar[1]); |
| 758 | } |
| 759 | } |
| 760 | } |
| 761 | } //end track loop |
| 762 | |
| 763 | |
| 764 | // Clear the V0 PID arrays |
| 765 | ClearV0PIDlist(); |
| 766 | } |
| 767 | |
| 768 | |
| 769 | //________________________________________________________________________ |
| 770 | Int_t AliTPCcalibResidualPID::CompareFloat(Float_t f1, Float_t f2) const |
| 771 | { |
| 772 | //compares if the Float_t f1 is equal with f2 and returns 1 if true and 0 if false |
| 773 | Float_t precision = 0.00001; |
| 774 | if (((f1 - precision) < f2) && |
| 775 | ((f1 + precision) > f2)) |
| 776 | { |
| 777 | return 1; |
| 778 | } |
| 779 | else |
| 780 | { |
| 781 | return 0; |
| 782 | } |
| 783 | } |
| 784 | |
| 785 | |
| 786 | //________________________________________________________________________ |
| 787 | void AliTPCcalibResidualPID::Terminate(const Option_t *) |
| 788 | { |
| 789 | |
| 790 | |
| 791 | } |
| 792 | |
| 793 | |
| 794 | |
| 795 | //________________________________________________________________________ |
| 796 | void AliTPCcalibResidualPID::BinLogAxis(const THnSparseF *h, Int_t axisNumber) { |
| 797 | // |
| 798 | // Method for the correct logarithmic binning of histograms |
| 799 | // |
| 800 | TAxis *axis = h->GetAxis(axisNumber); |
| 801 | int bins = axis->GetNbins(); |
| 802 | |
| 803 | Double_t from = axis->GetXmin(); |
| 804 | Double_t to = axis->GetXmax(); |
| 805 | Double_t *newBins = new Double_t[bins + 1]; |
| 806 | |
| 807 | newBins[0] = from; |
| 808 | Double_t factor = pow(to/from, 1./bins); |
| 809 | |
| 810 | for (int i = 1; i <= bins; i++) { |
| 811 | newBins[i] = factor * newBins[i-1]; |
| 812 | } |
| 813 | axis->Set(bins, newBins); |
| 814 | delete [] newBins; |
| 815 | |
| 816 | } |
| 817 | |
| 818 | |
| 819 | //________________________________________________________________________ |
| 820 | Double_t* AliTPCcalibResidualPID::ExtractResidualPID(THnSparseF * histPidQA, const Bool_t useV0s, const Char_t * outFile, |
| 821 | const Char_t * type, const Char_t * period, const Char_t * pass, |
| 822 | const Char_t * system, const Double_t * initialParameters, |
| 823 | const Char_t *dedxtype, |
| 824 | AliTPCcalibResidualPID::FitType fitType) { |
| 825 | // |
| 826 | // (1.) obtain residual graphs |
| 827 | // |
| 828 | TObjArray * arr = 0x0; |
| 829 | |
| 830 | Bool_t isMC = kFALSE; |
| 831 | if (strcmp(type, "MC") == 0) { |
| 832 | arr = GetResidualGraphsMC(histPidQA, system); |
| 833 | isMC = kTRUE; |
| 834 | } |
| 835 | else if (strcmp(type, "DATA") == 0) { |
| 836 | arr = GetResidualGraphs(histPidQA, system, useV0s); |
| 837 | } |
| 838 | else { |
| 839 | Printf("ERROR - ExtractResidualPID: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type); |
| 840 | |
| 841 | return 0x0; |
| 842 | } |
| 843 | |
| 844 | Bool_t isPPb = kFALSE; |
| 845 | if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { |
| 846 | Printf("p-Pb/Pb-p detected - Applying special handling!"); |
| 847 | isPPb = kTRUE; |
| 848 | } |
| 849 | // |
| 850 | // (2.) get old-style Bethe-Bloch parameters |
| 851 | // |
| 852 | TF1* parametrisation = FitBB(arr, isMC, isPPb, useV0s, initialParameters, fitType); |
| 853 | // |
| 854 | // (3.) obtain response functions to OADB |
| 855 | // |
| 856 | TObjArray* arrResponse = GetResponseFunctions(parametrisation, arr, type, period, pass, system, dedxtype); |
| 857 | // |
| 858 | // (4.) write results to file and exit |
| 859 | // |
| 860 | if (arrResponse) { |
| 861 | TFile outputFile(outFile,"RECREATE"); |
| 862 | arrResponse->Write(); |
| 863 | outputFile.Close(); |
| 864 | } |
| 865 | // |
| 866 | |
| 867 | return parametrisation->GetParameters(); |
| 868 | } |
| 869 | |
| 870 | |
| 871 | //________________________________________________________________________ |
| 872 | TObjArray * AliTPCcalibResidualPID::GetSeparation(THnSparseF * histPidQA, Int_t kParticle1, Int_t kParticle2) { |
| 873 | |
| 874 | //get THnSparse |
| 875 | THnSparse * hist = histPidQA; |
| 876 | |
| 877 | |
| 878 | Float_t nSigmaPlus1 = 0, nSigmaMinus1 = 0; //sigma of first particle in the TPC |
| 879 | Float_t nSigmaPlus2 = 0, nSigmaMinus2 = 0; //sigma of second particle in the TPC |
| 880 | |
| 881 | if(kParticle1 == 0){ // electron |
| 882 | nSigmaMinus1 = -1.999; |
| 883 | nSigmaPlus1 = 2.999; |
| 884 | } else if(kParticle1 == 1){ //pion |
| 885 | nSigmaMinus1 = -1.999; |
| 886 | nSigmaPlus1 = 2.999; |
| 887 | } else if(kParticle1 == 2){ //kaons |
| 888 | nSigmaMinus1 = -2.999; |
| 889 | nSigmaPlus1 = 2.999; |
| 890 | } else if(kParticle1 == 3){ //protons |
| 891 | nSigmaMinus1 = -2.999; |
| 892 | nSigmaPlus1 = 2.999; |
| 893 | } |
| 894 | |
| 895 | // |
| 896 | // 1. select and fit first particle |
| 897 | // |
| 898 | hist->GetAxis(3)->SetRangeUser(kParticle1,kParticle1); |
| 899 | hist->GetAxis(4)->SetRangeUser(nSigmaMinus1,nSigmaPlus1); |
| 900 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 901 | TH2D * histPart1 = (TH2D*) hist->Projection(1,0)->Clone("histPart1"); |
| 902 | histPart1->SetTitle(";p /(GeV/c) ; TPCSignal /(a.u.)"); |
| 903 | histPart1->RebinX(4); |
| 904 | |
| 905 | TObjArray arr; |
| 906 | histPart1->FitSlicesY(0,0,-1,10,"QNR",&arr); |
| 907 | TH1D * PartMean1 = (TH1D *) arr.At(1)->Clone("PartMean1"); |
| 908 | TH1D * PartSigma1 = (TH1D *) arr.At(2)->Clone("PartSigma1"); |
| 909 | PartMean1->SetTitle(";p /(GeV/c) ; Mean (Gauss)"); |
| 910 | PartSigma1->SetTitle(";p /(GeV/c) ; Sigma (Gauss)"); |
| 911 | |
| 912 | histPart1->SetMarkerStyle(22); |
| 913 | PartMean1->SetMarkerStyle(21); |
| 914 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 915 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 916 | |
| 917 | if(kParticle2==0){ // electron |
| 918 | nSigmaMinus2 = -1.999; |
| 919 | nSigmaPlus2 = 2.999; |
| 920 | } else if(kParticle2 == 1){ //pion |
| 921 | nSigmaMinus2 = -1.999; |
| 922 | nSigmaPlus2 = 2.999; |
| 923 | } else if(kParticle2 == 2){ //kaons |
| 924 | nSigmaMinus2 = -2.999; |
| 925 | nSigmaPlus2 = 2.999; |
| 926 | } else if(kParticle2 == 3){ //protons |
| 927 | nSigmaMinus2 = -2.999; |
| 928 | nSigmaPlus2 = 2.999; |
| 929 | } |
| 930 | |
| 931 | // |
| 932 | // 2. select and fit second particle |
| 933 | // |
| 934 | hist->GetAxis(3)->SetRangeUser(kParticle2,kParticle2); |
| 935 | hist->GetAxis(4)->SetRangeUser(nSigmaMinus2,nSigmaPlus2); |
| 936 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 937 | TH2D * histPart2 = (TH2D*)hist->Projection(1,0)->Clone("histPart2"); |
| 938 | histPart2->RebinX(4); |
| 939 | histPart2->FitSlicesY(0,0,-1,10,"QNR",&arr); |
| 940 | TH1D * PartMean2 = (TH1D *) arr.At(1)->Clone("PartMean2"); |
| 941 | TH1D * PartSigma2 = (TH1D *) arr.At(2)->Clone("PartSigma2"); |
| 942 | PartMean2->SetTitle(";p /(GeV/c) ; Mean (Gauss)"); |
| 943 | PartSigma2->SetTitle(";p /(GeV/c) ; Sigma (Gauss)"); |
| 944 | PartMean2->SetMarkerStyle(20); |
| 945 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 946 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 947 | |
| 948 | // |
| 949 | // 3. separation |
| 950 | // |
| 951 | |
| 952 | TH1F *fHistSeparation=(TH1F*)PartMean1->Clone("fHistSeparation"); //to get same binning |
| 953 | fHistSeparation->SetMarkerStyle(22); |
| 954 | const Int_t Nbins = PartMean1->GetNbinsX(); |
| 955 | fHistSeparation->SetTitle(";p /(GeV/c) ; Separation"); |
| 956 | |
| 957 | Float_t DeltaMean[Nbins] ; |
| 958 | Float_t DeltaSigma[Nbins]; |
| 959 | |
| 960 | for(Int_t i=0 ; i<Nbins; i++){ |
| 961 | |
| 962 | DeltaMean[i] = TMath::Abs(PartMean1->GetBinContent(i) - PartMean2->GetBinContent(i)); |
| 963 | DeltaSigma[i] = TMath::Abs((PartSigma1->GetBinContent(i) + PartSigma2->GetBinContent(i)))/2.; |
| 964 | |
| 965 | if(!(TMath::Abs(DeltaSigma[i])<0.000001))fHistSeparation->SetBinContent(i,DeltaMean[i]/DeltaSigma[i]); |
| 966 | }//for(Int_t i=0 ; i<Nbins ; i++) |
| 967 | |
| 968 | TObjArray *array = new TObjArray(); |
| 969 | array->Add(histPart1); |
| 970 | array->Add(histPart2); |
| 971 | array->Add(PartMean1); |
| 972 | array->Add(PartSigma1); |
| 973 | array->Add(PartMean2); |
| 974 | array->Add(PartSigma2); |
| 975 | array->Add(fHistSeparation); |
| 976 | return array; |
| 977 | |
| 978 | } |
| 979 | |
| 980 | |
| 981 | //________________________________________________________________________ |
| 982 | TObjArray * AliTPCcalibResidualPID::GetResidualGraphs(THnSparseF * histPidQA, const Char_t * system, Bool_t useV0s) { |
| 983 | // |
| 984 | // Extracts the residual graphs from THnSparse created from data (NOT MC!) |
| 985 | // |
| 986 | |
| 987 | const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})"; |
| 988 | const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)"; |
| 989 | |
| 990 | Bool_t isPPb = kFALSE; |
| 991 | if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { |
| 992 | Printf("p-Pb/Pb-p detected - Applying special handling!"); |
| 993 | isPPb = kTRUE; |
| 994 | } |
| 995 | |
| 996 | // the following parameters have to be extracted from the data itself |
| 997 | // |
| 998 | |
| 999 | Float_t nSigmaPionMinus = -3.999; |
| 1000 | Float_t nSigmaPionPlus = 3.999; |
| 1001 | Float_t nSigmaKaonMinus = -4.199; |
| 1002 | Float_t nSigmaKaonPlus = 4.7999; |
| 1003 | Float_t nSigmaElectronMinus = -1.999; |
| 1004 | Float_t nSigmaElectronPlus = 4.999; |
| 1005 | |
| 1006 | Int_t cutForFitting = 10; |
| 1007 | Double_t heightFractionForFittingRange = 0.1; |
| 1008 | // |
| 1009 | THnSparse * hist = histPidQA; |
| 1010 | // |
| 1011 | TCanvas * canvasQAtpc = new TCanvas("canvasQAtpcResGraph","Control canvas for residual graphs (TPC)",100,10,1380,800); |
| 1012 | canvasQAtpc->Divide(2,2); |
| 1013 | |
| 1014 | TCanvas * canvasQAtof = new TCanvas("canvasQAtofResGraph","Control canvas for residual graphs (TOF)",100,10,1380,800); |
| 1015 | canvasQAtof->Divide(2,2); |
| 1016 | |
| 1017 | TCanvas * canvasQAv0 = new TCanvas("canvasQAv0ResGraph","Control canvas for residual graphs (V0)",100,10,1380,800); |
| 1018 | canvasQAv0->Divide(2,2); |
| 1019 | |
| 1020 | TCanvas * canvasQAv0plusTOF = new TCanvas("canvasQAv0plusTOFResGraph","Control canvas for residual graphs (V0+TOF)",100,10,1380,800); |
| 1021 | canvasQAv0plusTOF->Divide(2,2); |
| 1022 | |
| 1023 | |
| 1024 | TCanvas * canvasQAv0DeDxPurityEl = new TCanvas("canvasQAv0DeDxPurityEl","Control canvas for residual graphs (V0 el)",100,10,600,400); |
| 1025 | TCanvas * canvasQAv0DeDxPurityPi = new TCanvas("canvasQAv0DeDxPurityPi","Control canvas for residual graphs (V0 pi)",100,10,600,400); |
| 1026 | TCanvas * canvasQAv0DeDxPurityPr = new TCanvas("canvasQAv0DeDxPurityPr","Control canvas for residual graphs (V0 pr)",100,10,600,400); |
| 1027 | |
| 1028 | canvasQAv0DeDxPurityEl->SetLogx(); |
| 1029 | canvasQAv0DeDxPurityPi->SetLogx(); |
| 1030 | canvasQAv0DeDxPurityPr->SetLogx(); |
| 1031 | |
| 1032 | canvasQAv0DeDxPurityEl->SetLogz(); |
| 1033 | canvasQAv0DeDxPurityPi->SetLogz(); |
| 1034 | canvasQAv0DeDxPurityPr->SetLogz(); |
| 1035 | |
| 1036 | canvasQAv0DeDxPurityEl->SetTopMargin(0.03); |
| 1037 | canvasQAv0DeDxPurityPi->SetTopMargin(0.03); |
| 1038 | canvasQAv0DeDxPurityPr->SetTopMargin(0.03); |
| 1039 | |
| 1040 | for (Int_t i = 1; i <= 4; i++) { |
| 1041 | canvasQAtpc->GetPad(i)->SetGrid(1, 1); |
| 1042 | canvasQAtof->GetPad(i)->SetGrid(1, 1); |
| 1043 | canvasQAv0->GetPad(i)->SetGrid(1, 1); |
| 1044 | canvasQAv0plusTOF->GetPad(i)->SetGrid(1, 1); |
| 1045 | |
| 1046 | canvasQAtpc->GetPad(i)->SetLogz(); |
| 1047 | canvasQAtof->GetPad(i)->SetLogz(); |
| 1048 | canvasQAv0->GetPad(i)->SetLogz(); |
| 1049 | canvasQAv0plusTOF->GetPad(i)->SetLogz(); |
| 1050 | |
| 1051 | canvasQAtpc->GetPad(i)->SetLogx(); |
| 1052 | canvasQAtof->GetPad(i)->SetLogx(); |
| 1053 | canvasQAv0->GetPad(i)->SetLogx(); |
| 1054 | canvasQAv0plusTOF->GetPad(i)->SetLogx(); |
| 1055 | } |
| 1056 | // |
| 1057 | // 1. select and fit electrons |
| 1058 | // |
| 1059 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s |
| 1060 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons |
| 1061 | hist->GetAxis(4)->SetRangeUser(nSigmaElectronMinus,nSigmaElectronPlus); // 3sigma in TPC |
| 1062 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1063 | TH2D * histElectron = hist->Projection(1,0); |
| 1064 | histElectron->SetName("histElectron"); |
| 1065 | histElectron->RebinX(4); |
| 1066 | histElectron->GetXaxis()->SetRangeUser(0.2,3.0); |
| 1067 | TObjArray arr; |
| 1068 | FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1069 | TH1D * electronPoints = (TH1D *) arr.At(1)->Clone(); |
| 1070 | // |
| 1071 | TGraphErrors * electronGraph = new TGraphErrors(electronPoints); |
| 1072 | electronGraph->SetName("electronGraph"); |
| 1073 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { |
| 1074 | Bool_t removePoint = electronGraph->GetY()[ip] < 10 || |
| 1075 | electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.05 || |
| 1076 | electronGraph->GetX()[ip] > 2.0 || |
| 1077 | electronGraph->GetX()[ip] < 0.5; |
| 1078 | if (removePoint) { |
| 1079 | electronGraph->RemovePoint(ip); |
| 1080 | ip--; |
| 1081 | continue; |
| 1082 | } |
| 1083 | } |
| 1084 | // |
| 1085 | canvasQAtof->cd(1); |
| 1086 | histElectron->SetTitle("electrons"); |
| 1087 | histElectron->GetYaxis()->SetRangeUser(50, 120); |
| 1088 | histElectron->GetXaxis()->SetTitle(momTitle.Data()); |
| 1089 | histElectron->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1090 | histElectron->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1091 | histElectron->GetXaxis()->SetNoExponent(kTRUE); |
| 1092 | histElectron->Draw("colz"); |
| 1093 | electronPoints->SetMarkerStyle(24); |
| 1094 | electronPoints->Draw("same"); |
| 1095 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 1096 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 1097 | // |
| 1098 | // 2. protons |
| 1099 | // |
| 1100 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s |
| 1101 | hist->GetAxis(3)->SetRangeUser(3,3); // protons |
| 1102 | // |
| 1103 | //hist->GetAxis(4)->SetRangeUser(nSigmaProtonMinus,nSigmaProtonPlus); // protons --> not reliable, use PATTERN RECOGNITION |
| 1104 | TH2D * histProtonTPC = hist->Projection(1,0); |
| 1105 | histProtonTPC->SetName("histProtonTPC"); |
| 1106 | histProtonTPC->GetXaxis()->SetRangeUser(0.15,0.7); |
| 1107 | histProtonTPC->GetYaxis()->SetRangeUser(50, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins())); |
| 1108 | |
| 1109 | // PATTERN RECOGNITION |
| 1110 | //OLD TF1 betaSq("betaSq","50./TMath::Power(x,1.3)",0.1,10); |
| 1111 | // Add some additional Erf - cuts away kaons for pPb and PbPb and does not hurt in pp |
| 1112 | TF1 betaSq("betaSq","50./TMath::Power(x,1.3) + (1-TMath::Erf((x-0.2) / 0.075)) * 100",0.1,10); |
| 1113 | for(Int_t ix = 1; ix <= histProtonTPC->GetXaxis()->GetNbins(); ix++) { |
| 1114 | for(Int_t jy = 1; jy <= histProtonTPC->GetYaxis()->GetNbins(); jy++) { |
| 1115 | Float_t yPos = histProtonTPC->GetYaxis()->GetBinCenter(jy); |
| 1116 | Float_t xPos = histProtonTPC->GetXaxis()->GetBinCenter(ix); |
| 1117 | Int_t bin = histProtonTPC->GetBin(ix,jy); |
| 1118 | if (yPos < betaSq.Eval(xPos)) histProtonTPC->SetBinContent(bin,0); |
| 1119 | } |
| 1120 | } |
| 1121 | // |
| 1122 | //TODO Use likelihood option -> Inside fitting range ~no contamination, but little statistics at low momenta requires L |
| 1123 | FitSlicesY(histProtonTPC, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); |
| 1124 | TH1D * protonPointsTPC = (TH1D *) arr.At(1)->Clone(); |
| 1125 | // |
| 1126 | TGraphErrors * protonTpcGraph = new TGraphErrors(protonPointsTPC); |
| 1127 | protonTpcGraph->SetName("protonTpcGraph"); |
| 1128 | for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) { |
| 1129 | Bool_t removePoint = protonTpcGraph->GetY()[ip] < 10 || protonTpcGraph->GetEY()[ip]/protonTpcGraph->GetY()[ip] > 0.05 // TODO Larger tolerance, since this is only for the correction function, where 5% are still better than no data point |
| 1130 | || protonTpcGraph->GetX()[ip] > (useV0s ? (isPPb ? 0.499 : 0.349) : 0.649); // If V0s are to be used, don't use TPC only protons to too high momenta; for pPb low statistics for the moment, therefore, go a bit further with TPC protons |
| 1131 | //TODO NOW -> Changed range of TPC-signal -> Hopefully, sufficient statistics now for very low momenta!|| protonTpcGraph->GetX()[ip] < 0.19; // Added this cut because almost always bad statistics below 0.19 |
| 1132 | if (removePoint) { |
| 1133 | protonTpcGraph->RemovePoint(ip); |
| 1134 | ip--; |
| 1135 | continue; |
| 1136 | } |
| 1137 | } |
| 1138 | // |
| 1139 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1140 | TH2D * histProtonTOF = hist->Projection(1,0); |
| 1141 | histProtonTOF->SetName("histProtonTOF"); |
| 1142 | histProtonTOF->GetXaxis()->SetRangeUser(0.6,2.5); |
| 1143 | |
| 1144 | // Same pattern recognition as before |
| 1145 | for(Int_t ix = 1; ix <= histProtonTOF->GetXaxis()->GetNbins(); ix++) { |
| 1146 | for(Int_t jy = 1; jy <= histProtonTOF->GetYaxis()->GetNbins(); jy++) { |
| 1147 | Float_t yPos = histProtonTOF->GetYaxis()->GetBinCenter(jy); |
| 1148 | Float_t xPos = histProtonTOF->GetXaxis()->GetBinCenter(ix); |
| 1149 | Int_t bin = histProtonTOF->GetBin(ix,jy); |
| 1150 | if (yPos < betaSq.Eval(xPos)) histProtonTOF->SetBinContent(bin,0); |
| 1151 | } |
| 1152 | } |
| 1153 | |
| 1154 | FitSlicesY(histProtonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1155 | |
| 1156 | TH1D * protonPointsTOF = (TH1D *)arr.At(1)->Clone(); |
| 1157 | // |
| 1158 | TGraphErrors * protonTofGraph = new TGraphErrors(protonPointsTOF); |
| 1159 | protonTofGraph->SetName("protonTofGraph"); |
| 1160 | for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) { |
| 1161 | // If V0s are to be used, TPC+TOF protons are only for reference and can be plotted to higher momenta |
| 1162 | Bool_t removePoint = protonTofGraph->GetY()[ip] < 10 || protonTofGraph->GetEY()[ip]/protonTofGraph->GetY()[ip] > 0.02 || protonTofGraph->GetX()[ip] > (useV0s ? 3 : 2) || protonTofGraph->GetX()[ip] < 0.65; |
| 1163 | if (removePoint) { |
| 1164 | protonTofGraph->RemovePoint(ip); |
| 1165 | ip--; |
| 1166 | continue; |
| 1167 | } |
| 1168 | } |
| 1169 | // |
| 1170 | canvasQAtpc->cd(2); |
| 1171 | histProtonTPC->SetTitle("protons"); |
| 1172 | histProtonTPC->GetXaxis()->SetTitle(momTitle.Data()); |
| 1173 | histProtonTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1174 | histProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1175 | histProtonTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 1176 | histProtonTPC->Draw("colz"); |
| 1177 | betaSq.DrawCopy("same"); |
| 1178 | protonPointsTPC->SetMarkerStyle(20); |
| 1179 | protonPointsTOF->SetMarkerStyle(24); |
| 1180 | protonPointsTOF->Draw("same"); |
| 1181 | protonPointsTPC->Draw("same"); |
| 1182 | // |
| 1183 | protonTpcGraph->SetMarkerStyle(26); |
| 1184 | protonTpcGraph->SetMarkerColor(kMagenta); |
| 1185 | protonTpcGraph->DrawClone("p"); |
| 1186 | protonTofGraph->SetMarkerStyle(25); |
| 1187 | protonTofGraph->SetMarkerColor(kMagenta); |
| 1188 | protonTofGraph->DrawClone("p"); |
| 1189 | |
| 1190 | canvasQAtof->cd(2); |
| 1191 | histProtonTOF->SetTitle("protons"); |
| 1192 | histProtonTOF->GetYaxis()->SetRangeUser(30, 250); |
| 1193 | histProtonTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1194 | histProtonTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1195 | histProtonTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1196 | histProtonTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1197 | histProtonTOF->Draw("colz"); |
| 1198 | betaSq.DrawCopy("same"); |
| 1199 | protonPointsTOF->Draw("same"); |
| 1200 | protonPointsTPC->Draw("same"); |
| 1201 | // |
| 1202 | protonTpcGraph->DrawClone("p"); |
| 1203 | protonTofGraph->DrawClone("p"); |
| 1204 | |
| 1205 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 1206 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 1207 | // |
| 1208 | // 3. pions |
| 1209 | // |
| 1210 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s |
| 1211 | hist->GetAxis(3)->SetRangeUser(1,1); // pions |
| 1212 | // |
| 1213 | Double_t lowerPionThreshold = 0.3; |
| 1214 | |
| 1215 | hist->GetAxis(4)->SetRangeUser(nSigmaPionMinus,nSigmaPionPlus); // pions |
| 1216 | TH2D * histPionTPC = hist->Projection(1,0); |
| 1217 | histPionTPC->SetName("histPionTPC"); |
| 1218 | histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7); |
| 1219 | FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1220 | |
| 1221 | // In case of really bad splines comment last 5 lines and use the following 6 lines: |
| 1222 | //TH2D * histPionTPC = hist->Projection(1,0); |
| 1223 | //histPionTPC->SetName("histPionTPC"); |
| 1224 | //histPionTPC->GetYaxis()->SetRangeUser(30,100); |
| 1225 | //histPionTPC->GetXaxis()->SetRangeUser(0.15,0.7); |
| 1226 | //FitSlicesY(histPionTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1227 | //lowerPionThreshold = 0.15; |
| 1228 | |
| 1229 | |
| 1230 | TH1D * pionPointsTPC = (TH1D *) arr.At(1)->Clone(); |
| 1231 | // |
| 1232 | TGraphErrors * pionTpcGraph = new TGraphErrors(pionPointsTPC); |
| 1233 | pionTpcGraph->SetName("pionTpcGraph"); |
| 1234 | for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) { |
| 1235 | Bool_t removePoint = pionTpcGraph->GetY()[ip] < 10 || pionTpcGraph->GetEY()[ip]/pionTpcGraph->GetY()[ip] > 0.02 || pionTpcGraph->GetX()[ip] > 0.5 |
| 1236 | || pionTpcGraph->GetX()[ip] < lowerPionThreshold; // Exclude contamination by electrons |
| 1237 | if (removePoint) { |
| 1238 | pionTpcGraph->RemovePoint(ip); |
| 1239 | ip--; |
| 1240 | continue; |
| 1241 | } |
| 1242 | } |
| 1243 | // |
| 1244 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1245 | TH2D * histPionTOF = hist->Projection(1,0); |
| 1246 | histPionTOF->SetName("histPionTOF"); |
| 1247 | histPionTOF->GetXaxis()->SetRangeUser(0.5,1.1); |
| 1248 | FitSlicesY(histPionTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1249 | TH1D * pionPointsTOF = (TH1D *) arr.At(1)->Clone(); |
| 1250 | TGraphErrors * pionTofGraph = new TGraphErrors(pionPointsTOF); |
| 1251 | pionTofGraph->SetName("pionTofGraph"); |
| 1252 | for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) { |
| 1253 | Bool_t removePoint = pionTofGraph->GetY()[ip] < 10 || pionTofGraph->GetEY()[ip]/pionTofGraph->GetY()[ip] > 0.02 || pionTofGraph->GetX()[ip] > 1.1 || pionTofGraph->GetX()[ip] < 0.5; // Changed by Ben (was 0.35) to avoid problems with TOF efficiency/systematics |
| 1254 | if (removePoint) { |
| 1255 | pionTofGraph->RemovePoint(ip); |
| 1256 | ip--; |
| 1257 | continue; |
| 1258 | } |
| 1259 | } |
| 1260 | canvasQAtpc->cd(3); |
| 1261 | histPionTPC->GetYaxis()->SetRangeUser(30, 90); |
| 1262 | histPionTPC->SetTitle("pions"); |
| 1263 | histPionTPC->GetXaxis()->SetTitle(momTitle.Data()); |
| 1264 | histPionTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1265 | histPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1266 | histPionTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 1267 | histPionTPC->Draw("colz"); |
| 1268 | pionPointsTPC->SetMarkerStyle(20); |
| 1269 | pionPointsTOF->SetMarkerStyle(24); |
| 1270 | pionPointsTOF->Draw("same"); |
| 1271 | pionPointsTPC->Draw("same"); |
| 1272 | // |
| 1273 | pionTpcGraph->SetMarkerStyle(26); |
| 1274 | pionTpcGraph->SetMarkerColor(kMagenta); |
| 1275 | pionTpcGraph->DrawClone("p"); |
| 1276 | pionTofGraph->SetMarkerStyle(25); |
| 1277 | pionTofGraph->SetMarkerColor(kMagenta); |
| 1278 | pionTofGraph->DrawClone("p"); |
| 1279 | |
| 1280 | canvasQAtof->cd(3); |
| 1281 | histPionTOF->GetYaxis()->SetRangeUser(30, 80); |
| 1282 | histPionTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1283 | histPionTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1284 | histPionTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1285 | histPionTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1286 | histPionTOF->Draw("colz"); |
| 1287 | histPionTOF->SetTitle("pions"); |
| 1288 | pionPointsTOF->Draw("same"); |
| 1289 | pionPointsTPC->Draw("same"); |
| 1290 | // |
| 1291 | pionTpcGraph->DrawClone("p"); |
| 1292 | pionTofGraph->DrawClone("p"); |
| 1293 | |
| 1294 | |
| 1295 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 1296 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 1297 | // |
| 1298 | // 4. kaons |
| 1299 | // |
| 1300 | hist->GetAxis(2)->SetRangeUser(0,0); // Non-V0s |
| 1301 | hist->GetAxis(3)->SetRangeUser(2,2); // kaons |
| 1302 | // |
| 1303 | hist->GetAxis(4)->SetRangeUser(nSigmaKaonMinus,nSigmaKaonPlus); // kaons |
| 1304 | TH2D * histKaonTPC = hist->Projection(1,0); |
| 1305 | histKaonTPC->SetName("histKaonTPC"); |
| 1306 | histKaonTPC->GetXaxis()->SetRangeUser(0.12,0.6999); |
| 1307 | FitSlicesY(histKaonTPC, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1308 | TH1D * kaonPointsTPC = (TH1D*) arr.At(1)->Clone(); |
| 1309 | // |
| 1310 | TGraphErrors * kaonTpcGraph = new TGraphErrors(kaonPointsTPC); |
| 1311 | kaonTpcGraph->SetName("kaonTpcGraph"); |
| 1312 | for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) { |
| 1313 | Bool_t removePoint = kaonTpcGraph->GetY()[ip] < 10 || kaonTpcGraph->GetEY()[ip]/kaonTpcGraph->GetY()[ip] > 0.02 |
| 1314 | || kaonTpcGraph->GetX()[ip] < 0.12 || kaonTpcGraph->GetX()[ip] > 0.319; //TODO BEN was > 0.399 |
| 1315 | if (removePoint) { |
| 1316 | kaonTpcGraph->RemovePoint(ip); |
| 1317 | ip--; |
| 1318 | continue; |
| 1319 | } |
| 1320 | } |
| 1321 | // |
| 1322 | |
| 1323 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // sigma in TOF |
| 1324 | |
| 1325 | /* |
| 1326 | hist->GetAxis(5)->SetRangeUser(-1.999,1.999); // sigma in TOF |
| 1327 | |
| 1328 | if (hist->GetAxis(3)->GetNbins() >= 5) { |
| 1329 | printf("Using restricted eta range for TPC+TOF Kaons!\n"); |
| 1330 | hist->GetAxis(3)->SetRangeUser(4,4); |
| 1331 | } |
| 1332 | else { |
| 1333 | printf("WARNING: Data points for restricted eta range for TPC+TOF Kaons not available! Using full eta range (worse w.r.t. contamination)\n"); |
| 1334 | }*/ |
| 1335 | |
| 1336 | TH2D * histKaonTOF = hist->Projection(1,0); |
| 1337 | histKaonTOF->SetName("histKaonTOF"); |
| 1338 | histKaonTOF->GetXaxis()->SetRangeUser(0.3,1.5); |
| 1339 | FitSlicesY(histKaonTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1340 | TH1D * kaonPointsTOF = (TH1D*) arr.At(1)->Clone(); |
| 1341 | // |
| 1342 | TGraphErrors * kaonTofGraph = new TGraphErrors(kaonPointsTOF); |
| 1343 | kaonTofGraph->SetName("kaonTofGraph"); |
| 1344 | // |
| 1345 | for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) { |
| 1346 | Bool_t removePoint = kaonTofGraph->GetY()[ip] < 10 || kaonTofGraph->GetEY()[ip]/kaonTofGraph->GetY()[ip] > 0.02 || kaonTofGraph->GetX()[ip] > 1.0 |
| 1347 | || kaonTofGraph->GetX()[ip] < 0.36; //TODO BEN NOW was 0.5 |
| 1348 | if (removePoint) { |
| 1349 | kaonTofGraph->RemovePoint(ip); |
| 1350 | ip--; |
| 1351 | continue; |
| 1352 | } |
| 1353 | } |
| 1354 | // |
| 1355 | canvasQAtpc->cd(4); |
| 1356 | histKaonTPC->SetTitle("kaons"); |
| 1357 | histKaonTPC->GetYaxis()->SetRangeUser(30, 500); |
| 1358 | histKaonTPC->GetXaxis()->SetTitle(momTitle.Data()); |
| 1359 | histKaonTPC->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1360 | histKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1361 | histKaonTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 1362 | histKaonTPC->Draw("colz"); |
| 1363 | kaonPointsTPC->SetMarkerStyle(20); |
| 1364 | kaonPointsTOF->SetMarkerStyle(24); |
| 1365 | kaonPointsTOF->Draw("same"); |
| 1366 | kaonPointsTPC->Draw("same"); |
| 1367 | // |
| 1368 | kaonTpcGraph->SetMarkerStyle(26); |
| 1369 | kaonTpcGraph->SetMarkerColor(kMagenta); |
| 1370 | kaonTpcGraph->DrawClone("p"); |
| 1371 | kaonTofGraph->SetMarkerStyle(25); |
| 1372 | kaonTofGraph->SetMarkerColor(kMagenta); |
| 1373 | kaonTofGraph->DrawClone("p"); |
| 1374 | |
| 1375 | |
| 1376 | canvasQAtof->cd(4); |
| 1377 | histKaonTOF->GetYaxis()->SetRangeUser(30, 250); |
| 1378 | histKaonTOF->SetTitle("kaons"); |
| 1379 | histKaonTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1380 | histKaonTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1381 | histKaonTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1382 | histKaonTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1383 | histKaonTOF->Draw("colz"); |
| 1384 | kaonPointsTOF->Draw("same"); |
| 1385 | kaonPointsTPC->Draw("same"); |
| 1386 | kaonTpcGraph->DrawClone("p"); |
| 1387 | kaonTofGraph->DrawClone("p"); |
| 1388 | |
| 1389 | hist->GetAxis(4)->SetRange(0,-1); // RESET RANGES |
| 1390 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 1391 | |
| 1392 | |
| 1393 | // |
| 1394 | // 5. V0 electrons |
| 1395 | // |
| 1396 | hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons |
| 1397 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons |
| 1398 | |
| 1399 | // |
| 1400 | TH2D * histElectronV0 = hist->Projection(1,0); |
| 1401 | histElectronV0->SetName("histElectronV0"); |
| 1402 | //histElectronV0->RebinX(2); |
| 1403 | histElectronV0->GetXaxis()->SetRangeUser(0.15,5.0); |
| 1404 | FitSlicesY(histElectronV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1405 | TH1D * electronPointsV0 = (TH1D*) arr.At(1)->Clone(); |
| 1406 | // |
| 1407 | TGraphErrors * electronV0Graph = new TGraphErrors(electronPointsV0); |
| 1408 | electronV0Graph->SetName("electronV0Graph"); |
| 1409 | for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) { |
| 1410 | Bool_t removePoint = electronV0Graph->GetY()[ip] < 10 || electronV0Graph->GetEY()[ip]/electronV0Graph->GetY()[ip] > 0.02; |
| 1411 | |
| 1412 | if (removePoint) { |
| 1413 | electronV0Graph->RemovePoint(ip); |
| 1414 | ip--; |
| 1415 | continue; |
| 1416 | } |
| 1417 | } |
| 1418 | // |
| 1419 | canvasQAv0->cd(1); |
| 1420 | histElectronV0->SetTitle("V0 electrons"); |
| 1421 | histElectronV0->GetYaxis()->SetRangeUser(50, 120); |
| 1422 | histElectronV0->GetXaxis()->SetTitle(momTitle.Data()); |
| 1423 | histElectronV0->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1424 | histElectronV0->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1425 | histElectronV0->GetXaxis()->SetNoExponent(kTRUE); |
| 1426 | histElectronV0->SetStats(0); |
| 1427 | histElectronV0->Draw("colz"); |
| 1428 | electronPointsV0->SetMarkerStyle(34); |
| 1429 | electronPointsV0->Draw("same"); |
| 1430 | electronGraph->SetMarkerStyle(25); |
| 1431 | electronGraph->SetMarkerColor(kMagenta); |
| 1432 | electronGraph->DrawClone("p"); |
| 1433 | |
| 1434 | canvasQAv0DeDxPurityEl->cd(); |
| 1435 | gStyle->SetOptTitle(0); |
| 1436 | TH1D* hElNew = (TH1D*)histElectronV0->DrawClone("colz"); |
| 1437 | hElNew->GetXaxis()->SetTitleOffset(1.0); |
| 1438 | hElNew->SetTitle(""); |
| 1439 | electronPointsV0->DrawClone("same"); |
| 1440 | gStyle->SetOptTitle(1); |
| 1441 | |
| 1442 | |
| 1443 | |
| 1444 | canvasQAtof->cd(1); |
| 1445 | electronV0Graph->SetMarkerStyle(24); |
| 1446 | electronV0Graph->SetMarkerColor(kMagenta); |
| 1447 | electronV0Graph->DrawClone("p"); |
| 1448 | |
| 1449 | // |
| 1450 | // 6. V0 pions |
| 1451 | // |
| 1452 | hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions |
| 1453 | hist->GetAxis(3)->SetRangeUser(1,1); // pions |
| 1454 | // |
| 1455 | TH2D * histPionV0 = hist->Projection(1,0); |
| 1456 | histPionV0->SetName("histPionV0"); |
| 1457 | histPionV0->GetXaxis()->SetRangeUser(0.15,10.); |
| 1458 | FitSlicesY(histPionV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1459 | TH1D * pionPointsV0 = (TH1D*) arr.At(1)->Clone(); |
| 1460 | // |
| 1461 | TGraphErrors * pionV0Graph = new TGraphErrors(pionPointsV0); |
| 1462 | pionV0Graph->SetName("pionV0Graph"); |
| 1463 | for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) { |
| 1464 | Bool_t removePoint = pionV0Graph->GetY()[ip] < 10 || pionV0Graph->GetEY()[ip]/pionV0Graph->GetY()[ip] > 0.02; |
| 1465 | |
| 1466 | if (removePoint) { |
| 1467 | pionV0Graph->RemovePoint(ip); |
| 1468 | ip--; |
| 1469 | continue; |
| 1470 | } |
| 1471 | } |
| 1472 | // |
| 1473 | canvasQAv0->cd(3); |
| 1474 | histPionV0->SetTitle("V0 pions"); |
| 1475 | histPionV0->GetYaxis()->SetRangeUser(30, 100); |
| 1476 | histPionV0->GetXaxis()->SetTitle(momTitle.Data()); |
| 1477 | histPionV0->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1478 | histPionV0->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1479 | histPionV0->GetXaxis()->SetNoExponent(kTRUE); |
| 1480 | histPionV0->Draw("colz"); |
| 1481 | pionPointsV0->SetMarkerStyle(34); |
| 1482 | pionPointsV0->Draw("same"); |
| 1483 | pionTofGraph->DrawClone("p"); |
| 1484 | pionTpcGraph->DrawClone("p"); |
| 1485 | |
| 1486 | canvasQAv0DeDxPurityPi->cd(); |
| 1487 | gStyle->SetOptTitle(0); |
| 1488 | TH1D* hPiNew = (TH1D*)histPionV0->DrawClone("colz"); |
| 1489 | hPiNew->GetXaxis()->SetTitleOffset(1.0); |
| 1490 | hPiNew->SetTitle(""); |
| 1491 | pionPointsV0->DrawClone("same"); |
| 1492 | gStyle->SetOptTitle(1); |
| 1493 | |
| 1494 | canvasQAtof->cd(3); |
| 1495 | pionV0Graph->SetMarkerStyle(24); |
| 1496 | pionV0Graph->SetMarkerColor(kMagenta); |
| 1497 | pionV0Graph->DrawClone("p"); |
| 1498 | |
| 1499 | canvasQAtpc->cd(3); |
| 1500 | pionV0Graph->DrawClone("p"); |
| 1501 | |
| 1502 | // |
| 1503 | // 6. V0 protons |
| 1504 | // |
| 1505 | hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons |
| 1506 | hist->GetAxis(3)->SetRangeUser(3,3); // protons |
| 1507 | // |
| 1508 | TH2D * histProtonV0 = hist->Projection(1,0); |
| 1509 | histProtonV0->SetName("histProtonV0"); |
| 1510 | histProtonV0->GetXaxis()->SetRangeUser(0.15,10.); |
| 1511 | |
| 1512 | // Remove contamination due to el and pions and low momenta because there the statistics |
| 1513 | // for the V0 protons goes down and becomes comparable with the contamination |
| 1514 | // -> This spoils the fit, if the contamination is not removed |
| 1515 | |
| 1516 | // PATTERN RECOGNITION |
| 1517 | Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.6); |
| 1518 | |
| 1519 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { |
| 1520 | for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) { |
| 1521 | Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy); |
| 1522 | Float_t xPos = histProtonV0->GetXaxis()->GetBinCenter(ix); |
| 1523 | Int_t bin = histProtonV0->GetBin(ix,jy); |
| 1524 | if (yPos < betaSq.Eval(xPos)) histProtonV0->SetBinContent(bin,0); |
| 1525 | } |
| 1526 | } |
| 1527 | |
| 1528 | /* |
| 1529 | Int_t correctionUpToBin = histProtonV0->GetXaxis()->FindBin(0.4); |
| 1530 | Double_t threshold = 120.; |
| 1531 | |
| 1532 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { |
| 1533 | for(Int_t jy = 1; jy <= histProtonV0->GetYaxis()->GetNbins(); jy++) { |
| 1534 | Float_t yPos = histProtonV0->GetYaxis()->GetBinCenter(jy); |
| 1535 | Int_t bin = histProtonV0->GetBin(ix,jy); |
| 1536 | if (yPos < threshold) histProtonV0->SetBinContent(bin,0); |
| 1537 | else break; |
| 1538 | } |
| 1539 | } |
| 1540 | */ |
| 1541 | |
| 1542 | FitSlicesY(histProtonV0, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1543 | TH1D * protonPointsV0 = (TH1D*) arr.At(1)->Clone(); |
| 1544 | // |
| 1545 | TGraphErrors * protonV0Graph = new TGraphErrors(protonPointsV0); |
| 1546 | protonV0Graph->SetName("protonV0Graph"); |
| 1547 | for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) { |
| 1548 | Bool_t removePoint = protonV0Graph->GetY()[ip] < 10 || protonV0Graph->GetEY()[ip]/protonV0Graph->GetY()[ip] > 0.02; |
| 1549 | |
| 1550 | if (removePoint) { |
| 1551 | protonV0Graph->RemovePoint(ip); |
| 1552 | ip--; |
| 1553 | continue; |
| 1554 | } |
| 1555 | } |
| 1556 | // |
| 1557 | canvasQAv0->cd(2); |
| 1558 | histProtonV0->SetTitle("V0 protons"); |
| 1559 | histProtonV0->GetXaxis()->SetTitle(momTitle.Data()); |
| 1560 | histProtonV0->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1561 | histProtonV0->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1562 | histProtonV0->GetXaxis()->SetNoExponent(kTRUE); |
| 1563 | histProtonV0->Draw("colz"); |
| 1564 | protonPointsV0->SetMarkerStyle(34); |
| 1565 | protonPointsV0->Draw("same"); |
| 1566 | protonTofGraph->DrawClone("p"); |
| 1567 | protonTpcGraph->DrawClone("p"); |
| 1568 | |
| 1569 | canvasQAv0DeDxPurityPr->cd(); |
| 1570 | gStyle->SetOptTitle(0); |
| 1571 | TH1D* hPrNew = (TH1D*)histProtonV0->DrawClone("colz"); |
| 1572 | hPrNew->GetXaxis()->SetTitleOffset(1.0); |
| 1573 | hPrNew->SetTitle(""); |
| 1574 | protonPointsV0->DrawClone("same"); |
| 1575 | gStyle->SetOptTitle(1); |
| 1576 | |
| 1577 | canvasQAtof->cd(2); |
| 1578 | protonV0Graph->SetMarkerStyle(24); |
| 1579 | protonV0Graph->SetMarkerColor(kMagenta); |
| 1580 | protonV0Graph->DrawClone("p"); |
| 1581 | |
| 1582 | canvasQAtpc->cd(2); |
| 1583 | protonV0Graph->DrawClone("p"); |
| 1584 | |
| 1585 | |
| 1586 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES |
| 1587 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES |
| 1588 | |
| 1589 | |
| 1590 | // |
| 1591 | // 5. V0 electrons + TOF |
| 1592 | // |
| 1593 | hist->GetAxis(2)->SetRangeUser(1,1); // V0 electrons |
| 1594 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons |
| 1595 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1596 | // |
| 1597 | TH2D * histElectronV0plusTOF = hist->Projection(1,0); |
| 1598 | histElectronV0plusTOF->SetName("histElectronV0plusTOF"); |
| 1599 | histElectronV0plusTOF->RebinX(2); |
| 1600 | histElectronV0plusTOF->GetXaxis()->SetRangeUser(0.2,5.0); |
| 1601 | FitSlicesY(histElectronV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1602 | TH1D * electronPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); |
| 1603 | // |
| 1604 | TGraphErrors * electronV0plusTOFGraph = new TGraphErrors(electronPointsV0plusTOF); |
| 1605 | electronV0plusTOFGraph->SetName("electronV0plusTOFGraph"); |
| 1606 | for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) { |
| 1607 | Bool_t removePoint = electronV0plusTOFGraph->GetY()[ip] < 10 || electronV0plusTOFGraph->GetEY()[ip]/electronV0plusTOFGraph->GetY()[ip] > 0.02; |
| 1608 | |
| 1609 | if (removePoint) { |
| 1610 | electronV0plusTOFGraph->RemovePoint(ip); |
| 1611 | ip--; |
| 1612 | continue; |
| 1613 | } |
| 1614 | } |
| 1615 | // |
| 1616 | canvasQAv0plusTOF->cd(1); |
| 1617 | histElectronV0plusTOF->SetTitle("V0+TOF electrons"); |
| 1618 | histElectronV0plusTOF->GetYaxis()->SetRangeUser(50, 120); |
| 1619 | histElectronV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1620 | histElectronV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1621 | histElectronV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1622 | histElectronV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1623 | histElectronV0plusTOF->Draw("colz"); |
| 1624 | electronPointsV0plusTOF->SetMarkerStyle(29); |
| 1625 | electronPointsV0plusTOF->Draw("same"); |
| 1626 | electronGraph->DrawClone("p"); |
| 1627 | electronV0Graph->DrawClone("p"); |
| 1628 | |
| 1629 | canvasQAv0->cd(1); |
| 1630 | electronV0plusTOFGraph->SetMarkerStyle(30); |
| 1631 | electronV0plusTOFGraph->SetMarkerColor(kMagenta); |
| 1632 | electronV0plusTOFGraph->DrawClone("p"); |
| 1633 | |
| 1634 | canvasQAtof->cd(1); |
| 1635 | electronV0plusTOFGraph->DrawClone("p"); |
| 1636 | |
| 1637 | // |
| 1638 | // 6. V0 pions + TOF |
| 1639 | // |
| 1640 | hist->GetAxis(2)->SetRangeUser(2,2); // V0 pions |
| 1641 | hist->GetAxis(3)->SetRangeUser(1,1); // pions |
| 1642 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1643 | // |
| 1644 | TH2D * histPionV0plusTOF = hist->Projection(1,0); |
| 1645 | histPionV0plusTOF->SetName("histPionV0plusTOF"); |
| 1646 | histPionV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.); |
| 1647 | FitSlicesY(histPionV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1648 | TH1D * pionPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); |
| 1649 | // |
| 1650 | TGraphErrors * pionV0plusTOFGraph = new TGraphErrors(pionPointsV0plusTOF); |
| 1651 | pionV0plusTOFGraph->SetName("pionV0plusTOFGraph"); |
| 1652 | for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) { |
| 1653 | Bool_t removePoint = pionV0plusTOFGraph->GetY()[ip] < 10 || pionV0plusTOFGraph->GetEY()[ip]/pionV0plusTOFGraph->GetY()[ip] > 0.02; |
| 1654 | |
| 1655 | if (removePoint) { |
| 1656 | pionV0plusTOFGraph->RemovePoint(ip); |
| 1657 | ip--; |
| 1658 | continue; |
| 1659 | } |
| 1660 | } |
| 1661 | // |
| 1662 | canvasQAv0plusTOF->cd(3); |
| 1663 | histPionV0plusTOF->SetTitle("V0+TOF pions"); |
| 1664 | histPionV0plusTOF->GetYaxis()->SetRangeUser(30, 100); |
| 1665 | histPionV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1666 | histPionV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1667 | histPionV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1668 | histPionV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1669 | histPionV0plusTOF->Draw("colz"); |
| 1670 | pionPointsV0plusTOF->SetMarkerStyle(29); |
| 1671 | pionPointsV0plusTOF->Draw("same"); |
| 1672 | pionV0Graph->DrawClone("p"); |
| 1673 | pionTofGraph->DrawClone("p"); |
| 1674 | pionTpcGraph->DrawClone("p"); |
| 1675 | |
| 1676 | canvasQAv0->cd(3); |
| 1677 | pionV0plusTOFGraph->SetMarkerStyle(30); |
| 1678 | pionV0plusTOFGraph->SetMarkerColor(kMagenta); |
| 1679 | pionV0plusTOFGraph->DrawClone("p"); |
| 1680 | |
| 1681 | canvasQAtof->cd(3); |
| 1682 | pionV0plusTOFGraph->DrawClone("p"); |
| 1683 | |
| 1684 | canvasQAtpc->cd(3); |
| 1685 | pionV0plusTOFGraph->DrawClone("p"); |
| 1686 | |
| 1687 | // |
| 1688 | // 6. V0 protons + TOF |
| 1689 | // |
| 1690 | hist->GetAxis(2)->SetRangeUser(3,3); // V0 protons |
| 1691 | hist->GetAxis(3)->SetRangeUser(3,3); // protons |
| 1692 | hist->GetAxis(5)->SetRangeUser(-2.999,2.999); // 3sigma in TOF |
| 1693 | // |
| 1694 | TH2D * histProtonV0plusTOF = hist->Projection(1,0); |
| 1695 | histProtonV0plusTOF->SetName("histProtonV0plusTOF"); |
| 1696 | histProtonV0plusTOF->GetXaxis()->SetRangeUser(0.15,10.); |
| 1697 | |
| 1698 | // Remove contamination due to el and pions and low momenta because there the statistics |
| 1699 | // for the V0 protons goes down and becomes comparable with the contamination |
| 1700 | // -> This spoils the fit, if the contamination is not removed |
| 1701 | |
| 1702 | // PATTERN RECOGNITION |
| 1703 | correctionUpToBin = histProtonV0plusTOF->GetXaxis()->FindBin(0.6); |
| 1704 | |
| 1705 | for(Int_t ix = 1; ix <= correctionUpToBin; ix++) { |
| 1706 | for(Int_t jy = 1; jy <= histProtonV0plusTOF->GetYaxis()->GetNbins(); jy++) { |
| 1707 | Float_t yPos = histProtonV0plusTOF->GetYaxis()->GetBinCenter(jy); |
| 1708 | Float_t xPos = histProtonV0plusTOF->GetXaxis()->GetBinCenter(ix); |
| 1709 | Int_t bin = histProtonV0plusTOF->GetBin(ix,jy); |
| 1710 | if (yPos < betaSq.Eval(xPos)) histProtonV0plusTOF->SetBinContent(bin,0); |
| 1711 | } |
| 1712 | } |
| 1713 | |
| 1714 | FitSlicesY(histProtonV0plusTOF, heightFractionForFittingRange, cutForFitting, "QNR", &arr); |
| 1715 | TH1D * protonPointsV0plusTOF = (TH1D*) arr.At(1)->Clone(); |
| 1716 | // |
| 1717 | TGraphErrors * protonV0plusTOFGraph = new TGraphErrors(protonPointsV0plusTOF); |
| 1718 | protonV0plusTOFGraph->SetName("protonV0plusTOFGraph"); |
| 1719 | for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) { |
| 1720 | Bool_t removePoint = protonV0plusTOFGraph->GetY()[ip] < 10 || protonV0plusTOFGraph->GetEY()[ip]/protonV0plusTOFGraph->GetY()[ip] > 0.02; |
| 1721 | |
| 1722 | if (removePoint) { |
| 1723 | protonV0plusTOFGraph->RemovePoint(ip); |
| 1724 | ip--; |
| 1725 | continue; |
| 1726 | } |
| 1727 | } |
| 1728 | // |
| 1729 | canvasQAv0plusTOF->cd(2); |
| 1730 | histProtonV0plusTOF->SetTitle("V0+TOF protons"); |
| 1731 | histProtonV0plusTOF->GetXaxis()->SetTitle(momTitle.Data()); |
| 1732 | histProtonV0plusTOF->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 1733 | histProtonV0plusTOF->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 1734 | histProtonV0plusTOF->GetXaxis()->SetNoExponent(kTRUE); |
| 1735 | histProtonV0plusTOF->Draw("colz"); |
| 1736 | protonPointsV0plusTOF->SetMarkerStyle(29); |
| 1737 | protonPointsV0plusTOF->Draw("same"); |
| 1738 | protonV0Graph->DrawClone("p"); |
| 1739 | protonTofGraph->DrawClone("p"); |
| 1740 | protonTpcGraph->DrawClone("p"); |
| 1741 | |
| 1742 | canvasQAv0->cd(2); |
| 1743 | protonV0plusTOFGraph->SetMarkerStyle(30); |
| 1744 | protonV0plusTOFGraph->SetMarkerColor(kMagenta); |
| 1745 | protonV0plusTOFGraph->DrawClone("p"); |
| 1746 | |
| 1747 | canvasQAtof->cd(2); |
| 1748 | protonV0plusTOFGraph->DrawClone("p"); |
| 1749 | |
| 1750 | canvasQAtpc->cd(2); |
| 1751 | protonV0plusTOFGraph->DrawClone("p"); |
| 1752 | |
| 1753 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES |
| 1754 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES |
| 1755 | hist->GetAxis(5)->SetRange(0,-1); // RESET RANGES |
| 1756 | |
| 1757 | |
| 1758 | // Clone (V0, V0+TOF) graphs: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that |
| 1759 | // these graphs can be used for the BB fit. The original graph will keep all point in order to determine the correction and response |
| 1760 | // functions |
| 1761 | TGraphErrors * electronV0GraphForBBfit = new TGraphErrors(*electronV0Graph); |
| 1762 | electronV0GraphForBBfit->SetName("electronV0GraphForBBfit"); |
| 1763 | |
| 1764 | for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) { |
| 1765 | Bool_t removePoint = electronV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.4);// || electronV0GraphForBBfit->GetX()[ip] > 2.2; |
| 1766 | |
| 1767 | if (removePoint) { |
| 1768 | electronV0GraphForBBfit->RemovePoint(ip); |
| 1769 | ip--; |
| 1770 | continue; |
| 1771 | } |
| 1772 | } |
| 1773 | |
| 1774 | TGraphErrors * pionV0GraphForBBfit = new TGraphErrors(*pionV0Graph); |
| 1775 | pionV0GraphForBBfit->SetName("pionV0GraphForBBfit"); |
| 1776 | |
| 1777 | for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) { |
| 1778 | Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < (isPPb ? 1.0 : 0.76);// || pionV0GraphForBBfit->GetX()[ip] > 6.; //TODO NOW was < 0.4 before |
| 1779 | //Bool_t removePoint = pionV0GraphForBBfit->GetX()[ip] < 0.4 || pionV0GraphForBBfit->GetX()[ip] > 0.399;// In this case NOT used for BB fit |
| 1780 | if (removePoint) { |
| 1781 | pionV0GraphForBBfit->RemovePoint(ip); |
| 1782 | ip--; |
| 1783 | continue; |
| 1784 | } |
| 1785 | } |
| 1786 | |
| 1787 | TGraphErrors * protonV0GraphForBBfit = new TGraphErrors(*protonV0Graph); |
| 1788 | protonV0GraphForBBfit->SetName("protonV0GraphForBBfit"); |
| 1789 | |
| 1790 | for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) { |
| 1791 | Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 5.0; //TODO NOW was 2.5 before |
| 1792 | //Bool_t removePoint = protonV0GraphForBBfit->GetX()[ip] < 0.9 || protonV0GraphForBBfit->GetX()[ip] > 0.599;// In this case NOT used for BB fit |
| 1793 | if (removePoint) { |
| 1794 | protonV0GraphForBBfit->RemovePoint(ip); |
| 1795 | ip--; |
| 1796 | continue; |
| 1797 | } |
| 1798 | } |
| 1799 | |
| 1800 | |
| 1801 | // V0 plus TOF |
| 1802 | TGraphErrors * electronV0plusTOFGraphForBBfit = new TGraphErrors(*electronV0plusTOFGraph); |
| 1803 | electronV0plusTOFGraphForBBfit->SetName("electronV0plusTOFGraphForBBfit"); |
| 1804 | |
| 1805 | for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1806 | Bool_t removePoint = electronV0plusTOFGraphForBBfit->GetX()[ip] < (isPPb ? 1.2 : 0.6);//0.4 || electronV0plusTOFGraphForBBfit->GetX()[ip] > 1.3799; |
| 1807 | if (removePoint) { |
| 1808 | electronV0plusTOFGraphForBBfit->RemovePoint(ip); |
| 1809 | ip--; |
| 1810 | continue; |
| 1811 | } |
| 1812 | } |
| 1813 | |
| 1814 | TGraphErrors * pionV0plusTOFGraphForBBfit = new TGraphErrors(*pionV0plusTOFGraph); |
| 1815 | pionV0plusTOFGraphForBBfit->SetName("pionV0plusTOFGraphForBBfit"); |
| 1816 | |
| 1817 | for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1818 | Bool_t removePoint = pionV0plusTOFGraphForBBfit->GetX()[ip] < 0.4;// || pionV0plusTOFGraphForBBfit->GetX()[ip] > 6.; |
| 1819 | if (removePoint) { |
| 1820 | pionV0plusTOFGraphForBBfit->RemovePoint(ip); |
| 1821 | ip--; |
| 1822 | continue; |
| 1823 | } |
| 1824 | } |
| 1825 | |
| 1826 | TGraphErrors * protonV0plusTOFGraphForBBfit = new TGraphErrors(*protonV0plusTOFGraph); |
| 1827 | protonV0plusTOFGraphForBBfit->SetName("protonV0plusTOFGraphForBBfit"); |
| 1828 | |
| 1829 | for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1830 | Bool_t removePoint = protonV0plusTOFGraphForBBfit->GetX()[ip] < 0.9 || protonV0plusTOFGraphForBBfit->GetX()[ip] > 2.5; |
| 1831 | if (removePoint) { |
| 1832 | protonV0plusTOFGraphForBBfit->RemovePoint(ip); |
| 1833 | ip--; |
| 1834 | continue; |
| 1835 | } |
| 1836 | } |
| 1837 | |
| 1838 | |
| 1839 | // |
| 1840 | // rescale graphs to betaGamma |
| 1841 | // |
| 1842 | kaonTofGraph->SetMarkerStyle(21); |
| 1843 | kaonTpcGraph->SetMarkerStyle(22); |
| 1844 | for(Int_t ip = 0; ip < kaonTofGraph->GetN(); ip ++) { |
| 1845 | kaonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 1846 | kaonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 1847 | } |
| 1848 | for(Int_t ip = 0; ip < kaonTpcGraph->GetN(); ip ++) { |
| 1849 | kaonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 1850 | kaonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 1851 | } |
| 1852 | // |
| 1853 | electronGraph->SetMarkerStyle(22); |
| 1854 | electronV0Graph->SetMarkerStyle(34); |
| 1855 | electronV0GraphForBBfit->SetMarkerStyle(34); |
| 1856 | electronV0plusTOFGraph->SetMarkerStyle(30); |
| 1857 | electronV0plusTOFGraphForBBfit->SetMarkerStyle(30); |
| 1858 | electronGraph->SetMarkerColor(2); |
| 1859 | electronV0Graph->SetMarkerColor(2); |
| 1860 | electronV0GraphForBBfit->SetMarkerColor(2); |
| 1861 | electronV0plusTOFGraph->SetMarkerColor(2); |
| 1862 | electronV0plusTOFGraphForBBfit->SetMarkerColor(2); |
| 1863 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { |
| 1864 | electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1865 | electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1866 | } |
| 1867 | for(Int_t ip = 0; ip < electronV0Graph->GetN(); ip ++) { |
| 1868 | electronV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1869 | electronV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1870 | } |
| 1871 | for(Int_t ip = 0; ip < electronV0GraphForBBfit->GetN(); ip ++) { |
| 1872 | electronV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1873 | electronV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1874 | } |
| 1875 | for(Int_t ip = 0; ip < electronV0plusTOFGraph->GetN(); ip ++) { |
| 1876 | electronV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1877 | electronV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1878 | } |
| 1879 | for(Int_t ip = 0; ip < electronV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1880 | electronV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1881 | electronV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 1882 | } |
| 1883 | // |
| 1884 | pionTofGraph->SetMarkerStyle(21); |
| 1885 | pionTpcGraph->SetMarkerStyle(22); |
| 1886 | pionV0Graph->SetMarkerStyle(34); |
| 1887 | pionV0GraphForBBfit->SetMarkerStyle(34); |
| 1888 | pionV0plusTOFGraph->SetMarkerStyle(30); |
| 1889 | pionV0plusTOFGraphForBBfit->SetMarkerStyle(30); |
| 1890 | pionTofGraph->SetMarkerColor(3); |
| 1891 | pionTpcGraph->SetMarkerColor(3); |
| 1892 | pionV0Graph->SetMarkerColor(3); |
| 1893 | pionV0GraphForBBfit->SetMarkerColor(3); |
| 1894 | pionV0plusTOFGraph->SetMarkerColor(3); |
| 1895 | pionV0plusTOFGraphForBBfit->SetMarkerColor(3); |
| 1896 | for(Int_t ip = 0; ip < pionTofGraph->GetN(); ip ++) { |
| 1897 | pionTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1898 | pionTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1899 | } |
| 1900 | for(Int_t ip = 0; ip < pionTpcGraph->GetN(); ip ++) { |
| 1901 | pionTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1902 | pionTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1903 | } |
| 1904 | for(Int_t ip = 0; ip < pionV0Graph->GetN(); ip ++) { |
| 1905 | pionV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1906 | pionV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1907 | } |
| 1908 | for(Int_t ip = 0; ip < pionV0GraphForBBfit->GetN(); ip ++) { |
| 1909 | pionV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1910 | pionV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1911 | } |
| 1912 | for(Int_t ip = 0; ip < pionV0plusTOFGraph->GetN(); ip ++) { |
| 1913 | pionV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1914 | pionV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1915 | } |
| 1916 | for(Int_t ip = 0; ip < pionV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1917 | pionV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1918 | pionV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 1919 | } |
| 1920 | // |
| 1921 | protonTofGraph->SetMarkerStyle(21); |
| 1922 | protonTpcGraph->SetMarkerStyle(22); |
| 1923 | protonV0Graph->SetMarkerStyle(34); |
| 1924 | protonV0GraphForBBfit->SetMarkerStyle(34); |
| 1925 | protonV0plusTOFGraph->SetMarkerStyle(30); |
| 1926 | protonV0plusTOFGraphForBBfit->SetMarkerStyle(30); |
| 1927 | protonTofGraph->SetMarkerColor(4); |
| 1928 | protonTpcGraph->SetMarkerColor(4); |
| 1929 | protonV0Graph->SetMarkerColor(4); |
| 1930 | protonV0GraphForBBfit->SetMarkerColor(4); |
| 1931 | protonV0plusTOFGraph->SetMarkerColor(4); |
| 1932 | protonV0plusTOFGraphForBBfit->SetMarkerColor(4); |
| 1933 | for(Int_t ip = 0; ip < protonTofGraph->GetN(); ip ++) { |
| 1934 | protonTofGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1935 | protonTofGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1936 | } |
| 1937 | for(Int_t ip = 0; ip < protonTpcGraph->GetN(); ip ++) { |
| 1938 | protonTpcGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1939 | protonTpcGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1940 | } |
| 1941 | for(Int_t ip = 0; ip < protonV0Graph->GetN(); ip ++) { |
| 1942 | protonV0Graph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1943 | protonV0Graph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1944 | } |
| 1945 | for(Int_t ip = 0; ip < protonV0GraphForBBfit->GetN(); ip ++) { |
| 1946 | protonV0GraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1947 | protonV0GraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1948 | } |
| 1949 | for(Int_t ip = 0; ip < protonV0plusTOFGraph->GetN(); ip ++) { |
| 1950 | protonV0plusTOFGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1951 | protonV0plusTOFGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1952 | } |
| 1953 | for(Int_t ip = 0; ip < protonV0plusTOFGraphForBBfit->GetN(); ip ++) { |
| 1954 | protonV0plusTOFGraphForBBfit->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1955 | protonV0plusTOFGraphForBBfit->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 1956 | } |
| 1957 | // |
| 1958 | // |
| 1959 | TGraphErrors * graphAll = new TGraphErrors(); |
| 1960 | TList * listColl = new TList(); |
| 1961 | if (!useV0s) { |
| 1962 | //listColl->Add(protonTpcGraph); |
| 1963 | //listColl->Add(kaonTpcGraph); |
| 1964 | listColl->Add(pionTpcGraph); |
| 1965 | listColl->Add(electronGraph); |
| 1966 | listColl->Add(protonTofGraph); |
| 1967 | listColl->Add(kaonTofGraph); |
| 1968 | listColl->Add(pionTofGraph); |
| 1969 | } |
| 1970 | else { |
| 1971 | listColl->Add(pionV0GraphForBBfit); |
| 1972 | //listColl->Add(electronV0GraphForBBfit); |
| 1973 | listColl->Add(protonV0GraphForBBfit); |
| 1974 | |
| 1975 | //listColl->Add(pionV0plusTOFGraphForBBfit); |
| 1976 | listColl->Add(electronV0plusTOFGraphForBBfit); |
| 1977 | //listColl->Add(protonV0plusTOFGraphForBBfit); |
| 1978 | } |
| 1979 | MergeGraphErrors(graphAll, listColl); |
| 1980 | graphAll->SetNameTitle("beamDataPoints","beamDataPoints"); |
| 1981 | |
| 1982 | // Create graph out of TPC only, V0, V0+TOF and TPC+TOF to determine correction functions |
| 1983 | TList * listCollPr = new TList(); |
| 1984 | if (useV0s) { |
| 1985 | listCollPr->Add(protonTpcGraph); |
| 1986 | //listCollPr->Add(protonTofGraph); |
| 1987 | |
| 1988 | TGraphErrors * protonV0GraphCut = new TGraphErrors(*protonV0Graph); |
| 1989 | protonV0GraphCut->SetName("protonV0GraphCut"); |
| 1990 | for(Int_t ip = 0; ip < protonV0GraphCut->GetN(); ip ++) { |
| 1991 | Double_t mom = protonV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton); |
| 1992 | //Bool_t removePoint = mom >= 0.6 || mom < 0.35; // Will take TPC protons instead (statistics) for very low momenta |
| 1993 | Bool_t removePoint = mom < (isPPb ? 0.6 : 0.35); // Will take TPC protons instead (statistics) for very low momenta; |
| 1994 | // for pPb low statistics for the moment, therefore, go a bit further with TPC protons |
| 1995 | |
| 1996 | |
| 1997 | if (removePoint) { |
| 1998 | protonV0GraphCut->RemovePoint(ip); |
| 1999 | ip--; |
| 2000 | continue; |
| 2001 | } |
| 2002 | } |
| 2003 | /* |
| 2004 | TGraphErrors * protonV0plusTOFGraphCut = new TGraphErrors(*protonV0plusTOFGraph); |
| 2005 | protonV0plusTOFGraphCut->SetName("protonV0plusTOFGraphCut"); |
| 2006 | for(Int_t ip = 0; ip < protonV0plusTOFGraphCut->GetN(); ip ++) { |
| 2007 | Double_t mom = protonV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kProton); |
| 2008 | Bool_t removePoint = mom < 0.6; |
| 2009 | |
| 2010 | if (removePoint) { |
| 2011 | protonV0plusTOFGraphCut->RemovePoint(ip); |
| 2012 | ip--; |
| 2013 | continue; |
| 2014 | } |
| 2015 | }*/ |
| 2016 | |
| 2017 | listCollPr->Add(protonV0GraphCut); |
| 2018 | //listCollPr->Add(protonV0plusTOFGraphCut); |
| 2019 | } |
| 2020 | else { |
| 2021 | listCollPr->Add(protonTpcGraph); |
| 2022 | listCollPr->Add(protonTofGraph); |
| 2023 | } |
| 2024 | TGraphErrors * graphPrCombined = new TGraphErrors(); |
| 2025 | MergeGraphErrors(graphPrCombined, listCollPr); |
| 2026 | graphPrCombined->SetMarkerStyle(22); |
| 2027 | graphPrCombined->SetMarkerColor(4); |
| 2028 | graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined"); |
| 2029 | |
| 2030 | TList * listCollPi = new TList(); |
| 2031 | if (useV0s) { |
| 2032 | //listCollPi->Add(pionTpcGraph); |
| 2033 | //listCollPi->Add(pionTofGraph); |
| 2034 | |
| 2035 | TGraphErrors * pionV0GraphCut = new TGraphErrors(*pionV0Graph); |
| 2036 | pionV0GraphCut->SetName("pionV0GraphCut"); |
| 2037 | for(Int_t ip = 0; ip < pionV0GraphCut->GetN(); ip ++) { |
| 2038 | //Double_t mom = pionV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion); |
| 2039 | Bool_t removePoint = kFALSE;//mom >= 0.4; |
| 2040 | //|| mom < 0.2;//TODO NOW NEEDED? |
| 2041 | if (removePoint) { |
| 2042 | pionV0GraphCut->RemovePoint(ip); |
| 2043 | ip--; |
| 2044 | continue; |
| 2045 | } |
| 2046 | } |
| 2047 | /* |
| 2048 | TGraphErrors * pionV0plusTOFGraphCut = new TGraphErrors(*pionV0plusTOFGraph); |
| 2049 | pionV0plusTOFGraphCut->SetName("pionV0plusTOFGraphCut"); |
| 2050 | for(Int_t ip = 0; ip < pionV0plusTOFGraphCut->GetN(); ip ++) { |
| 2051 | Double_t mom = pionV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kPion); |
| 2052 | Bool_t removePoint = mom < 0.4; |
| 2053 | |
| 2054 | if (removePoint) { |
| 2055 | pionV0plusTOFGraphCut->RemovePoint(ip); |
| 2056 | ip--; |
| 2057 | continue; |
| 2058 | } |
| 2059 | }*/ |
| 2060 | |
| 2061 | listCollPi->Add(pionV0GraphCut); |
| 2062 | //listCollPi->Add(pionV0plusTOFGraphCut); |
| 2063 | } |
| 2064 | else { |
| 2065 | listCollPi->Add(pionTpcGraph); |
| 2066 | listCollPi->Add(pionTofGraph); |
| 2067 | } |
| 2068 | TGraphErrors * graphPiCombined = new TGraphErrors(); |
| 2069 | MergeGraphErrors(graphPiCombined, listCollPi); |
| 2070 | graphPiCombined->SetMarkerStyle(22); |
| 2071 | graphPiCombined->SetMarkerColor(3); |
| 2072 | graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined"); |
| 2073 | |
| 2074 | TList * listCollKa = new TList(); |
| 2075 | listCollKa->Add(kaonTpcGraph); |
| 2076 | listCollKa->Add(kaonTofGraph); |
| 2077 | TGraphErrors * graphKaCombined = new TGraphErrors(); |
| 2078 | MergeGraphErrors(graphKaCombined, listCollKa); |
| 2079 | graphKaCombined->SetMarkerStyle(22); |
| 2080 | graphKaCombined->SetMarkerColor(5); |
| 2081 | graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined"); |
| 2082 | |
| 2083 | TList * listCollEl = new TList(); |
| 2084 | if (useV0s) { |
| 2085 | //listCollEl->Add(electronGraph); |
| 2086 | |
| 2087 | TGraphErrors * electronV0GraphCut = new TGraphErrors(*electronV0Graph); |
| 2088 | electronV0GraphCut->SetName("electronV0GraphCut"); |
| 2089 | /* |
| 2090 | for(Int_t ip = 0; ip < electronV0GraphCut->GetN(); ip ++) { |
| 2091 | Double_t mom = electronV0GraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron); |
| 2092 | Bool_t removePoint = mom > 0.3; |
| 2093 | if (removePoint) { |
| 2094 | electronV0GraphCut->RemovePoint(ip); |
| 2095 | ip--; |
| 2096 | continue; |
| 2097 | } |
| 2098 | } |
| 2099 | |
| 2100 | TGraphErrors * electronV0plusTOFGraphCut = new TGraphErrors(*electronV0plusTOFGraph); |
| 2101 | electronV0plusTOFGraphCut->SetName("electronV0plusTOFGraphCut"); |
| 2102 | for(Int_t ip = 0; ip < electronV0plusTOFGraphCut->GetN(); ip ++) { |
| 2103 | Double_t mom = electronV0plusTOFGraphCut->GetX()[ip] * AliPID::ParticleMass(AliPID::kElectron); |
| 2104 | Bool_t removePoint = mom <= 0.4; |
| 2105 | |
| 2106 | if (removePoint) { |
| 2107 | electronV0plusTOFGraphCut->RemovePoint(ip); |
| 2108 | ip--; |
| 2109 | continue; |
| 2110 | } |
| 2111 | } |
| 2112 | */ |
| 2113 | listCollEl->Add(electronV0GraphCut); |
| 2114 | //listCollEl->Add(electronV0plusTOFGraphCut); |
| 2115 | } |
| 2116 | else { |
| 2117 | listCollEl->Add(electronGraph); |
| 2118 | } |
| 2119 | TGraphErrors * graphElCombined = new TGraphErrors(); |
| 2120 | MergeGraphErrors(graphElCombined, listCollEl); |
| 2121 | graphElCombined->SetMarkerStyle(22); |
| 2122 | graphElCombined->SetMarkerColor(2); |
| 2123 | graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined"); |
| 2124 | |
| 2125 | |
| 2126 | |
| 2127 | // |
| 2128 | // return array with all graphs |
| 2129 | // |
| 2130 | TObjArray * arrGraphs = new TObjArray(); |
| 2131 | arrGraphs->Add(graphAll); |
| 2132 | // |
| 2133 | arrGraphs->Add(electronGraph); |
| 2134 | arrGraphs->Add(electronV0Graph); |
| 2135 | arrGraphs->Add(electronV0plusTOFGraph); |
| 2136 | arrGraphs->Add(electronV0GraphForBBfit); |
| 2137 | arrGraphs->Add(pionTpcGraph); |
| 2138 | arrGraphs->Add(pionTofGraph); |
| 2139 | arrGraphs->Add(pionV0Graph); |
| 2140 | arrGraphs->Add(pionV0plusTOFGraph); |
| 2141 | arrGraphs->Add(pionV0GraphForBBfit); |
| 2142 | arrGraphs->Add(kaonTpcGraph); |
| 2143 | arrGraphs->Add(kaonTofGraph); |
| 2144 | arrGraphs->Add(protonTpcGraph); |
| 2145 | arrGraphs->Add(protonTofGraph); |
| 2146 | arrGraphs->Add(protonV0Graph); |
| 2147 | arrGraphs->Add(protonV0plusTOFGraph); |
| 2148 | arrGraphs->Add(protonV0GraphForBBfit); |
| 2149 | // |
| 2150 | arrGraphs->Add(graphElCombined); |
| 2151 | arrGraphs->Add(graphPiCombined); |
| 2152 | arrGraphs->Add(graphKaCombined); |
| 2153 | arrGraphs->Add(graphPrCombined); |
| 2154 | // |
| 2155 | // |
| 2156 | |
| 2157 | canvasQAtpc->SaveAs("splines_QA_ResidualGraphsTPC.root"); |
| 2158 | canvasQAtof->SaveAs("splines_QA_ResidualGraphsTOF.root"); |
| 2159 | canvasQAv0->SaveAs("splines_QA_ResidualGraphsV0.root"); |
| 2160 | canvasQAv0plusTOF->SaveAs("splines_QA_ResidualGraphsV0plusTOF.root"); |
| 2161 | |
| 2162 | |
| 2163 | canvasQAv0DeDxPurityEl->SaveAs("V0_dEdx_purity_el.root"); |
| 2164 | canvasQAv0DeDxPurityPi->SaveAs("V0_dEdx_purity_Pi.root"); |
| 2165 | canvasQAv0DeDxPurityPr->SaveAs("V0_dEdx_purity_Pr.root"); |
| 2166 | |
| 2167 | return arrGraphs; |
| 2168 | |
| 2169 | } |
| 2170 | |
| 2171 | |
| 2172 | //________________________________________________________________________ |
| 2173 | TObjArray * AliTPCcalibResidualPID::GetResidualGraphsMC(THnSparseF * histPidQA, const Char_t * /*system*/) { |
| 2174 | // |
| 2175 | // Extracts the residual graphs from THnSparse created from MC (NOT DATA!) |
| 2176 | // |
| 2177 | |
| 2178 | const TString momTitle = "#it{p}_{TPC} (GeV/#it{c})"; |
| 2179 | const TString dEdxTitle = "d#it{E}/d#it{x} (arb. unit)"; |
| 2180 | |
| 2181 | Int_t cutForFitting = 10; |
| 2182 | Double_t heightFractionForFittingRange = 0.1; |
| 2183 | |
| 2184 | // |
| 2185 | THnSparse * hist = histPidQA; |
| 2186 | // |
| 2187 | TCanvas * canvasQAmc = new TCanvas("canvasQAmcResGraph","Control canvas for residual graphs (MC)",100,10,1380,800); |
| 2188 | canvasQAmc->Divide(2,2); |
| 2189 | |
| 2190 | for (Int_t i = 1; i <= 4; i++) { |
| 2191 | canvasQAmc->GetPad(i)->SetGrid(1, 1); |
| 2192 | canvasQAmc->GetPad(i)->SetLogz(); |
| 2193 | canvasQAmc->GetPad(i)->SetLogx(); |
| 2194 | } |
| 2195 | // |
| 2196 | // 1. select and fit electrons |
| 2197 | // |
| 2198 | // In the following, axis 3 will also be limited in range, although nsigma itself is not used. This is to avoid multiple counting of entries |
| 2199 | hist->GetAxis(2)->SetRangeUser(0,0); // electrons |
| 2200 | hist->GetAxis(3)->SetRangeUser(0,0); // electrons (used for nsigma and for the eta correction of the signal) |
| 2201 | TH2D * histElectron = hist->Projection(1,0); |
| 2202 | histElectron->SetName("histElectron"); |
| 2203 | histElectron->RebinX(2); |
| 2204 | histElectron->GetXaxis()->SetRangeUser(0.15,8.0); |
| 2205 | TObjArray arr; |
| 2206 | FitSlicesY(histElectron, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); |
| 2207 | TH1D * electronPoints = (TH1D *) arr.At(1)->Clone(); |
| 2208 | // |
| 2209 | TGraphErrors * electronGraph = new TGraphErrors(electronPoints); |
| 2210 | electronGraph->SetName("electronGraphMC"); |
| 2211 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { |
| 2212 | Bool_t removePoint = electronGraph->GetY()[ip] < 10 || electronGraph->GetEY()[ip]/electronGraph->GetY()[ip] > 0.03 |
| 2213 | || electronGraph->GetX()[ip] < 0.15;// || electronGraph->GetX()[ip] > 3.1; |
| 2214 | if (removePoint) { |
| 2215 | electronGraph->RemovePoint(ip); |
| 2216 | ip--; |
| 2217 | continue; |
| 2218 | } |
| 2219 | } |
| 2220 | // |
| 2221 | canvasQAmc->cd(1); |
| 2222 | histElectron->SetTitle("electrons"); |
| 2223 | histElectron->GetYaxis()->SetRangeUser(50, 120); |
| 2224 | histElectron->GetXaxis()->SetTitle(momTitle.Data()); |
| 2225 | histElectron->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 2226 | histElectron->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2227 | histElectron->GetXaxis()->SetNoExponent(kTRUE); |
| 2228 | histElectron->Draw("colz"); |
| 2229 | electronPoints->SetMarkerStyle(24); |
| 2230 | electronPoints->Draw("same"); |
| 2231 | // |
| 2232 | // 2. protons |
| 2233 | // |
| 2234 | hist->GetAxis(2)->SetRangeUser(3,3); // protons |
| 2235 | hist->GetAxis(3)->SetRangeUser(3,3); // protons (used for nsigma and for the eta correction of the signal) |
| 2236 | // |
| 2237 | TH2D * histProton = hist->Projection(1,0); |
| 2238 | histProton->SetName("histProton"); |
| 2239 | histProton->GetXaxis()->SetRangeUser(0.15,8); |
| 2240 | histProton->GetYaxis()->SetRangeUser(10, hist->GetAxis(1)->GetBinUpEdge(hist->GetAxis(1)->GetNbins())); |
| 2241 | FitSlicesY(histProton, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); |
| 2242 | TH1D * protonPoints = (TH1D *) arr.At(1)->Clone(); |
| 2243 | // |
| 2244 | TGraphErrors * protonGraph = new TGraphErrors(protonPoints); |
| 2245 | protonGraph->SetName("protonGraphMC"); |
| 2246 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { |
| 2247 | Bool_t removePoint = protonGraph->GetY()[ip] < 10 || protonGraph->GetEY()[ip]/protonGraph->GetY()[ip] > 0.02 || |
| 2248 | protonGraph->GetX()[ip] < 0.15 || protonGraph->GetX()[ip] > 6; |
| 2249 | if (removePoint) { |
| 2250 | protonGraph->RemovePoint(ip); |
| 2251 | ip--; |
| 2252 | continue; |
| 2253 | } |
| 2254 | } |
| 2255 | // |
| 2256 | canvasQAmc->cd(2); |
| 2257 | histProton->SetTitle("protons"); |
| 2258 | histProton->GetXaxis()->SetTitle(momTitle.Data()); |
| 2259 | histProton->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 2260 | histProton->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2261 | histProton->GetXaxis()->SetNoExponent(kTRUE); |
| 2262 | histProton->Draw("colz"); |
| 2263 | protonPoints->SetMarkerStyle(20); |
| 2264 | protonPoints->Draw("same"); |
| 2265 | |
| 2266 | // |
| 2267 | // 3. pions |
| 2268 | // |
| 2269 | hist->GetAxis(2)->SetRangeUser(1,1); // pions |
| 2270 | hist->GetAxis(3)->SetRangeUser(1,1); // pions (used for nsigma and for the eta correction of the signal) |
| 2271 | // |
| 2272 | TH2D * histPion = hist->Projection(1,0); |
| 2273 | histPion->SetName("histPion"); |
| 2274 | histPion->GetXaxis()->SetRangeUser(0.15,50); |
| 2275 | FitSlicesY(histPion, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); |
| 2276 | TH1D * pionPoints = (TH1D *) arr.At(1)->Clone(); |
| 2277 | // |
| 2278 | TGraphErrors * pionGraph = new TGraphErrors(pionPoints); |
| 2279 | pionGraph->SetName("pionGraphMC"); |
| 2280 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { |
| 2281 | Bool_t removePoint = pionGraph->GetY()[ip] < 10 || pionGraph->GetEY()[ip]/pionGraph->GetY()[ip] > 0.005 |
| 2282 | || pionGraph->GetX()[ip] < 0.15 || pionGraph->GetX()[ip] > 30; |
| 2283 | if (removePoint) { |
| 2284 | pionGraph->RemovePoint(ip); |
| 2285 | ip--; |
| 2286 | continue; |
| 2287 | } |
| 2288 | } |
| 2289 | // |
| 2290 | canvasQAmc->cd(3); |
| 2291 | histPion->GetYaxis()->SetRangeUser(30, 90); |
| 2292 | histPion->GetXaxis()->SetTitle(momTitle.Data()); |
| 2293 | histPion->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 2294 | histPion->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2295 | histPion->GetXaxis()->SetNoExponent(kTRUE); |
| 2296 | histPion->Draw("colz"); |
| 2297 | histPion->SetTitle("pions"); |
| 2298 | pionPoints->SetMarkerStyle(20); |
| 2299 | pionPoints->Draw("same"); |
| 2300 | |
| 2301 | // |
| 2302 | // 4. kaons |
| 2303 | // |
| 2304 | hist->GetAxis(2)->SetRangeUser(2,2); // kaons |
| 2305 | hist->GetAxis(3)->SetRangeUser(2,2); // kaons (used for nsigma and for the eta correction of the signal) |
| 2306 | // |
| 2307 | TH2D * histKaon = hist->Projection(1,0); |
| 2308 | histKaon->SetName("histKaon"); |
| 2309 | histKaon->GetXaxis()->SetRangeUser(0.15,8); |
| 2310 | FitSlicesY(histKaon, heightFractionForFittingRange, cutForFitting, "QNRL", &arr); |
| 2311 | TH1D * kaonPoints = (TH1D*) arr.At(1)->Clone(); |
| 2312 | // |
| 2313 | TGraphErrors * kaonGraph = new TGraphErrors(kaonPoints); |
| 2314 | kaonGraph->SetName("kaonGraphMC"); |
| 2315 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { |
| 2316 | Bool_t removePoint = kaonGraph->GetY()[ip] < 10 || kaonGraph->GetEY()[ip]/kaonGraph->GetY()[ip] > 0.02 |
| 2317 | || kaonGraph->GetX()[ip] < 0.15; |
| 2318 | if (removePoint) { |
| 2319 | kaonGraph->RemovePoint(ip); |
| 2320 | ip--; |
| 2321 | continue; |
| 2322 | } |
| 2323 | } |
| 2324 | // |
| 2325 | canvasQAmc->cd(4); |
| 2326 | histKaon->SetTitle("kaons"); |
| 2327 | histKaon->GetYaxis()->SetRangeUser(30, 500); |
| 2328 | histKaon->GetXaxis()->SetTitle(momTitle.Data()); |
| 2329 | histKaon->GetYaxis()->SetTitle(dEdxTitle.Data()); |
| 2330 | histKaon->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2331 | histKaon->GetXaxis()->SetNoExponent(kTRUE); |
| 2332 | histKaon->Draw("colz"); |
| 2333 | kaonPoints->SetMarkerStyle(20); |
| 2334 | kaonPoints->Draw("same"); |
| 2335 | |
| 2336 | |
| 2337 | |
| 2338 | hist->GetAxis(2)->SetRange(0,-1); // RESET RANGES |
| 2339 | hist->GetAxis(3)->SetRange(0,-1); // RESET RANGES |
| 2340 | |
| 2341 | |
| 2342 | // Clone graphs to have graphs with the same names as for the data case -> same subsequent functions can be used to determine |
| 2343 | // the correction functions and, finally, the response functions. |
| 2344 | // In addition: Remove some points, which are expected to deviate from the Bethe-Bloch curve, from the original graphs, so that |
| 2345 | // these graphs can be used for the BB fit |
| 2346 | TGraphErrors * graphPrCombined = new TGraphErrors(*protonGraph); |
| 2347 | graphPrCombined->SetMarkerStyle(22); |
| 2348 | graphPrCombined->SetMarkerColor(4); |
| 2349 | graphPrCombined->SetNameTitle("Graph_Protons_Combined", "Graph_Protons_Combined"); |
| 2350 | |
| 2351 | TGraphErrors * graphPiCombined = new TGraphErrors(*pionGraph); |
| 2352 | graphPiCombined->SetMarkerStyle(22); |
| 2353 | graphPiCombined->SetMarkerColor(3); |
| 2354 | graphPiCombined->SetNameTitle("Graph_Pions_Combined", "Graph_Pions_Combined"); |
| 2355 | |
| 2356 | TGraphErrors * graphKaCombined = new TGraphErrors(*kaonGraph); |
| 2357 | graphKaCombined->SetMarkerStyle(22); |
| 2358 | graphKaCombined->SetMarkerColor(5); |
| 2359 | graphKaCombined->SetNameTitle("Graph_Kaons_Combined", "Graph_Kaons_Combined"); |
| 2360 | |
| 2361 | TGraphErrors * graphElCombined = new TGraphErrors(*electronGraph); |
| 2362 | graphElCombined->SetMarkerStyle(22); |
| 2363 | graphElCombined->SetMarkerColor(2); |
| 2364 | graphElCombined->SetNameTitle("Graph_Electrons_Combined", "Graph_Electrons_Combined"); |
| 2365 | |
| 2366 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { |
| 2367 | Bool_t removePoint = electronGraph->GetX()[ip] < 2.5;// || electronGraph->GetX()[ip] > 5.0; |
| 2368 | //Bool_t removePoint = electronGraph->GetX()[ip] < 1.2 || electronGraph->GetX()[ip] > 5.0; |
| 2369 | if (removePoint) { |
| 2370 | electronGraph->RemovePoint(ip); |
| 2371 | ip--; |
| 2372 | continue; |
| 2373 | } |
| 2374 | } |
| 2375 | |
| 2376 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { |
| 2377 | Bool_t removePoint = protonGraph->GetX()[ip] < 0.9 |
| 2378 | || protonGraph->GetX()[ip] > 2.5;//3.5; //TODO NEW in order not to get into conflict with pions/kaons |
| 2379 | if (removePoint) { |
| 2380 | protonGraph->RemovePoint(ip); |
| 2381 | ip--; |
| 2382 | continue; |
| 2383 | } |
| 2384 | } |
| 2385 | |
| 2386 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { |
| 2387 | Bool_t removePoint = pionGraph->GetX()[ip] < 2.1;//TODO APR18 0.8; |
| 2388 | if (removePoint) { |
| 2389 | pionGraph->RemovePoint(ip); |
| 2390 | ip--; |
| 2391 | continue; |
| 2392 | } |
| 2393 | } |
| 2394 | |
| 2395 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { |
| 2396 | Bool_t removePoint = kaonGraph->GetX()[ip] < 1.25 || kaonGraph->GetX()[ip] > 4.0;// < 0.7;// || kaonGraph->GetX()[ip] > 4; |
| 2397 | if (removePoint) { |
| 2398 | kaonGraph->RemovePoint(ip); |
| 2399 | ip--; |
| 2400 | continue; |
| 2401 | } |
| 2402 | } |
| 2403 | // |
| 2404 | // rescale graphs to betaGamma |
| 2405 | // |
| 2406 | kaonGraph->SetMarkerStyle(22); |
| 2407 | for(Int_t ip = 0; ip < kaonGraph->GetN(); ip ++) { |
| 2408 | kaonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 2409 | kaonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 2410 | } |
| 2411 | for(Int_t ip = 0; ip < graphKaCombined->GetN(); ip ++) { |
| 2412 | graphKaCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 2413 | graphKaCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kKaon); |
| 2414 | } |
| 2415 | // |
| 2416 | electronGraph->SetMarkerStyle(22); |
| 2417 | electronGraph->SetMarkerColor(2); |
| 2418 | for(Int_t ip = 0; ip < electronGraph->GetN(); ip ++) { |
| 2419 | electronGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 2420 | electronGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 2421 | } |
| 2422 | for(Int_t ip = 0; ip < graphElCombined->GetN(); ip ++) { |
| 2423 | graphElCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 2424 | graphElCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kElectron); |
| 2425 | } |
| 2426 | // |
| 2427 | pionGraph->SetMarkerStyle(22); |
| 2428 | pionGraph->SetMarkerColor(3); |
| 2429 | for(Int_t ip = 0; ip < pionGraph->GetN(); ip ++) { |
| 2430 | pionGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 2431 | pionGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 2432 | } |
| 2433 | for(Int_t ip = 0; ip < graphPiCombined->GetN(); ip ++) { |
| 2434 | graphPiCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 2435 | graphPiCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kPion); |
| 2436 | } |
| 2437 | // |
| 2438 | protonGraph->SetMarkerStyle(22); |
| 2439 | protonGraph->SetMarkerColor(4); |
| 2440 | for(Int_t ip = 0; ip < protonGraph->GetN(); ip ++) { |
| 2441 | protonGraph->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 2442 | protonGraph->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 2443 | } |
| 2444 | for(Int_t ip = 0; ip < graphPrCombined->GetN(); ip ++) { |
| 2445 | graphPrCombined->GetX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 2446 | graphPrCombined->GetEX()[ip] /= AliPID::ParticleMass(AliPID::kProton); |
| 2447 | } |
| 2448 | // |
| 2449 | // |
| 2450 | |
| 2451 | // Store graphs without further restriction |
| 2452 | TGraphErrors * graphPrAll = new TGraphErrors(*graphPrCombined); |
| 2453 | graphPrAll->SetNameTitle("Protons_MC_all", "Protons_MC_all"); |
| 2454 | |
| 2455 | TGraphErrors * graphPiAll = new TGraphErrors(*graphPiCombined); |
| 2456 | graphPiAll->SetNameTitle("Pions_MC_all", "Pions_MC_all"); |
| 2457 | |
| 2458 | TGraphErrors * graphKaAll = new TGraphErrors(*graphKaCombined); |
| 2459 | graphKaAll->SetNameTitle("Kaons_MC_all", "Kaons_MC_all"); |
| 2460 | |
| 2461 | TGraphErrors * graphElAll = new TGraphErrors(*graphElCombined); |
| 2462 | graphElAll->SetNameTitle("Electrons_MC_all", "Electrons_MC_all"); |
| 2463 | |
| 2464 | // |
| 2465 | // |
| 2466 | TGraphErrors * graphAll = new TGraphErrors(); |
| 2467 | TList * listColl = new TList(); |
| 2468 | listColl->Add(electronGraph); |
| 2469 | listColl->Add(protonGraph); |
| 2470 | listColl->Add(kaonGraph); |
| 2471 | listColl->Add(pionGraph); |
| 2472 | MergeGraphErrors(graphAll, listColl); |
| 2473 | graphAll->SetNameTitle("beamDataPoints","beamDataPoints"); |
| 2474 | |
| 2475 | // |
| 2476 | // return array with all graphs |
| 2477 | // |
| 2478 | TObjArray * arrGraphs = new TObjArray(); |
| 2479 | arrGraphs->Add(graphAll); |
| 2480 | // |
| 2481 | arrGraphs->Add(electronGraph); |
| 2482 | arrGraphs->Add(pionGraph); |
| 2483 | arrGraphs->Add(kaonGraph); |
| 2484 | arrGraphs->Add(protonGraph); |
| 2485 | // |
| 2486 | arrGraphs->Add(graphElAll); |
| 2487 | arrGraphs->Add(graphPiAll); |
| 2488 | arrGraphs->Add(graphKaAll); |
| 2489 | arrGraphs->Add(graphPrAll); |
| 2490 | // |
| 2491 | arrGraphs->Add(graphElCombined); |
| 2492 | arrGraphs->Add(graphPiCombined); |
| 2493 | arrGraphs->Add(graphKaCombined); |
| 2494 | arrGraphs->Add(graphPrCombined); |
| 2495 | // |
| 2496 | // |
| 2497 | |
| 2498 | canvasQAmc->SaveAs("splines_QA_ResidualGraphsTPC.root"); |
| 2499 | |
| 2500 | return arrGraphs; |
| 2501 | |
| 2502 | } |
| 2503 | |
| 2504 | |
| 2505 | //________________________________________________________________________ |
| 2506 | TObjArray * AliTPCcalibResidualPID::GetResponseFunctions(TF1* parametrisation, TObjArray* inputGraphs, const Char_t * type, const Char_t * period, const Char_t * pass, const Char_t * system, const Char_t* dedxtype) { |
| 2507 | // |
| 2508 | // |
| 2509 | |
| 2510 | Bool_t isMC = kFALSE; |
| 2511 | if (strcmp(type, "MC") == 0) { |
| 2512 | isMC = kTRUE; |
| 2513 | } |
| 2514 | else if (strcmp(type, "DATA") == 0) { |
| 2515 | isMC = kFALSE; |
| 2516 | } |
| 2517 | else { |
| 2518 | Printf("ERROR - GetResponseFunctions: Unknown type \"%s\" - must be \"MC\" or \"DATA\"!", type); |
| 2519 | |
| 2520 | return 0x0; |
| 2521 | } |
| 2522 | |
| 2523 | Bool_t isPbPb = kFALSE; |
| 2524 | Bool_t isPPb = kFALSE; |
| 2525 | if (strcmp(system, "PBPB") == 0) { |
| 2526 | Printf("PbPb detected - Applying special handling!"); |
| 2527 | isPbPb = kTRUE; |
| 2528 | } |
| 2529 | else if (strcmp(system, "PPB") == 0 || strcmp(system, "PBP") == 0) { |
| 2530 | Printf("p-Pb/Pb-p detected - Applying special handling!"); |
| 2531 | isPPb = kTRUE; |
| 2532 | } |
| 2533 | |
| 2534 | TCanvas* canvasQA[4] = { 0x0, }; |
| 2535 | for (Int_t i = 0; i < 4; i++) { |
| 2536 | canvasQA[i] = new TCanvas(Form("canvasQAres_%d", i), "Control canvas for residual polynomial",100,10,1380,800); |
| 2537 | canvasQA[i]->SetGrid(1, 1); |
| 2538 | canvasQA[i]->SetLogx(); |
| 2539 | } |
| 2540 | /* |
| 2541 | TCanvas * canvasQA = new TCanvas("canvasQAres","Control canvas for residual polynomial",100,10,1380,800); |
| 2542 | canvasQA->Divide(2,2); |
| 2543 | |
| 2544 | for (Int_t i = 1; i <= 4; i++) { |
| 2545 | canvasQA->GetPad(i)->SetGrid(1, 1); |
| 2546 | } |
| 2547 | */ |
| 2548 | // |
| 2549 | // smallest needed (100MeV proton) bg = 0.1/0.938(=AliPID::ParticleMass(AliPID::kProton)) = 0.1, |
| 2550 | // largest needed 100 GeV electron, bg = 100/0.000511(=AliPID::ParticleMass(AliPID::kElectron)) = 200 000 |
| 2551 | // |
| 2552 | Double_t from = 0.1; |
| 2553 | Double_t to = 200000; |
| 2554 | |
| 2555 | TF1* funcBB = new TF1(*parametrisation); |
| 2556 | funcBB->SetLineColor(2); |
| 2557 | funcBB->SetLineWidth(1); |
| 2558 | |
| 2559 | TString fitFuncString = "[0] + [1] / x + [2] / (x**2) + [3] / (x**3) + [4] / (x**4) + [5] / (x**5) + [6] * x"; |
| 2560 | TString fitFuncString2 = "pol6"; |
| 2561 | // |
| 2562 | // 1. extract proton corrections |
| 2563 | // |
| 2564 | TGraphErrors * graphProtonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Protons_Combined"); |
| 2565 | TGraphErrors * graphProtonTPCfinal = new TGraphErrors(*graphProtonTPC); |
| 2566 | graphProtonTPCfinal->SetName("graphProtonTPCfinal"); |
| 2567 | for(Int_t ip = 0; ip < graphProtonTPC->GetN(); ip ++) { |
| 2568 | graphProtonTPC->GetY()[ip] -= funcBB->Eval(graphProtonTPC->GetX()[ip]); |
| 2569 | graphProtonTPC->GetY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]); |
| 2570 | graphProtonTPC->GetEY()[ip] /= funcBB->Eval(graphProtonTPC->GetX()[ip]);; |
| 2571 | } |
| 2572 | canvasQA[0]->cd(); |
| 2573 | //canvasQA->cd(1); |
| 2574 | //gPad->SetLogx(); |
| 2575 | graphProtonTPC->GetXaxis()->SetTitle("#beta#gamma"); |
| 2576 | graphProtonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2577 | graphProtonTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 2578 | graphProtonTPC->GetYaxis()->SetLabelSize(0.04); |
| 2579 | graphProtonTPC->GetYaxis()->SetTitleOffset(0.85); |
| 2580 | graphProtonTPC->GetXaxis()->SetTitleOffset(1.0); |
| 2581 | graphProtonTPC->GetYaxis()->SetTitle("(data - fit) / fit"); |
| 2582 | graphProtonTPC->Draw("ap"); |
| 2583 | |
| 2584 | if (graphProtonTPC->GetN() <= 0) { |
| 2585 | Printf("ERROR - GetResponseFunctions: Proton graph has no data points!"); |
| 2586 | |
| 2587 | return 0x0; |
| 2588 | } |
| 2589 | graphProtonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function |
| 2590 | TF1 * funcCorrProton = new TF1("funcCorrProton", fitFuncString2.Data(), TMath::Max(graphProtonTPC->GetX()[0], 0.15),1.0); // TODO BEN was 0.18 - 0.85 |
| 2591 | graphProtonTPC->Fit(funcCorrProton, "QREX0M"); |
| 2592 | // |
| 2593 | // 2. extract kaon corrections |
| 2594 | // |
| 2595 | TGraphErrors * graphKaonTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Kaons_Combined"); |
| 2596 | TGraphErrors * graphKaonTPCfinal = new TGraphErrors(*graphKaonTPC); |
| 2597 | graphKaonTPCfinal->SetName("graphKaonTPCfinal"); |
| 2598 | for(Int_t ip = 0; ip < graphKaonTPC->GetN(); ip ++) { |
| 2599 | graphKaonTPC->GetY()[ip] -= funcBB->Eval(graphKaonTPC->GetX()[ip]); |
| 2600 | graphKaonTPC->GetY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]); |
| 2601 | graphKaonTPC->GetEY()[ip] /= funcBB->Eval(graphKaonTPC->GetX()[ip]);; |
| 2602 | } |
| 2603 | canvasQA[1]->cd(); |
| 2604 | //canvasQA->cd(2); |
| 2605 | //gPad->SetLogx(); |
| 2606 | graphKaonTPC->GetXaxis()->SetTitle("#beta#gamma"); |
| 2607 | graphKaonTPC->GetYaxis()->SetTitle("(data - fit) / fit"); |
| 2608 | graphKaonTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2609 | graphKaonTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 2610 | graphKaonTPC->GetYaxis()->SetLabelSize(0.04); |
| 2611 | graphKaonTPC->GetYaxis()->SetTitleOffset(0.85); |
| 2612 | graphKaonTPC->GetXaxis()->SetTitleOffset(1.0); |
| 2613 | graphKaonTPC->Draw("ap"); |
| 2614 | |
| 2615 | if (graphKaonTPC->GetN() <= 0) { |
| 2616 | Printf("ERROR - GetResponseFunctions: Kaon graph has no data points!"); |
| 2617 | |
| 2618 | return 0x0; |
| 2619 | } |
| 2620 | graphKaonTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function |
| 2621 | |
| 2622 | TF1 * funcCorrKaon = 0x0; |
| 2623 | |
| 2624 | if (isMC) { |
| 2625 | funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), isPPb ? 3.44 : 1.981); |
| 2626 | graphKaonTPC->Fit(funcCorrKaon, "QREX0M"); |
| 2627 | } |
| 2628 | else { |
| 2629 | // In case of data there are sometimes problems to fit the shape with one function (could describe the overall deviation, |
| 2630 | // but will not cover all the details). |
| 2631 | // Nevertheless, this shape (including most of the "details") can be fitted with the following approach with two functions |
| 2632 | TF1 * funcCorrKaon1 = new TF1("funcCorrKaon1", fitFuncString.Data(), |
| 2633 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); |
| 2634 | graphKaonTPC->Fit(funcCorrKaon1, "QREX0M", "same", TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.0); |
| 2635 | |
| 2636 | TF1 * funcCorrKaon2 = new TF1("funcCorrKaon2", fitFuncString2.Data(), TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); |
| 2637 | graphKaonTPC->Fit(funcCorrKaon2, "QREX0M", "same", (isMC ? 1.981 : 1.0), 1.981); |
| 2638 | |
| 2639 | funcCorrKaon = new TF1("funcCorrKaon", |
| 2640 | "funcCorrKaon1 * 0.5*(1.+TMath::Erf((1 - x) / 0.1)) + ([0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x+[5]*x*x*x*x*x+[6]*x*x*x*x*x*x) * 0.5*(1.+TMath::Erf((x - 1) / 0.1))", |
| 2641 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), 1.981); |
| 2642 | |
| 2643 | for (Int_t i = funcCorrKaon1->GetNpar(), j = 0; i < funcCorrKaon1->GetNpar() + funcCorrKaon2->GetNpar(); i++, j++) { |
| 2644 | funcCorrKaon->SetParameter(j, funcCorrKaon2->GetParameter(j)); |
| 2645 | } |
| 2646 | funcCorrKaon->SetLineColor(kRed); |
| 2647 | funcCorrKaon->GetHistogram()->DrawClone("csame"); |
| 2648 | //funcCorrKaon->Draw("same"); |
| 2649 | } |
| 2650 | /*TODO |
| 2651 | TF1 * funcCorrKaon = new TF1("funcCorrKaon", fitFuncString.Data(),//TODO BEN was fitFuncString2 |
| 2652 | TMath::Max(graphKaonTPC->GetX()[0], 0.25), (isMC ? 1.981 : 1.0)); //TODO BEN was 0.79 for data and 1.45 for MC |
| 2653 | graphKaonTPC->Fit(funcCorrKaon, "QREX0M"); |
| 2654 | */ |
| 2655 | // |
| 2656 | // 3. extract pion corrections |
| 2657 | // |
| 2658 | TGraphErrors * graphPionTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Pions_Combined"); |
| 2659 | TGraphErrors * graphPionTPCfinal = new TGraphErrors(*graphPionTPC); |
| 2660 | graphPionTPCfinal->SetName("graphPionTPCfinal"); |
| 2661 | for(Int_t ip = 0; ip < graphPionTPC->GetN(); ip ++) { |
| 2662 | graphPionTPC->GetY()[ip] -= funcBB->Eval(graphPionTPC->GetX()[ip]); |
| 2663 | graphPionTPC->GetY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]); |
| 2664 | graphPionTPC->GetEY()[ip] /= funcBB->Eval(graphPionTPC->GetX()[ip]); |
| 2665 | } |
| 2666 | canvasQA[2]->cd(); |
| 2667 | //canvasQA->cd(3); |
| 2668 | //gPad->SetLogx(); |
| 2669 | graphPionTPC->GetXaxis()->SetTitle("#beta#gamma"); |
| 2670 | graphPionTPC->GetYaxis()->SetTitle("(data - fit) / fit"); |
| 2671 | graphPionTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2672 | graphPionTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 2673 | graphPionTPC->GetYaxis()->SetLabelSize(0.04); |
| 2674 | graphPionTPC->GetYaxis()->SetTitleOffset(0.85); |
| 2675 | graphPionTPC->GetXaxis()->SetTitleOffset(1.0); |
| 2676 | graphPionTPC->Draw("ap"); |
| 2677 | |
| 2678 | if (graphPionTPC->GetN() <= 0) { |
| 2679 | Printf("ERROR - GetResponseFunctions: Pion graph has no data points!"); |
| 2680 | |
| 2681 | return 0x0; |
| 2682 | } |
| 2683 | graphPionTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function |
| 2684 | // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to electrons |
| 2685 | TF1 * funcCorrPion = new TF1("funcCorrPion", fitFuncString.Data(), ((isPbPb && !isMC) ? (0.2 / AliPID::ParticleMass(AliPID::kPion)) : graphPionTPC->GetX()[0]), isMC ? (isPPb ? 12.8 : (isPbPb ? 12.8 : 7.1)) : (isPPb ? 7.68 : 7.1)); |
| 2686 | graphPionTPC->Fit(funcCorrPion, "QREX0M"); |
| 2687 | // |
| 2688 | // 4. extract electron corrections |
| 2689 | // |
| 2690 | TGraphErrors * graphElectronTPC = (TGraphErrors *) inputGraphs->FindObject("Graph_Electrons_Combined"); |
| 2691 | TGraphErrors * graphElectronTPCfinal = new TGraphErrors(*graphElectronTPC); |
| 2692 | graphElectronTPCfinal->SetName("graphElectronTPCfinal"); |
| 2693 | for(Int_t ip = 0; ip < graphElectronTPC->GetN(); ip ++) { |
| 2694 | graphElectronTPC->GetY()[ip] -= funcBB->Eval(graphElectronTPC->GetX()[ip]); |
| 2695 | graphElectronTPC->GetY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]); |
| 2696 | graphElectronTPC->GetEY()[ip] /= funcBB->Eval(graphElectronTPC->GetX()[ip]);; |
| 2697 | } |
| 2698 | canvasQA[3]->cd(); |
| 2699 | //canvasQA->cd(4); |
| 2700 | //gPad->SetLogx(); |
| 2701 | graphElectronTPC->GetXaxis()->SetTitle("#beta#gamma"); |
| 2702 | graphElectronTPC->GetYaxis()->SetTitle("(data - fit) / fit"); |
| 2703 | graphElectronTPC->GetXaxis()->SetMoreLogLabels(kTRUE); |
| 2704 | graphElectronTPC->GetXaxis()->SetNoExponent(kTRUE); |
| 2705 | graphElectronTPC->GetYaxis()->SetLabelSize(0.04); |
| 2706 | graphElectronTPC->GetYaxis()->SetTitleOffset(0.85); |
| 2707 | graphElectronTPC->GetXaxis()->SetTitleOffset(1.0); |
| 2708 | graphElectronTPC->Draw("ap"); |
| 2709 | |
| 2710 | if (graphElectronTPC->GetN() <= 0) { |
| 2711 | Printf("ERROR - GetResponseFunctions: Electron graph has no data points!"); |
| 2712 | |
| 2713 | return 0x0; |
| 2714 | } |
| 2715 | graphElectronTPC->Sort(); // Sort points along x. Next, the very first point will be used to determine the starting point of the correction function |
| 2716 | // In case of PbPb (data, not MC) only correct down to 0.2 GeV/c; otherwise: contamination due to pions |
| 2717 | TF1 * funcCorrElectron = new TF1("funcCorrElectron", fitFuncString.Data(), (!isMC && isPbPb ? (0.2 / AliPID::ParticleMass(AliPID::kElectron)) :graphElectronTPC->GetX()[0]), (isMC ? 3565 : (isPPb ? 2900 : 1920/*970*/)));// TODO was 1800 for pp data |
| 2718 | // NOTE: For data, the results are almost the same for fitFuncString and fitFuncString2. Maybe, fitFuncString2 is slightly better. |
| 2719 | graphElectronTPC->Fit(funcCorrElectron, "QREX0M"); |
| 2720 | // |
| 2721 | // EXTRACT GRAPHS AND PUT THEM TO THE TOBJARRAY |
| 2722 | // |
| 2723 | const Int_t nBins = 500; |
| 2724 | Double_t xBetaGamma[nBins]; |
| 2725 | Double_t yProton[nBins]; |
| 2726 | Double_t yKaon[nBins]; |
| 2727 | Double_t yPion[nBins]; |
| 2728 | Double_t yElectron[nBins]; |
| 2729 | Double_t yDefault[nBins]; |
| 2730 | // |
| 2731 | // |
| 2732 | // |
| 2733 | xBetaGamma[0] = from; |
| 2734 | Double_t factor = pow(to/from, 1./nBins); |
| 2735 | |
| 2736 | // |
| 2737 | for(Int_t kk = 0; kk < nBins; kk++) { |
| 2738 | if (kk > 0) xBetaGamma[kk] = factor * xBetaGamma[kk-1]; |
| 2739 | yProton[kk] = funcBB->Eval(xBetaGamma[kk])/50.; |
| 2740 | yPion[kk] = funcBB->Eval(xBetaGamma[kk])/50.; |
| 2741 | yKaon[kk] = funcBB->Eval(xBetaGamma[kk])/50.; |
| 2742 | yElectron[kk] = funcBB->Eval(xBetaGamma[kk])/50.; |
| 2743 | yDefault[kk] = funcBB->Eval(xBetaGamma[kk])/50.; |
| 2744 | |
| 2745 | // Added by Ben |
| 2746 | Double_t widthFactor = 0.020; |
| 2747 | Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kProton)/widthFactor) + 1); |
| 2748 | Double_t smoothKaon = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kKaon)/widthFactor) + 1); |
| 2749 | Double_t smoothPion = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kPion)/widthFactor) + 1); |
| 2750 | Double_t smoothElectron = 0.5*(TMath::Erf((funcCorrElectron->GetXmax()-xBetaGamma[kk])*AliPID::ParticleMass(AliPID::kElectron)/widthFactor) + 1); |
| 2751 | |
| 2752 | if (xBetaGamma[kk] > funcCorrProton->GetXmax()) |
| 2753 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmax())); |
| 2754 | // Correction is so large that one runs into trouble at low bg for Protons. |
| 2755 | // The deviation is smaller if one takes the lower bound of the correction function |
| 2756 | else if (xBetaGamma[kk] < funcCorrProton->GetXmin()) |
| 2757 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(funcCorrProton->GetXmin())); |
| 2758 | else |
| 2759 | yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk])); |
| 2760 | |
| 2761 | if (xBetaGamma[kk] > funcCorrKaon->GetXmax()) |
| 2762 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmax())); |
| 2763 | // Correction is so large that one runs into trouble at low bg for Kaons. |
| 2764 | // The deviation is smaller if one takes the lower bound of the correction function |
| 2765 | else if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) |
| 2766 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(funcCorrKaon->GetXmin())); |
| 2767 | else |
| 2768 | yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk])); |
| 2769 | |
| 2770 | if (xBetaGamma[kk] > funcCorrElectron->GetXmax()) |
| 2771 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmax())); |
| 2772 | else if (xBetaGamma[kk] < funcCorrElectron->GetXmin()) |
| 2773 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(funcCorrElectron->GetXmin())); |
| 2774 | else |
| 2775 | yElectron[kk] *= (1 + smoothElectron*funcCorrElectron->Eval(xBetaGamma[kk])); |
| 2776 | |
| 2777 | // Only true for LHC10d.pass?: Seems not to be needed because any (very small) deviations are most likely due to electron contamination(BEN) |
| 2778 | if (xBetaGamma[kk] > funcCorrPion->GetXmax()) |
| 2779 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmax())); |
| 2780 | else if (xBetaGamma[kk] < funcCorrPion->GetXmin()) |
| 2781 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(funcCorrPion->GetXmin())); |
| 2782 | else |
| 2783 | yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk])); |
| 2784 | |
| 2785 | |
| 2786 | /* Removed by Ben |
| 2787 | Double_t smoothProton = 0.5*(TMath::Erf((funcCorrProton->GetXmax()-xBetaGamma[kk])/0.002) + 1); |
| 2788 | Double_t smoothKaon = 0.5*(TMath::Erf((funcCorrKaon->GetXmax()-xBetaGamma[kk])/0.002) + 1); |
| 2789 | Double_t smoothPion = 0.5*(TMath::Erf((funcCorrPion->GetXmax()-xBetaGamma[kk])/0.002) + 1); |
| 2790 | |
| 2791 | if (xBetaGamma[kk] < funcCorrProton->GetXmax() && xBetaGamma[kk] > funcCorrProton->GetXmin()) yProton[kk] *= (1 + smoothProton*funcCorrProton->Eval(xBetaGamma[kk])); |
| 2792 | if (xBetaGamma[kk] < funcCorrKaon->GetXmax() && xBetaGamma[kk] > funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + smoothKaon*funcCorrKaon->Eval(xBetaGamma[kk])); |
| 2793 | if (xBetaGamma[kk] < funcCorrPion->GetXmax() && xBetaGamma[kk] > funcCorrPion->GetXmin()) yPion[kk] *= (1 + smoothPion*funcCorrPion->Eval(xBetaGamma[kk])); |
| 2794 | //if (xBetaGamma[kk] < funcCorrElectron->GetXmax()) yElectron[kk] *= (1 + funcCorrElectron->Eval(xBetaGamma[kk])); |
| 2795 | // |
| 2796 | if (xBetaGamma[kk] < funcCorrProton->GetXmin()) yProton[kk] *= (1 + funcCorrProton->Eval(funcCorrProton->GetXmin())); |
| 2797 | if (xBetaGamma[kk] < funcCorrKaon->GetXmin()) yKaon[kk] *= (1 + funcCorrKaon->Eval(funcCorrKaon->GetXmin())); |
| 2798 | if (xBetaGamma[kk] < funcCorrPion->GetXmin()) yPion[kk] *= (1 + funcCorrPion->Eval(funcCorrPion->GetXmin())); |
| 2799 | */ |
| 2800 | } |
| 2801 | // |
| 2802 | TGraph * graphProton = new TGraph(nBins, xBetaGamma, yProton); |
| 2803 | TGraph * graphElectron = new TGraph(nBins, xBetaGamma, yElectron); |
| 2804 | TGraph * graphKaon = new TGraph(nBins, xBetaGamma, yKaon); |
| 2805 | TGraph * graphPion = new TGraph(nBins, xBetaGamma, yPion); |
| 2806 | TGraph * graphDefault = new TGraph(nBins, xBetaGamma, yDefault); |
| 2807 | // |
| 2808 | // |
| 2809 | TSpline3 * splineProton = new TSpline3("splineProton", graphProton); |
| 2810 | TSpline3 * splineElectron = new TSpline3("splineElectron", graphElectron); |
| 2811 | TSpline3 * splineKaon = new TSpline3("splineKaon", graphKaon); |
| 2812 | TSpline3 * splinePion = new TSpline3("splinePion", graphPion); |
| 2813 | TSpline3 * splineDefault = new TSpline3("splineDefault", graphDefault); |
| 2814 | // |
| 2815 | // |
| 2816 | splineProton->SetNameTitle(Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), |
| 2817 | Form("TSPLINE3_%s_PROTON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); |
| 2818 | splineElectron->SetNameTitle(Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), |
| 2819 | Form("TSPLINE3_%s_ELECTRON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); |
| 2820 | splineKaon->SetNameTitle(Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), |
| 2821 | Form("TSPLINE3_%s_KAON_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); |
| 2822 | splinePion->SetNameTitle(Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), |
| 2823 | Form("TSPLINE3_%s_PION_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); |
| 2824 | splineDefault->SetNameTitle(Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype), |
| 2825 | Form("TSPLINE3_%s_ALL_%s_%s_%s_%sMEAN",type,period,pass,system,dedxtype)); |
| 2826 | // |
| 2827 | TObjArray * arrResponse = new TObjArray(); |
| 2828 | arrResponse->SetName("TPCpidResponseFunctions"); |
| 2829 | arrResponse->AddLast(splineProton); |
| 2830 | arrResponse->AddLast(splineElectron); |
| 2831 | arrResponse->AddLast(splineKaon); |
| 2832 | arrResponse->AddLast(splinePion); |
| 2833 | arrResponse->AddLast(splineDefault); |
| 2834 | // |
| 2835 | |
| 2836 | // Draw deviation from final results |
| 2837 | for(Int_t ip = 0; ip < graphProtonTPCfinal->GetN(); ip ++) { |
| 2838 | graphProtonTPCfinal->GetY()[ip] -= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]); |
| 2839 | graphProtonTPCfinal->GetY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]); |
| 2840 | graphProtonTPCfinal->GetEY()[ip] /= 50. * splineProton->Eval(graphProtonTPCfinal->GetX()[ip]);; |
| 2841 | } |
| 2842 | canvasQA[0]->cd(); |
| 2843 | //canvasQA->cd(1); |
| 2844 | graphProtonTPCfinal->SetMarkerStyle(26); |
| 2845 | graphProtonTPCfinal->Draw("psame"); |
| 2846 | |
| 2847 | for(Int_t ip = 0; ip < graphKaonTPCfinal->GetN(); ip ++) { |
| 2848 | graphKaonTPCfinal->GetY()[ip] -= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]); |
| 2849 | graphKaonTPCfinal->GetY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]); |
| 2850 | graphKaonTPCfinal->GetEY()[ip] /= 50. * splineKaon->Eval(graphKaonTPCfinal->GetX()[ip]);; |
| 2851 | } |
| 2852 | canvasQA[1]->cd(); |
| 2853 | //canvasQA->cd(2); |
| 2854 | graphKaonTPCfinal->SetMarkerStyle(26); |
| 2855 | graphKaonTPCfinal->Draw("psame"); |
| 2856 | |
| 2857 | for(Int_t ip = 0; ip < graphPionTPCfinal->GetN(); ip ++) { |
| 2858 | graphPionTPCfinal->GetY()[ip] -= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]); |
| 2859 | graphPionTPCfinal->GetY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]); |
| 2860 | graphPionTPCfinal->GetEY()[ip] /= 50. * splinePion->Eval(graphPionTPCfinal->GetX()[ip]);; |
| 2861 | } |
| 2862 | canvasQA[2]->cd(); |
| 2863 | //canvasQA->cd(3); |
| 2864 | graphPionTPCfinal->SetMarkerStyle(26); |
| 2865 | graphPionTPCfinal->Draw("psame"); |
| 2866 | |
| 2867 | for(Int_t ip = 0; ip < graphElectronTPCfinal->GetN(); ip ++) { |
| 2868 | graphElectronTPCfinal->GetY()[ip] -= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]); |
| 2869 | graphElectronTPCfinal->GetY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]); |
| 2870 | graphElectronTPCfinal->GetEY()[ip] /= 50. * splineElectron->Eval(graphElectronTPCfinal->GetX()[ip]);; |
| 2871 | } |
| 2872 | canvasQA[3]->cd(); |
| 2873 | //canvasQA->cd(4); |
| 2874 | graphElectronTPCfinal->SetMarkerStyle(26); |
| 2875 | graphElectronTPCfinal->Draw("psame"); |
| 2876 | |
| 2877 | TFile* fSave = TFile::Open("splines_QA_ResidualPolynomials.root", "RECREATE"); |
| 2878 | fSave->cd(); |
| 2879 | for (Int_t i = 0; i < 4; i++) |
| 2880 | canvasQA[i]->Write(); |
| 2881 | |
| 2882 | fSave->Close(); |
| 2883 | //canvasQA->SaveAs("splines_QA_ResidualPolynomials.root"); |
| 2884 | |
| 2885 | |
| 2886 | delete funcBB; |
| 2887 | |
| 2888 | return arrResponse; |
| 2889 | } |
| 2890 | |
| 2891 | |
| 2892 | //________________________________________________________________________ |
| 2893 | TF1* AliTPCcalibResidualPID::FitBB(TObjArray* inputGraphs, Bool_t isMC, Bool_t isPPb, const Bool_t useV0s, |
| 2894 | const Double_t * initialParameters, AliTPCcalibResidualPID::FitType fitType) { |
| 2895 | // |
| 2896 | // Fit Bethe-Bloch parametrisation to data points (inputGraphs) |
| 2897 | // |
| 2898 | const Int_t nPar = 6; |
| 2899 | TGraphErrors * graphAll = (TGraphErrors *) inputGraphs->FindObject("beamDataPoints"); |
| 2900 | // |
| 2901 | Float_t from = 0.9; //TODO ADJUST -> Very important |
| 2902 | Float_t to = graphAll->GetXaxis()->GetXmax() * 2; |
| 2903 | |
| 2904 | |
| 2905 | |
| 2906 | TF1* funcBB = 0x0; |
| 2907 | |
| 2908 | Double_t parametersBBForward[nPar] = { 0, }; |
| 2909 | |
| 2910 | TVirtualFitter::SetMaxIterations(5e6); |
| 2911 | |
| 2912 | if (fitType == AliTPCcalibResidualPID::kSaturatedLund) { |
| 2913 | printf("Fit function: Saturated Lund\n"); |
| 2914 | |
| 2915 | funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar); |
| 2916 | //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values |
| 2917 | //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation |
| 2918 | Double_t parametersBB[nPar] = {61.0, 8.7, 1.86, 0.85, 113.4, -38}; // Yields reasonable results for data and MC ~ all periods |
| 2919 | |
| 2920 | if (isPPb && !isMC) { |
| 2921 | parametersBB[0] = 51.6; |
| 2922 | parametersBB[1] = 9.7; |
| 2923 | parametersBB[2] = 1.62; |
| 2924 | parametersBB[3] = 0.99; |
| 2925 | parametersBB[4] = 104.4; |
| 2926 | parametersBB[5] = -27.0; |
| 2927 | } |
| 2928 | if (isMC) { |
| 2929 | parametersBB[0] = 41.4; |
| 2930 | parametersBB[1] = 8.6; |
| 2931 | parametersBB[2] = 2.2; |
| 2932 | parametersBB[3] = 0.92; |
| 2933 | parametersBB[4] = 90.4; |
| 2934 | parametersBB[5] = -20.0; |
| 2935 | } |
| 2936 | funcBB->SetParameters(parametersBB); |
| 2937 | |
| 2938 | Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10}; |
| 2939 | funcBB->SetParErrors(parameterErrorsBB); |
| 2940 | |
| 2941 | for (Int_t i = 0; i < nPar; i++) |
| 2942 | parametersBBForward[i] = parametersBB[i]; |
| 2943 | } |
| 2944 | else if (fitType == AliTPCcalibResidualPID::kLund) { |
| 2945 | printf("Fit function: Lund\n"); |
| 2946 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); |
| 2947 | |
| 2948 | funcBB = new TF1("SaturatedLund", SaturatedLund, from, to, nPar); |
| 2949 | //Double_t parametersBB[nPar] = {34.0446, 8.42221, 4.16724, 1.29473, 80.6663, 0}; //Xianguos values |
| 2950 | //Double_t parametersBB[nPar] = {35.5, 8.7, 2.0, 1.09, 75.6, 0}; // No saturation |
| 2951 | Double_t parametersBB[nPar] = {35.0, 7., 1.86, 1., 75., 0}; // Yields reasonable results for data and MC ~ all periods |
| 2952 | |
| 2953 | funcBB->SetParameters(parametersBB); |
| 2954 | |
| 2955 | Double_t parameterErrorsBB[nPar] = {5., 0.5, 0.2, 0.05, 10, 10}; |
| 2956 | funcBB->SetParErrors(parameterErrorsBB); |
| 2957 | |
| 2958 | funcBB->FixParameter(5, 0.0); // No saturation |
| 2959 | |
| 2960 | for (Int_t i = 0; i < nPar; i++) |
| 2961 | parametersBBForward[i] = parametersBB[i]; |
| 2962 | } |
| 2963 | else if (fitType == kAlephWithAdditionalParam) { |
| 2964 | printf("Fit function: Aleph with additional parameter\n"); |
| 2965 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); |
| 2966 | |
| 2967 | // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability. |
| 2968 | // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus, |
| 2969 | // fixed to unity |
| 2970 | funcBB = new TF1("funcAleph", |
| 2971 | "[0]/TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])*([1]-[2]-[5]*TMath::Power(TMath::Sqrt(1. + x*x)/x , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to); |
| 2972 | //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); |
| 2973 | |
| 2974 | Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06}; |
| 2975 | //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1}; |
| 2976 | funcBB->SetParameters(parametersBB); |
| 2977 | |
| 2978 | for (Int_t i = 0; i < nPar; i++) |
| 2979 | parametersBBForward[i] = parametersBB[i]; |
| 2980 | } |
| 2981 | else if (fitType == AliTPCcalibResidualPID::kAleph) { |
| 2982 | printf("Fit function: Aleph\n"); |
| 2983 | printf("****WARNING: Fit settings not tuned for this fit function, only for saturated Lund!\n"); |
| 2984 | |
| 2985 | // NOTE: This form is equivalent to the original form, but with parameter [2] redefined for better numerical stability. |
| 2986 | // The additional parameter [5] has been introduced later and is unity originally. It seems not to be needed and is, thus, |
| 2987 | // fixed to unity |
| 2988 | funcBB = new TF1("funcAleph", |
| 2989 | "[0]/TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])*([1]-[2]-[5]*TMath::Power(x/TMath::Sqrt(1. + x*x) , [3])-TMath::Log(1. + TMath::Power(x, -[4])*TMath::Exp(-[2])))", from, to); |
| 2990 | //OLD"[0]*([1]*TMath::Power(TMath::Sqrt(1 + x*x)/x , [3]) - [5] - TMath::Power(TMath::Sqrt(1 + x*x)/x , [3])*TMath::Log(TMath::Exp([2]) + 1/TMath::Power(x, [4])))", from, to); |
| 2991 | //TEST Double_t parametersBB[nPar] = {1.2, 26.0, -30.0, -2.15, 5.6, 1}; |
| 2992 | Double_t parametersBB[nPar] = {1.25, 27.5, -29.0, 2.2, 5.2, 1}; |
| 2993 | |
| 2994 | // For [5] floating: Double_t parametersBB[nPar] = {2.42,15.2,-16,-2.24,2.8, 0.057}; |
| 2995 | //OLD with different sign for [3] and [5] not fix: Double_t parametersBB[nPar] = {2.6,14.3,-15,2.2,2.7, 0.06}; |
| 2996 | //OLD with different sign for [3]: Double_t parametersBB[nPar] = {0.0762*50.3,10.632,TMath::Log(1.34e-05),1.863,1.948, 1}; |
| 2997 | funcBB->SetParameters(parametersBB); |
| 2998 | |
| 2999 | // Fix parameter 5 to original value of unity |
| 3000 | funcBB->FixParameter(5, 1); |
| 3001 | |
| 3002 | for (Int_t i = 0; i < nPar; i++) |
| 3003 | parametersBBForward[i] = parametersBB[i]; |
| 3004 | } |
| 3005 | else { |
| 3006 | printf("Error: fitType not supported!\n"); |
| 3007 | return 0x0; |
| 3008 | } |
| 3009 | |
| 3010 | funcBB->SetLineColor(2); |
| 3011 | funcBB->SetLineWidth(1); |
| 3012 | |
| 3013 | // Override initial parameters, if user provides some |
| 3014 | if (initialParameters) { |
| 3015 | printf("Setting user initial parameters!\n"); |
| 3016 | funcBB->SetParameters(initialParameters); |
| 3017 | } |
| 3018 | |
| 3019 | // |
| 3020 | // |
| 3021 | // |
| 3022 | TGraphErrors * graphDelta = new TGraphErrors(*graphAll); |
| 3023 | graphDelta->SetName("graphDelta"); |
| 3024 | |
| 3025 | // In MC case: Remove errors from fit -> Some data points with extremely small errors otherwise completely dominate the fit, |
| 3026 | // but these points could still be wrong due to systematics (low p effects, e.g.) |
| 3027 | if (isMC) { |
| 3028 | for(Int_t ip = 0; ip < graphAll->GetN(); ip++) { |
| 3029 | graphAll->SetPointError(ip, 0, 0); |
| 3030 | } |
| 3031 | } |
| 3032 | else if (useV0s) { |
| 3033 | // Increase weight in plateau by reducing errors. Neccessary because errors are large there compared to other data points |
| 3034 | for(Int_t ip = 0; ip < graphAll->GetN(); ip++) { |
| 3035 | if (graphAll->GetX()[ip] >= 2500) { |
| 3036 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); |
| 3037 | } |
| 3038 | // Same for pions in the very rel. rise |
| 3039 | if (graphAll->GetX()[ip] >= 25 && graphAll->GetX()[ip] < 1000) { |
| 3040 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); |
| 3041 | } |
| 3042 | // Same for protons in the MIP region |
| 3043 | if (graphAll->GetX()[ip] >= 2 && graphAll->GetX()[ip] < 4) { |
| 3044 | graphAll->SetPointError(ip, 0, graphAll->GetEY()[ip] / 6.0); |
| 3045 | } |
| 3046 | } |
| 3047 | } |
| 3048 | graphAll->Fit(funcBB, "REX0M"); |
| 3049 | funcBB->SetRange(from, to); |
| 3050 | funcBB->GetParameters(parametersBBForward); |
| 3051 | // |
| 3052 | // |
| 3053 | // |
| 3054 | for(Int_t ip = 0; ip < graphDelta->GetN(); ip++) { |
| 3055 | graphDelta->GetY()[ip] -= funcBB->Eval(graphDelta->GetX()[ip]); |
| 3056 | graphDelta->GetY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]); |
| 3057 | graphDelta->GetEY()[ip] /= funcBB->Eval(graphDelta->GetX()[ip]); |
| 3058 | } |
| 3059 | TCanvas * canvDelta_1 = new TCanvas("canvDelta_1","control histogram for Bethe-Bloch fit 1",100,10,1380,800); |
| 3060 | TCanvas * canvDelta_2 = new TCanvas("canvDelta_2","control histogram for Bethe-Bloch fit 2",100,10,1380,800); |
| 3061 | canvDelta_1->SetGrid(1, 1); |
| 3062 | canvDelta_1->SetLogx(); |
| 3063 | canvDelta_2->SetGrid(1, 1); |
| 3064 | canvDelta_2->SetLogx(); |
| 3065 | /* |
| 3066 | canvDelta->Divide(2, 1); |
| 3067 | canvDelta->GetPad(1)->SetGrid(1, 1); |
| 3068 | canvDelta->GetPad(1)->SetLogx(); |
| 3069 | canvDelta->GetPad(2)->SetGrid(1, 1); |
| 3070 | canvDelta->GetPad(2)->SetLogx(); |
| 3071 | */ |
| 3072 | |
| 3073 | canvDelta_1->cd(); |
| 3074 | //canvDelta->cd(1); |
| 3075 | TH1F *hBBdummy=new TH1F("hBBdummy","BB fit;#beta#gamma;#LTdE/dx#GT (arb. unit)",100,.8,1e4); |
| 3076 | hBBdummy->SetMinimum(45); |
| 3077 | hBBdummy->SetMaximum(120); |
| 3078 | hBBdummy->GetXaxis()->SetTitleOffset(1.1); |
| 3079 | hBBdummy->SetStats(kFALSE); |
| 3080 | hBBdummy->Draw(); |
| 3081 | |
| 3082 | graphAll->SetTitle("BB multi-graph fit"); |
| 3083 | graphAll->SetMarkerStyle(22); |
| 3084 | graphAll->SetMarkerColor(kMagenta); |
| 3085 | graphAll->Draw("p"); |
| 3086 | |
| 3087 | TLegend *leg=new TLegend(.7,.12,.89,isMC?.4:.6); |
| 3088 | leg->SetFillColor(10); |
| 3089 | leg->SetBorderSize(1); |
| 3090 | |
| 3091 | //overlay individual graphs |
| 3092 | TGraph *gr=0x0; |
| 3093 | if (isMC) { |
| 3094 | gr=(TGraph*)inputGraphs->FindObject("Protons_MC_all"); |
| 3095 | gr->SetMarkerColor(kRed); |
| 3096 | gr->SetMarkerStyle(24); |
| 3097 | gr->Draw("p"); |
| 3098 | leg->AddEntry(gr,"p (MC)","p"); |
| 3099 | gr=(TGraph*)inputGraphs->FindObject("Pions_MC_all"); |
| 3100 | gr->SetMarkerColor(kBlue); |
| 3101 | gr->SetMarkerStyle(24); |
| 3102 | gr->Draw("p"); |
| 3103 | leg->AddEntry(gr,"#pi (MC)","p"); |
| 3104 | gr=(TGraph*)inputGraphs->FindObject("Kaons_MC_all"); |
| 3105 | gr->SetMarkerColor(kGreen); |
| 3106 | gr->SetMarkerStyle(24); |
| 3107 | gr->Draw("p"); |
| 3108 | leg->AddEntry(gr,"K (MC)","p"); |
| 3109 | gr=(TGraph*)inputGraphs->FindObject("Electrons_MC_all"); |
| 3110 | gr->SetMarkerColor(kBlack); |
| 3111 | gr->SetMarkerStyle(24); |
| 3112 | gr->Draw("p"); |
| 3113 | leg->AddEntry(gr,"e (MC)","p"); |
| 3114 | } |
| 3115 | else { |
| 3116 | gr=(TGraph*)inputGraphs->FindObject("protonTpcGraph"); |
| 3117 | gr->SetMarkerColor(kRed); |
| 3118 | gr->SetMarkerStyle(26); |
| 3119 | gr->Draw("p"); |
| 3120 | leg->AddEntry(gr,"p (TPC)","p"); |
| 3121 | gr=(TGraph*)inputGraphs->FindObject("protonTofGraph"); |
| 3122 | gr->SetMarkerColor(kRed); |
| 3123 | gr->SetMarkerStyle(25); |
| 3124 | gr->Draw("p"); |
| 3125 | leg->AddEntry(gr,"p (TOF)","p"); |
| 3126 | gr=(TGraph*)inputGraphs->FindObject("protonV0Graph");//ForBBfit"); |
| 3127 | gr->SetMarkerColor(kRed); |
| 3128 | gr->SetMarkerStyle(24); |
| 3129 | gr->Draw("p"); |
| 3130 | leg->AddEntry(gr,"p (V0)","p"); |
| 3131 | gr=(TGraph*)inputGraphs->FindObject("protonV0plusTOFGraph");//ForBBfit"); |
| 3132 | gr->SetMarkerColor(kRed); |
| 3133 | gr->SetMarkerStyle(30); |
| 3134 | gr->Draw("p"); |
| 3135 | leg->AddEntry(gr,"p (V0+TOF)","p"); |
| 3136 | gr=(TGraph*)inputGraphs->FindObject("pionTpcGraph"); |
| 3137 | gr->SetMarkerColor(kBlue); |
| 3138 | gr->SetMarkerStyle(26); |
| 3139 | gr->Draw("p"); |
| 3140 | leg->AddEntry(gr,"#pi (TPC)","p"); |
| 3141 | gr=(TGraph*)inputGraphs->FindObject("pionTofGraph"); |
| 3142 | gr->SetMarkerColor(kBlue); |
| 3143 | gr->SetMarkerStyle(25); |
| 3144 | gr->Draw("p"); |
| 3145 | leg->AddEntry(gr,"#pi (TOF)","p"); |
| 3146 | gr=(TGraph*)inputGraphs->FindObject("pionV0Graph");//ForBBfit"); |
| 3147 | gr->SetMarkerColor(kBlue); |
| 3148 | gr->SetMarkerStyle(24); |
| 3149 | gr->Draw("p"); |
| 3150 | leg->AddEntry(gr,"#pi (V0)","p"); |
| 3151 | gr=(TGraph*)inputGraphs->FindObject("pionV0plusTOFGraph");//ForBBfit"); |
| 3152 | gr->SetMarkerColor(kBlue); |
| 3153 | gr->SetMarkerStyle(30); |
| 3154 | gr->Draw("p"); |
| 3155 | leg->AddEntry(gr,"#pi (V0+TOF)","p"); |
| 3156 | gr=(TGraph*)inputGraphs->FindObject("kaonTpcGraph"); |
| 3157 | gr->SetMarkerColor(kGreen); |
| 3158 | gr->SetMarkerStyle(26); |
| 3159 | gr->Draw("p"); |
| 3160 | leg->AddEntry(gr,"K (TPC)","p"); |
| 3161 | gr=(TGraph*)inputGraphs->FindObject("kaonTofGraph"); |
| 3162 | gr->SetMarkerColor(kGreen); |
| 3163 | gr->SetMarkerStyle(25); |
| 3164 | gr->Draw("p"); |
| 3165 | leg->AddEntry(gr,"K (TOF)","p"); |
| 3166 | gr=(TGraph*)inputGraphs->FindObject("electronGraph"); |
| 3167 | gr->SetMarkerColor(kBlack); |
| 3168 | gr->SetMarkerStyle(25); |
| 3169 | gr->Draw("p"); |
| 3170 | leg->AddEntry(gr,"e","p"); |
| 3171 | gr=(TGraph*)inputGraphs->FindObject("electronV0Graph");//ForBBfit"); |
| 3172 | gr->SetMarkerColor(kBlack); |
| 3173 | gr->SetMarkerStyle(24); |
| 3174 | gr->Draw("p"); |
| 3175 | leg->AddEntry(gr,"e (V0)","p"); |
| 3176 | gr=(TGraph*)inputGraphs->FindObject("electronV0plusTOFGraph");//ForBBfit"); |
| 3177 | gr->SetMarkerColor(kBlack); |
| 3178 | gr->SetMarkerStyle(30); |
| 3179 | gr->Draw("p"); |
| 3180 | leg->AddEntry(gr,"e (V0+TOF)","p"); |
| 3181 | } |
| 3182 | |
| 3183 | graphAll->GetXaxis()->SetTitle("#beta#gamma"); |
| 3184 | graphAll->GetYaxis()->SetTitle("<dE/dx> (arb. unit)"); |
| 3185 | graphAll->Draw("p"); |
| 3186 | funcBB->GetHistogram()->DrawClone("csame"); |
| 3187 | |
| 3188 | leg->AddEntry(graphAll,"used","p"); |
| 3189 | leg->AddEntry(funcBB,"fit","l"); |
| 3190 | leg->Draw("same"); |
| 3191 | |
| 3192 | canvDelta_2->cd(); |
| 3193 | //canvDelta->cd(2); |
| 3194 | TH1F *hBBresdummy=new TH1F("hBBresdummy","residuals of BB fit;#beta#gamma;(data-fit)/fit",100,.8,1e4); |
| 3195 | hBBresdummy->SetMinimum(-0.04); |
| 3196 | hBBresdummy->SetMaximum(0.04); |
| 3197 | hBBresdummy->GetXaxis()->SetTitleOffset(1.1); |
| 3198 | hBBresdummy->GetYaxis()->SetTitleOffset(0.8); |
| 3199 | hBBresdummy->SetStats(kFALSE); |
| 3200 | hBBresdummy->Draw(); |
| 3201 | |
| 3202 | graphDelta->SetTitle("residuals of BB fit to multi-graph"); |
| 3203 | graphDelta->GetXaxis()->SetTitle("#beta#gamma"); |
| 3204 | graphDelta->GetYaxis()->SetTitle("(data - fit) / fit"); |
| 3205 | graphDelta->SetMarkerStyle(22); |
| 3206 | graphDelta->SetMarkerColor(4); |
| 3207 | graphDelta->Draw("p"); |
| 3208 | |
| 3209 | TFile* fSave = TFile::Open("splines_QA_BetheBlochFit.root", "RECREATE"); |
| 3210 | fSave->cd(); |
| 3211 | canvDelta_1->Write(); |
| 3212 | canvDelta_2->Write(); |
| 3213 | |
| 3214 | TString fitResults = "Fit results:\n"; |
| 3215 | for (Int_t i = 0; i < nPar; i++) { |
| 3216 | fitResults.Append(Form("par%d:\t%f +- %f\n", i, funcBB->GetParameters()[i], funcBB->GetParErrors()[i])); |
| 3217 | } |
| 3218 | |
| 3219 | TNamed* settings = new TNamed(fitResults.Data(), fitResults.Data()); |
| 3220 | settings->Write(); |
| 3221 | fSave->Close(); |
| 3222 | |
| 3223 | |
| 3224 | printf("\n\n%s", fitResults.Data()); |
| 3225 | |
| 3226 | return funcBB; |
| 3227 | } |
| 3228 | |
| 3229 | |
| 3230 | //________________________________________________________________________ |
| 3231 | Double_t AliTPCcalibResidualPID::Lund(Double_t* xx, Double_t* par) |
| 3232 | { |
| 3233 | // bg is beta*gamma |
| 3234 | const Double_t bg = xx[0]; |
| 3235 | |
| 3236 | const Double_t beta2 = bg*bg / (1.0 + bg*bg); |
| 3237 | |
| 3238 | const Double_t a = par[0]; |
| 3239 | const Double_t b = par[1]; |
| 3240 | const Double_t c = par[2]; |
| 3241 | const Double_t e = par[3]; |
| 3242 | const Double_t f = par[4]; |
| 3243 | |
| 3244 | const Double_t d = TMath::Exp(c*(a - f) / b); |
| 3245 | |
| 3246 | const Double_t powbg = TMath::Power(1.0 + bg, c); |
| 3247 | |
| 3248 | const Double_t value = a / TMath::Power(beta2,e) + |
| 3249 | b/c * TMath::Log(powbg / (1.0 + d*powbg)); |
| 3250 | |
| 3251 | return value; |
| 3252 | } |
| 3253 | |
| 3254 | |
| 3255 | //________________________________________________________________________ |
| 3256 | Double_t AliTPCcalibResidualPID::SaturatedLund(Double_t* xx, Double_t* par) |
| 3257 | { |
| 3258 | const Double_t qq = Lund(xx, par); |
| 3259 | return qq * TMath::Exp(par[5] / qq); |
| 3260 | } |
| 3261 | |
| 3262 | |
| 3263 | //________________________________________________________________________ |
| 3264 | void AliTPCcalibResidualPID::FitSlicesY(TH2 *hist, Double_t heightFractionForRange, Int_t cutThreshold, TString fitOption, TObjArray *arr) |
| 3265 | { |
| 3266 | //heightPercentageForRange |
| 3267 | // custom slices fit |
| 3268 | // |
| 3269 | |
| 3270 | if (!hist) return; |
| 3271 | if (!arr) return; |
| 3272 | |
| 3273 | // If old style is to be used |
| 3274 | /* |
| 3275 | hist->FitSlicesY(0, 0, -1, cutThreshold, fitOption.Data(), &arr); |
| 3276 | return; |
| 3277 | */ |
| 3278 | |
| 3279 | |
| 3280 | arr->Clear(); |
| 3281 | arr->SetOwner(); |
| 3282 | arr->Expand(4); |
| 3283 | |
| 3284 | TAxis *axis=hist->GetXaxis(); |
| 3285 | const TArrayD *bins = axis->GetXbins(); |
| 3286 | TH1D** hList = new TH1D*[4]; |
| 3287 | |
| 3288 | for (Int_t i = 0; i < 4; i++) { |
| 3289 | delete gDirectory->FindObject(Form("%s_%d", hist->GetName(), i)); |
| 3290 | |
| 3291 | if (bins->fN == 0) { |
| 3292 | hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", |
| 3293 | hist->GetNbinsX(), axis->GetXmin(), axis->GetXmax()); |
| 3294 | } else { |
| 3295 | hList[i] = new TH1D(Form("%s_%d", hist->GetName(), i), i < 3 ? Form("Fitted value of par[%d]", i) : "Chi2/NDF", hist->GetNbinsX(), bins->fArray); |
| 3296 | } |
| 3297 | |
| 3298 | (*arr)[i] = hList[i]; |
| 3299 | } |
| 3300 | |
| 3301 | for (Int_t ibin=axis->GetFirst(); ibin<=axis->GetLast(); ++ibin){ |
| 3302 | TH1 *h=hist->ProjectionY("_temp",ibin,ibin); |
| 3303 | if (!h) |
| 3304 | continue; |
| 3305 | |
| 3306 | if (h->GetEntries() < cutThreshold) { |
| 3307 | delete h; |
| 3308 | continue; |
| 3309 | } |
| 3310 | |
| 3311 | // Average around maximum bin -> Might catch some outliers |
| 3312 | Int_t maxBin = h->GetMaximumBin(); |
| 3313 | Double_t maxVal = h->GetBinContent(maxBin); |
| 3314 | |
| 3315 | if (maxVal < 2) { // It could happen that all entries are in overflow/underflow; don't fit in this case |
| 3316 | delete h; |
| 3317 | continue; |
| 3318 | } |
| 3319 | |
| 3320 | UChar_t usedBins = 1; |
| 3321 | if (maxBin > 1) { |
| 3322 | maxVal += h->GetBinContent(maxBin - 1); |
| 3323 | usedBins++; |
| 3324 | } |
| 3325 | if (maxBin < h->GetNbinsX()) { |
| 3326 | maxVal += h->GetBinContent(maxBin + 1); |
| 3327 | usedBins++; |
| 3328 | } |
| 3329 | maxVal /= usedBins; |
| 3330 | |
| 3331 | Double_t thresholdFraction = heightFractionForRange * maxVal; |
| 3332 | Int_t lowThrBin = TMath::Max(1, h->FindFirstBinAbove(thresholdFraction)); |
| 3333 | Int_t highThrBin = TMath::Min(h->GetNbinsX(), h->FindLastBinAbove(thresholdFraction)); |
| 3334 | |
| 3335 | Double_t lowThreshold = h->GetBinCenter(lowThrBin); |
| 3336 | Double_t highThreshold = h->GetBinCenter(highThrBin); |
| 3337 | |
| 3338 | TFitResultPtr res = h->Fit("gaus", Form("%sS", fitOption.Data()), "", lowThreshold, highThreshold); |
| 3339 | |
| 3340 | if ((Int_t)res == 0) { |
| 3341 | Int_t resBin = ibin; |
| 3342 | hList[0]->SetBinContent(resBin,res->GetParams()[0]); |
| 3343 | hList[0]->SetBinError(resBin,res->GetErrors()[0]); |
| 3344 | hList[1]->SetBinContent(resBin,res->GetParams()[1]); |
| 3345 | hList[1]->SetBinError(resBin,res->GetErrors()[1]); |
| 3346 | hList[2]->SetBinContent(resBin,res->GetParams()[2]); |
| 3347 | hList[2]->SetBinError(resBin,res->GetErrors()[2]); |
| 3348 | hList[3]->SetBinContent(resBin,res->Ndf()>0?res->Chi2()/res->Ndf():0); |
| 3349 | } |
| 3350 | |
| 3351 | delete h; |
| 3352 | } |
| 3353 | } |
| 3354 | |
| 3355 | |
| 3356 | //______________________________________________________________________________ |
| 3357 | Bool_t AliTPCcalibResidualPID::GetVertexIsOk(AliVEvent* event) const |
| 3358 | { |
| 3359 | AliAODEvent* aod = 0x0; |
| 3360 | AliESDEvent* esd = 0x0; |
| 3361 | |
| 3362 | aod = dynamic_cast<AliAODEvent*>(event); |
| 3363 | if (!aod) { |
| 3364 | esd = dynamic_cast<AliESDEvent*>(event); |
| 3365 | |
| 3366 | if (!esd) { |
| 3367 | AliError("Event seems to be neither AOD nor ESD!"); |
| 3368 | return kFALSE; |
| 3369 | } |
| 3370 | } |
| 3371 | |
| 3372 | if (fIsPbpOrpPb) { |
| 3373 | const AliVVertex* trkVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : |
| 3374 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex())); |
| 3375 | |
| 3376 | if (!trkVtx || trkVtx->GetNContributors() <= 0) |
| 3377 | return kFALSE; |
| 3378 | |
| 3379 | TString vtxTtl = trkVtx->GetTitle(); |
| 3380 | if (!vtxTtl.Contains("VertexerTracks")) |
| 3381 | return kFALSE; |
| 3382 | |
| 3383 | Float_t zvtx = trkVtx->GetZ(); |
| 3384 | const AliVVertex* spdVtx = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) : |
| 3385 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD())); |
| 3386 | if (spdVtx->GetNContributors() <= 0) |
| 3387 | return kFALSE; |
| 3388 | |
| 3389 | TString vtxTyp = spdVtx->GetTitle(); |
| 3390 | Double_t cov[6] = {0}; |
| 3391 | spdVtx->GetCovarianceMatrix(cov); |
| 3392 | Double_t zRes = TMath::Sqrt(cov[5]); |
| 3393 | if (vtxTyp.Contains("vertexer:Z") && (zRes > 0.25)) |
| 3394 | return kFALSE; |
| 3395 | |
| 3396 | if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ()) > 0.5) |
| 3397 | return kFALSE; |
| 3398 | |
| 3399 | if (TMath::Abs(zvtx) > fZvtxCutEvent) //Default: 10 cm |
| 3400 | return kFALSE; |
| 3401 | |
| 3402 | return kTRUE; |
| 3403 | } |
| 3404 | |
| 3405 | |
| 3406 | // pp and PbPb |
| 3407 | const AliVVertex* primaryVertex = (aod ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) : |
| 3408 | dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexTracks())); |
| 3409 | |
| 3410 | if (!primaryVertex || primaryVertex->GetNContributors() <= 0) |
| 3411 | return kFALSE; |
| 3412 | |
| 3413 | if (TMath::Abs(primaryVertex->GetZ()) >= fZvtxCutEvent) //Default: 10 cm |
| 3414 | return kFALSE; |
| 3415 | |
| 3416 | return kTRUE; |
| 3417 | } |
| 3418 | |
| 3419 | |
| 3420 | //______________________________________________________________________________ |
| 3421 | void AliTPCcalibResidualPID::FillV0PIDlist(AliESDEvent* event) |
| 3422 | { |
| 3423 | // |
| 3424 | // Fill the PID tag list |
| 3425 | // |
| 3426 | |
| 3427 | // If no event forwarded as parameter (default), cast current input event. |
| 3428 | // Dynamic cast to ESD events (DO NOTHING for AOD events) |
| 3429 | if (!event) |
| 3430 | event = dynamic_cast<AliESDEvent *>(InputEvent()); |
| 3431 | |
| 3432 | // If this fails, do nothing |
| 3433 | if (!event) |
| 3434 | return; |
| 3435 | |
| 3436 | if (!fV0KineCuts) { |
| 3437 | AliError("V0KineCuts not available!"); |
| 3438 | return; |
| 3439 | } |
| 3440 | |
| 3441 | TString beamType(event->GetBeamType()); |
| 3442 | |
| 3443 | if (beamType.CompareTo("Pb-Pb") == 0 || beamType.CompareTo("A-A") == 0) { |
| 3444 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPbPb); |
| 3445 | } |
| 3446 | else { |
| 3447 | fV0KineCuts->SetMode(AliESDv0KineCuts::kPurity, AliESDv0KineCuts::kPP); |
| 3448 | } |
| 3449 | |
| 3450 | // V0 selection |
| 3451 | // set event |
| 3452 | fV0KineCuts->SetEvent(event); |
| 3453 | |
| 3454 | const Int_t numTracks = event->GetNumberOfTracks(); |
| 3455 | fV0tags = new Char_t[numTracks]; |
| 3456 | for (Int_t i = 0; i < numTracks; i++) |
| 3457 | fV0tags[i] = 0; |
| 3458 | |
| 3459 | fV0motherIndex = new Int_t[numTracks]; |
| 3460 | for (Int_t i = 0; i < numTracks; i++) |
| 3461 | fV0motherIndex[i] = -1; |
| 3462 | |
| 3463 | fV0motherPDG = new Int_t[numTracks]; |
| 3464 | for (Int_t i = 0; i < numTracks; i++) |
| 3465 | fV0motherPDG[i] = 0; |
| 3466 | |
| 3467 | fNumTagsStored = numTracks; |
| 3468 | |
| 3469 | |
| 3470 | |
| 3471 | |
| 3472 | |
| 3473 | // loop over V0 particles |
| 3474 | for (Int_t iv0 = 0; iv0 < event->GetNumberOfV0s(); iv0++) { |
| 3475 | AliESDv0* v0 = (AliESDv0*)event->GetV0(iv0); |
| 3476 | |
| 3477 | if (!v0) |
| 3478 | continue; |
| 3479 | |
| 3480 | // Reject onFly V0's <-> Only take V0's from offline V0 finder |
| 3481 | if (v0->GetOnFlyStatus()) |
| 3482 | continue; |
| 3483 | |
| 3484 | // Get the particle selection |
| 3485 | Bool_t foundV0 = kFALSE; |
| 3486 | Int_t pdgV0 = 0, pdgP = 0, pdgN = 0; |
| 3487 | |
| 3488 | foundV0 = fV0KineCuts->ProcessV0(v0, pdgV0, pdgP, pdgN); |
| 3489 | if (!foundV0) |
| 3490 | continue; |
| 3491 | |
| 3492 | Int_t iTrackP = v0->GetPindex(); // positive track |
| 3493 | Int_t iTrackN = v0->GetNindex(); // negative track |
| 3494 | |
| 3495 | /* |
| 3496 | AliESDtrack* trackP = event->GetTrack(iTrackP); |
| 3497 | AliESDtrack* trackN = event->GetTrack(iTrackN); |
| 3498 | |
| 3499 | if (!trackP || !trackN) |
| 3500 | continue; |
| 3501 | |
| 3502 | |
| 3503 | |
| 3504 | Float_t xy = 999, z = 999; |
| 3505 | trackP->GetImpactParameters(xy, z); |
| 3506 | |
| 3507 | Bool_t reject = kFALSE; |
| 3508 | if (TMath::Abs(z) > 1) |
| 3509 | continue; |
| 3510 | |
| 3511 | trackN->GetImpactParameters(xy, z); |
| 3512 | if (TMath::Abs(z) > 1) |
| 3513 | continue; |
| 3514 | */ |
| 3515 | |
| 3516 | |
| 3517 | // Fill the Object arrays |
| 3518 | // positive particles |
| 3519 | if (pdgP == -11) { |
| 3520 | fV0tags[iTrackP] = 14; |
| 3521 | } |
| 3522 | else if (pdgP == 211) { |
| 3523 | fV0tags[iTrackP] = 15; |
| 3524 | } |
| 3525 | else if(pdgP == 2212) { |
| 3526 | fV0tags[iTrackP] = 16; |
| 3527 | } |
| 3528 | |
| 3529 | fV0motherIndex[iTrackP] = iv0; |
| 3530 | fV0motherPDG[iTrackP] = pdgV0; |
| 3531 | |
| 3532 | // negative particles |
| 3533 | if( pdgN == 11){ |
| 3534 | fV0tags[iTrackN] = -14; |
| 3535 | } |
| 3536 | else if( pdgN == -211){ |
| 3537 | fV0tags[iTrackN] = -15; |
| 3538 | } |
| 3539 | else if( pdgN == -2212){ |
| 3540 | fV0tags[iTrackN] = -16; |
| 3541 | } |
| 3542 | |
| 3543 | fV0motherIndex[iTrackN] = iv0; |
| 3544 | fV0motherPDG[iTrackN] = pdgV0; |
| 3545 | |
| 3546 | } |
| 3547 | } |
| 3548 | |
| 3549 | |
| 3550 | //______________________________________________________________________________ |
| 3551 | void AliTPCcalibResidualPID::ClearV0PIDlist() |
| 3552 | { |
| 3553 | // |
| 3554 | // Clear the PID tag list |
| 3555 | // |
| 3556 | |
| 3557 | delete fV0tags; |
| 3558 | fV0tags = 0; |
| 3559 | |
| 3560 | delete fV0motherIndex; |
| 3561 | fV0motherIndex = 0; |
| 3562 | |
| 3563 | delete fV0motherPDG; |
| 3564 | fV0motherPDG = 0; |
| 3565 | |
| 3566 | fNumTagsStored = 0; |
| 3567 | } |
| 3568 | |
| 3569 | |
| 3570 | //______________________________________________________________________________ |
| 3571 | Char_t AliTPCcalibResidualPID::GetV0tag(Int_t trackIndex) const |
| 3572 | { |
| 3573 | // |
| 3574 | // Get the tag for the corresponding trackIndex. Returns -99 in case of invalid index/tag list. |
| 3575 | // |
| 3576 | |
| 3577 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0tags) |
| 3578 | return -99; |
| 3579 | |
| 3580 | return fV0tags[trackIndex]; |
| 3581 | } |
| 3582 | |
| 3583 | |
| 3584 | //______________________________________________________________________________ |
| 3585 | Int_t AliTPCcalibResidualPID::GetV0motherIndex(Int_t trackIndex) const |
| 3586 | { |
| 3587 | // |
| 3588 | // Get the index of the V0 mother for the corresponding trackIndex. Returns -99 in case of invalid index/mother index list. |
| 3589 | // |
| 3590 | |
| 3591 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherIndex) |
| 3592 | return -99; |
| 3593 | |
| 3594 | return fV0motherIndex[trackIndex]; |
| 3595 | } |
| 3596 | |
| 3597 | |
| 3598 | //______________________________________________________________________________ |
| 3599 | Int_t AliTPCcalibResidualPID::GetV0motherPDG(Int_t trackIndex) const |
| 3600 | { |
| 3601 | // |
| 3602 | // Get the PDG code of the V0 mother for the corresponding trackIndex. Returns 0 in case of invalid index/mother index list. |
| 3603 | // |
| 3604 | |
| 3605 | if (trackIndex < 0 || trackIndex >= fNumTagsStored || !fV0motherPDG) |
| 3606 | return 0; |
| 3607 | |
| 3608 | return fV0motherPDG[trackIndex]; |
| 3609 | } |
| 3610 | |
| 3611 | |
| 3612 | //______________________________________________________________________________ |
| 3613 | Int_t AliTPCcalibResidualPID::MergeGraphErrors(TGraphErrors* mergedGraph, TCollection* li) |
| 3614 | { |
| 3615 | // Adds all graphs with errors from the collection to this graph. |
| 3616 | // Returns the total number of poins in the result or -1 in case of an error. |
| 3617 | |
| 3618 | // TODO "Normal" merge of latest root trunk will also take into account the error bars, |
| 3619 | // so this function can be replaced by the normal root function with the latest root version. |
| 3620 | |
| 3621 | TIter next(li); |
| 3622 | while (TObject* o = next()) { |
| 3623 | TGraph *g = dynamic_cast<TGraph*>(o); |
| 3624 | if (!g) { |
| 3625 | Printf("ERROR: Cannot merge an object which doesn't inherit from TGraph found in the list"); |
| 3626 | return -1; |
| 3627 | } |
| 3628 | int n0 = mergedGraph->GetN(); |
| 3629 | int n1 = n0+g->GetN(); |
| 3630 | mergedGraph->Set(n1); |
| 3631 | Double_t * x = g->GetX(); |
| 3632 | Double_t * y = g->GetY(); |
| 3633 | Double_t * ex = g->GetEX(); |
| 3634 | Double_t * ey = g->GetEY(); |
| 3635 | for (Int_t i = 0 ; i < g->GetN(); i++) { |
| 3636 | mergedGraph->SetPoint(n0+i, x[i], y[i]); |
| 3637 | Double_t exPoint = ex ? ex[i] : 0; |
| 3638 | Double_t eyPoint = ey ? ey[i] : 0; |
| 3639 | mergedGraph->SetPointError(n0+i, exPoint, eyPoint); |
| 3640 | } |
| 3641 | } |
| 3642 | return mergedGraph->GetN(); |
| 3643 | } |
| 3644 | |
| 3645 | //________________________________________________________________________ |
| 3646 | Bool_t AliTPCcalibResidualPID::TPCCutMIGeo(const AliVTrack* track, const AliVEvent* evt, TTreeStream* streamer) |
| 3647 | { |
| 3648 | // |
| 3649 | // TPC Cut MIGeo |
| 3650 | // |
| 3651 | |
| 3652 | if (!track || !evt) |
| 3653 | return kFALSE; |
| 3654 | |
| 3655 | const Short_t sign = track->Charge(); |
| 3656 | Double_t xyz[50]; |
| 3657 | Double_t pxpypz[50]; |
| 3658 | Double_t cv[100]; |
| 3659 | |
| 3660 | track->GetXYZ(xyz); |
| 3661 | track->GetPxPyPz(pxpypz); |
| 3662 | |
| 3663 | AliExternalTrackParam* par = new AliExternalTrackParam(xyz, pxpypz, cv, sign); |
| 3664 | const AliESDtrack dummy; |
| 3665 | |
| 3666 | const Double_t magField = evt->GetMagneticField(); |
| 3667 | Double_t varGeom = dummy.GetLengthInActiveZone(par, 3, 236, magField, 0, 0); |
| 3668 | Double_t varNcr = track->GetTPCClusterInfo(3, 1); |
| 3669 | Double_t varNcls = track->GetTPCsignalN(); |
| 3670 | |
| 3671 | const Double_t varEval = 130. - 5. * TMath::Abs(1. / track->Pt()); |
| 3672 | Bool_t cutGeom = varGeom > fgCutGeo * varEval; |
| 3673 | Bool_t cutNcr = varNcr > fgCutNcr * varEval; |
| 3674 | Bool_t cutNcls = varNcls > fgCutNcl * varEval; |
| 3675 | |
| 3676 | Bool_t kout = cutGeom && cutNcr && cutNcls; |
| 3677 | |
| 3678 | if (streamer) { |
| 3679 | Double_t dedx = track->GetTPCsignal(); |
| 3680 | |
| 3681 | (*streamer)<<"tree"<< |
| 3682 | "param.="<< par<< |
| 3683 | "varGeom="<<varGeom<< |
| 3684 | "varNcr="<<varNcr<< |
| 3685 | "varNcls="<<varNcls<< |
| 3686 | |
| 3687 | "cutGeom="<<cutGeom<< |
| 3688 | "cutNcr="<<cutNcr<< |
| 3689 | "cutNcls="<<cutNcls<< |
| 3690 | |
| 3691 | "kout="<<kout<< |
| 3692 | "dedx="<<dedx<< |
| 3693 | |
| 3694 | "\n"; |
| 3695 | } |
| 3696 | |
| 3697 | delete par; |
| 3698 | |
| 3699 | return kout; |
| 3700 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame