]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Pass0/merge.C CPass1/merge.C
[u/mrichter/AliRoot.git] / PWGPP / EVCHAR / AliAnalysisTaskSPDdNdEta.cxx
CommitLineData
f1fc7d61 1/*************************************************************************
2* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
e1366e5e 16//////////////////////////////////////////////////////////////////////////
17// Class AliAnalysisTaskSPDdNdEta //
18// Analysis task for dN/dEta reconstruction with the SPD //
19// //
20// Author: M. Nicassio (INFN Bari) //
21// Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it //
22//////////////////////////////////////////////////////////////////////////
f1fc7d61 23
ea23c361 24#include "TChain.h"
25#include "TTree.h"
ea23c361 26
27#include "TH1F.h"
28#include "TH2F.h"
29#include "TH3F.h"
39f606ae 30#include "TList.h"
86841765 31#include "TGeoGlobalMagField.h"
32#include "AliMagF.h"
f1fc7d61 33#include "AliAnalysisManager.h"
ea23c361 34
f1fc7d61 35#include "AliMultiplicity.h"
39f606ae 36#include "AliESDEvent.h"
f1fc7d61 37#include "AliESDInputHandler.h"
ea23c361 38
39f606ae 39#include "AliESDInputHandlerRP.h"
ccfafc80 40#include "AliESDVZERO.h"
39f606ae 41#include "../ITS/AliITSRecPoint.h"
42#include "AliCDBPath.h"
43#include "AliCDBManager.h"
44#include "AliCDBEntry.h"
45#include "AliCDBStorage.h"
46#include "AliGeomManager.h"
47
ea23c361 48#include "AliMCEventHandler.h"
49#include "AliMCEvent.h"
f1fc7d61 50#include "AliStack.h"
ea23c361 51#include "AliTrackReference.h"
52
39f606ae 53#include "AliGenHijingEventHeader.h"
ccfafc80 54#include "AliGenDPMjetEventHeader.h"
39f606ae 55
56#include "AliLog.h"
57
58#include "AliTriggerAnalysis.h"
59#include "AliPhysicsSelection.h"
39f606ae 60#include "AliTrackletAlg.h"
ea23c361 61#include "AliAnalysisTaskSPDdNdEta.h"
62
63
64ClassImp(AliAnalysisTaskSPDdNdEta)
ea23c361 65//________________________________________________________________________
66AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
39f606ae 67 : AliAnalysisTaskSE(name),
ea23c361 68
39f606ae 69 fmyESD(0),
ea23c361 70 fOutput(0),
4ab6aaf0 71
39f606ae 72 fUseMC(kFALSE),
39f606ae 73 fTrigger(AliTriggerAnalysis::kAcceptAll),
74 fTR(kFALSE),
75 fRecoTracklets(kFALSE),
ea23c361 76
86841765 77 fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall),
78 fCentrLowLim(0),
79 fCentrUpLim(0),
ccfafc80 80 fCentrEst(kFALSE),
86841765 81 fMinClMultLay2(0),
ccfafc80 82 fMaxClMultLay2(0),
83 fMinMultV0(0),
84 fVtxLim(0),
85 fPartSpecies(0),
86841765 86
87 fPhiWindow(0),
88 fThetaWindow(0),
89 fPhiShift(0),
90 fRemoveClustersFromOverlaps(0),
91 fPhiOverlapCut(0),
92 fZetaOverlapCut(0),
93 fPhiRotationAngle(0),
94 fPhiWindowAna(0),
95
96 fMultReco(0),
97
ccfafc80 98 fV0Ampl(0),
ea23c361 99 fHistSPDRAWMultvsZ(0),
39f606ae 100 fHistSPDRAWMultvsZTriggCentrEvts(0),
101 fHistSPDRAWMultvsZCentrEvts(0),
ea23c361 102 fHistSPDRAWEtavsZ(0),
39f606ae 103
ea23c361 104 fHistSPDmultEtacut(0),
105 fHistSPDmult(0),
39f606ae 106 fHistSPDmultcl1(0),
86841765 107 fHistSPDmultcl2(0),
108 fHistSPDmultcl1vscl2(0),
109 fHistSPDmultvscl1(0),
110 fHistSPDmultvscl2(0),
111
39f606ae 112 fHistSPDeta(0),
ea23c361 113 fHistSPDphi(0),
ea23c361 114 fHistSPDtheta(0),
ea23c361 115 fHistSPDdePhi(0),
39f606ae 116 fHistSPDphivsSPDeta(0),
4ab6aaf0 117 fHistSPDdeTheta(0),
86841765 118 fHistSPDvtx(0),
39f606ae 119 fHistSPDvtxAnalysis(0),
120 fHistSPDdePhideTheta(0),
ea23c361 121
86841765 122 fHistSPDphicl1(0),
123 fHistSPDphicl2(0),
124
ea23c361 125 fHistBkgCorrDen(0),
ccfafc80 126 fHistReconstructedProtons(0),
127 fHistReconstructedKaons(0),
128 fHistReconstructedPions(0),
129 fHistReconstructedOthers(0),
130 fHistReconstructedSec(0),
39f606ae 131 fHistBkgCorrDenPrimGen(0),
ccfafc80 132 fHistBkgCombLabels(0),
39f606ae 133 fHistBkgCorrNum(0),
ea23c361 134 fHistAlgEffNum(0),
ea23c361 135 fHistNonDetectableCorrDen(0),
136
39f606ae 137 fHistNonDetectableCorrNum(0),
ccfafc80 138 fHistProtons(0),
139 fHistKaons(0),
140 fHistPions(0),
141 fHistOthers(0),
39f606ae 142 fHistAllPrimaries(0),
143 fHistTrackCentrEvts(0),
144 fHistTrackTrigCentrEvts(0),
ea23c361 145
39f606ae 146 fHistAllEvts(0),
147 fHistCentrEvts(0),
148 fHistTrigCentrEvts(0),
149 fHistSelEvts(0),
4ab6aaf0 150
39f606ae 151 fHistMCmultEtacut(0),
152 fHistMCmultEtacutvsSPDmultEtacut(0),
86841765 153 fHistMCmultEtacutvsSPDmultcl1(0),
154 fHistMCmultEtacutvsSPDmultcl2(0),
39f606ae 155
156 fHistMCvtxx(0),
157 fHistMCvtxy(0),
158 fHistMCvtxz(0),
159
160 fHistRecvsGenImpactPar(0),
161 fHistMCNpart(0),
162
ccfafc80 163 fHistdPhidThetaPP(0),
164 fHistdPhidThetaSS(0),
165 fHistdPhidThetaComb(0),
39f606ae 166
167 fHistDeVtx(0)
4ab6aaf0 168
ea23c361 169{
170
171 // Constructor
172
173 // Define input and output slots here
39f606ae 174 // Input slot #0 works with a TChain
175// DefineInput(0, TChain::Class());
176 // Output slot #0 id reserved by the base class for AOD
177 DefineOutput(1, TList::Class());
ea23c361 178
179}
180//________________________________________________________________________
181AliAnalysisTaskSPDdNdEta::~AliAnalysisTaskSPDdNdEta()
182{
183
184 // Destructor
185
186 // histograms are in the output list and deleted when the output
187 // list is deleted by the TSelector dtor
86841765 188 delete fMultReco;
ea23c361 189
39f606ae 190 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
ea23c361 191 delete fOutput;
192 fOutput = 0;
193 }
194
195}
196//________________________________________________________________________
39f606ae 197void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects()
ea23c361 198{
39f606ae 199 if (fRecoTracklets) {
200 AliCDBManager *man = AliCDBManager::Instance();
86841765 201 man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
202 if (fUseMC) man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
203 else man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
ccfafc80 204 man->SetRun(137161);
39f606ae 205 AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
e37e09b4 206 if(!obj) AliFatal("Unable to load geometry from CDB!");
39f606ae 207 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
208 AliGeomManager::GetNalignable("ITS");
209 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
ea23c361 210 }
86841765 211
212 fMultReco = new AliTrackletAlg();
213
ea23c361 214 // Create histograms
39f606ae 215 fOutput = new TList();
216 fOutput->SetOwner();
ea23c361 217
ccfafc80 218 Int_t nBinVtx = 40;
86841765 219 Double_t MaxVtx = 20.;
220
39f606ae 221 Int_t nBinMultCorr = 200;
222 Float_t nMaxMult = 20000.;
4ab6aaf0 223 Double_t binLimMultCorr[nBinMultCorr+1];
39f606ae 224 binLimMultCorr[0] = 0.;
225 binLimMultCorr[1] = 1.;
226 for (Int_t i = 2; i<=nBinMultCorr;++i) {
227 binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
4ab6aaf0 228 }
f1fc7d61 229
39f606ae 230 Int_t nBinEtaCorr = 60;
231 Float_t etaMax = 3.;
232 Float_t etaMin = -3.;
233 Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
234 for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
235 binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
236 }
ea23c361 237
39f606ae 238 // Data to be corrected
239 // ...event level
ccfafc80 240
241 fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
242 fOutput->Add(fV0Ampl);
243
86841765 244 fHistSPDRAWMultvsZ = new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 245 fHistSPDRAWMultvsZ->GetXaxis()->SetTitle("Tracklet multiplicity");
246 fHistSPDRAWMultvsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
39f606ae 247 fOutput->Add(fHistSPDRAWMultvsZ);
248
86841765 249 fHistSPDRAWMultvsZTriggCentrEvts = new TH2F("fHistSPDRAWMultvsZTriggCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 250 fHistSPDRAWMultvsZTriggCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
251 fHistSPDRAWMultvsZTriggCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
39f606ae 252 fOutput->Add(fHistSPDRAWMultvsZTriggCentrEvts);
253
86841765 254 fHistSPDRAWMultvsZCentrEvts = new TH2F("fHistSPDRAWMultvsZCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 255 fHistSPDRAWMultvsZCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
256 fHistSPDRAWMultvsZCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
39f606ae 257 fOutput->Add(fHistSPDRAWMultvsZCentrEvts);
258
259 // ...track level
86841765 260 fHistSPDRAWEtavsZ = new TH2F("fHistSPDRAWEtavsZ", "Tracklet pseudorapidity distribution", nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
ea23c361 261 fHistSPDRAWEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
262 fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
ea23c361 263 fOutput->Add(fHistSPDRAWEtavsZ);
ea23c361 264
39f606ae 265 fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution",nBinMultCorr,binLimMultCorr);
266 fHistSPDmultEtacut->GetXaxis()->SetTitle("Tracklet multiplicity (|#eta|<1.4)");
ea23c361 267 fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
268 fOutput->Add(fHistSPDmultEtacut);
269
39f606ae 270 fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", nBinMultCorr,binLimMultCorr);
271 fHistSPDmult->GetXaxis()->SetTitle("Tracklet multiplicity");
ea23c361 272 fHistSPDmult->GetYaxis()->SetTitle("Entries");
273 fOutput->Add(fHistSPDmult);
274
86841765 275 fHistSPDmultcl1 = new TH1F("fHistSPDmultcl1", "Inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
39f606ae 276 fHistSPDmultcl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
277 fHistSPDmultcl1->GetYaxis()->SetTitle("Entries");
278 fOutput->Add(fHistSPDmultcl1);
ea23c361 279
86841765 280 fHistSPDmultcl2 = new TH1F("fHistSPDmultcl2", "Outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
281 fHistSPDmultcl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
282 fHistSPDmultcl2->GetYaxis()->SetTitle("Entries");
283 fOutput->Add(fHistSPDmultcl2);
284
285 fHistSPDmultcl1vscl2 = new TH2F("fHistSPDmultcl1vscl2", "Inner layer cluster vs outer multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
286 fHistSPDmultcl1vscl2->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
287 fHistSPDmultcl1vscl2->GetYaxis()->SetTitle("Outer layer cluster multiplicity");
288 fOutput->Add(fHistSPDmultcl1vscl2);
289
290 fHistSPDmultvscl1 = new TH2F("fHistSPDmultvscl1", "Tracklet vs inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
291 fHistSPDmultvscl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
292 fHistSPDmultvscl1->GetYaxis()->SetTitle("Tracklet multiplicity");
293 fOutput->Add(fHistSPDmultvscl1);
294
295 fHistSPDmultvscl2 = new TH2F("fHistSPDmultvscl2", "Tracklet vs outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
296 fHistSPDmultvscl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
297 fHistSPDmultvscl2->GetYaxis()->SetTitle("Tracklet multiplicity");
298 fOutput->Add(fHistSPDmultvscl2);
299
ea23c361 300 fHistSPDphi = new TH1F("fHistSPDphi", "Tracklet #phi distribution", 360, 0.,2*TMath::Pi());
301 fHistSPDphi->GetXaxis()->SetTitle("#varphi [rad]");
302 fHistSPDphi->GetYaxis()->SetTitle("Entries");
303 fOutput->Add(fHistSPDphi);
304
4ab6aaf0 305 fHistSPDtheta = new TH1F("fHistSPDtheta", "Tracklet #theta distribution", 180, 0.,TMath::Pi());
ea23c361 306 fHistSPDtheta->GetXaxis()->SetTitle("#theta [rad]");
307 fHistSPDtheta->GetYaxis()->SetTitle("Entries");
308 fOutput->Add(fHistSPDtheta);
309
ea23c361 310 fHistSPDdePhi= new TH1F("fHistSPDdePhi", "Tracklet #Delta#varphi distribution",400,-0.1,.1);
311 fHistSPDdePhi->GetXaxis()->SetTitle("#Delta#varphi [rad]");
4ab6aaf0 312 fHistSPDdePhi->GetYaxis()->SetTitle("Entries");
ea23c361 313 fOutput->Add(fHistSPDdePhi);
314
ea23c361 315 fHistSPDphivsSPDeta= new TH2F("fHistSPDphivsSPDeta", "Tracklets - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());
316 fHistSPDphivsSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");
317 fHistSPDphivsSPDeta->GetYaxis()->SetTitle("#varphi [rad]");
318 fOutput->Add(fHistSPDphivsSPDeta);
319
4ab6aaf0 320 fHistSPDdeTheta= new TH1F("fHistSPDdeTheta", "Tracklet #Delta#theta distribution",100,-0.05,.05);
321 fHistSPDdeTheta->GetXaxis()->SetTitle("#Delta#theta [rad]");
322 fHistSPDdeTheta->GetYaxis()->SetTitle("Entries");
323 fOutput->Add(fHistSPDdeTheta);
324
86841765 325 fHistSPDvtx = new TH1F("fHistSPDvtx", "SPD vertex distribution",nBinVtx,-MaxVtx,MaxVtx);
326 fHistSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
327 fHistSPDvtx->GetYaxis()->SetTitle("Entries");
328 fOutput->Add(fHistSPDvtx);
329
330 fHistSPDvtxAnalysis = new TH1F("fHistSPDvtxAnalysis", "SPD vertex distribution: events selected for the analysis",nBinVtx,-MaxVtx,MaxVtx);
4ab6aaf0 331 fHistSPDvtxAnalysis->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
332 fHistSPDvtxAnalysis->GetYaxis()->SetTitle("Entries");
333 fOutput->Add(fHistSPDvtxAnalysis);
ea23c361 334
39f606ae 335 fHistSPDdePhideTheta= new TH2F("fHistSPDdePhideTheta", "Tracklet #Delta#varphi distribution",2000,-1.,1.,1000,-0.25,.25);
336 fHistSPDdePhideTheta->GetXaxis()->SetTitle("#Delta#varphi [rad]");
337 fHistSPDdePhideTheta->GetYaxis()->SetTitle("#Delta#theta [rad]");
338 fOutput->Add(fHistSPDdePhideTheta);
4ab6aaf0 339
86841765 340 fHistSPDphicl1 = new TH1F("fHistSPDphicl1", "Tracklet #phi distribution", 360, 0.,2*TMath::Pi());
341 fHistSPDphicl1->GetXaxis()->SetTitle("#varphi [rad]");
342 fHistSPDphicl1->GetYaxis()->SetTitle("Entries");
343 fOutput->Add(fHistSPDphicl1);
344
345 fHistSPDphicl2 = new TH1F("fHistSPDphicl2", "Tracklet #phi distribution", 360, 0.,2*TMath::Pi());
346 fHistSPDphicl2->GetXaxis()->SetTitle("#varphi [rad]");
347 fHistSPDphicl2->GetYaxis()->SetTitle("Entries");
348 fOutput->Add(fHistSPDphicl2);
349
ea23c361 350
39f606ae 351 if (fUseMC) {
ea23c361 352
39f606ae 353 // Track level correction histograms
86841765 354 fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
ea23c361 355 fOutput->Add(fHistBkgCorrDen);
ccfafc80 356
357 fHistReconstructedProtons = new TH2F("fHistReconstructedProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
358 fOutput->Add(fHistReconstructedProtons);
359 fHistReconstructedKaons = new TH2F("fHistReconstructedKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
360 fOutput->Add(fHistReconstructedKaons);
361 fHistReconstructedPions = new TH2F("fHistReconstructedPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
362 fOutput->Add(fHistReconstructedPions);
363 fHistReconstructedOthers = new TH2F("fHistReconstructedOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
364 fOutput->Add(fHistReconstructedOthers);
365
86841765 366 fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 367 fOutput->Add(fHistBkgCorrDenPrimGen);
ea23c361 368
ccfafc80 369 fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
370 fOutput->Add(fHistBkgCombLabels);
371
39f606ae 372 if (fTR) {
86841765 373 fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 374 fOutput->Add(fHistBkgCorrNum);
375
86841765 376 fHistAlgEffNum = new TH2F("fHistAlgEffNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 377 fOutput->Add(fHistAlgEffNum);
ea23c361 378
86841765 379 fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 380 fOutput->Add(fHistNonDetectableCorrDen);
ea23c361 381
39f606ae 382 }
ea23c361 383
86841765 384 fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 385 fOutput->Add(fHistNonDetectableCorrNum);
ccfafc80 386
387 fHistProtons = new TH2F("fHistProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
388 fOutput->Add(fHistProtons);
389 fHistKaons = new TH2F("fHistKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
390 fOutput->Add(fHistKaons);
391 fHistPions = new TH2F("fHistPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
392 fOutput->Add(fHistPions);
393 fHistOthers = new TH2F("fHistOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
394 fOutput->Add(fHistOthers);
395 fHistReconstructedSec = new TH2F("fHistReconstructedSec","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
396 fOutput->Add(fHistReconstructedSec);
397
398
86841765 399 fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 400 fOutput->Add(fHistAllPrimaries);
ea23c361 401
86841765 402 fHistTrackCentrEvts = new TH2F("fHistTrackCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 403 fOutput->Add(fHistTrackCentrEvts);
4ab6aaf0 404
86841765 405 fHistTrackTrigCentrEvts = new TH2F("fHistTrackTrigCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 406 fOutput->Add(fHistTrackTrigCentrEvts);
4ab6aaf0 407
408
f1fc7d61 409 // Event level correction histograms
86841765 410 fHistAllEvts = new TH2F("fHistAllEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 411 fOutput->Add(fHistAllEvts);
ea23c361 412
86841765 413 fHistCentrEvts = new TH2F("fHistCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 414 fOutput->Add(fHistCentrEvts);
ea23c361 415
86841765 416 fHistTrigCentrEvts = new TH2F("fHistTrigCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 417 fOutput->Add(fHistTrigCentrEvts);
ea23c361 418
86841765 419 fHistSelEvts = new TH2F("fHistSelEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
39f606ae 420 fOutput->Add(fHistSelEvts);
4ab6aaf0 421
4ab6aaf0 422
ea23c361 423 // MC distributions
39f606ae 424 fHistMCmultEtacut = new TH1F("fHistMCmultEtacut","Generated multiplicity",nBinMultCorr,binLimMultCorr);
425 fHistMCmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
426 fHistMCmultEtacut->GetYaxis()->SetTitle("Entries");
427 fOutput->Add(fHistMCmultEtacut);
428
429 fHistMCmultEtacutvsSPDmultEtacut = new TH2F("fHistMCmultEtacutvsSPDmultEtacut","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
430 fHistMCmultEtacutvsSPDmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
431 fHistMCmultEtacutvsSPDmultEtacut->GetYaxis()->SetTitle("Tracklet multiplicity |#eta|<1.4");
432 fOutput->Add(fHistMCmultEtacutvsSPDmultEtacut);
433
86841765 434 fHistMCmultEtacutvsSPDmultcl1 = new TH2F("fHistMCmultEtacutvsSPDmultcl1","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
435 fHistMCmultEtacutvsSPDmultcl1->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
436 fHistMCmultEtacutvsSPDmultcl1->GetYaxis()->SetTitle("Cluster inner layer multiplicity");
437 fOutput->Add(fHistMCmultEtacutvsSPDmultcl1);
438
439 fHistMCmultEtacutvsSPDmultcl2 = new TH2F("fHistMCmultEtacutvsSPDmultcl2","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
440 fHistMCmultEtacutvsSPDmultcl2->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
441 fHistMCmultEtacutvsSPDmultcl2->GetYaxis()->SetTitle("Cluster outer layer multiplicity");
442 fOutput->Add(fHistMCmultEtacutvsSPDmultcl2);
443
444
39f606ae 445 fHistMCvtxx = new TH1F("fHistMCvtxx", "MC vertex distribution - x",100,-.5,.5);
446 fOutput->Add(fHistMCvtxx);
447 fHistMCvtxy = new TH1F("fHistMCvtxy", "MC vertex distribution - y",100,-.5,.5);
448 fOutput->Add(fHistMCvtxy);
449 fHistMCvtxz = new TH1F("fHistMCvtxz", "MC vertex distribution - z",500,-50.,50.);
450 fOutput->Add(fHistMCvtxz);
451
452 fHistRecvsGenImpactPar= new TH2F("fHistRecvsGenImpactPar","",20,0.,20.,20,0.,20.);
453 fHistRecvsGenImpactPar->GetXaxis()->SetTitle("b_{gen} [fm]");
454 fHistRecvsGenImpactPar->GetYaxis()->SetTitle("b_{rec} [fm]");
455 fOutput->Add(fHistRecvsGenImpactPar);
456
457 fHistMCNpart = new TH1F("fHistMCNpart","Number of participants",450,0,450);
458 fHistMCNpart->GetXaxis()->SetTitle("N_{part}");
459 fHistMCNpart->GetYaxis()->SetTitle("Entries");
460 fOutput->Add(fHistMCNpart);
ea23c361 461
ccfafc80 462 fHistdPhidThetaPP = new TH2F("fHistdPhidThetaPP","",2000,-1.,1.,1000,-0.25,.25);
463 fOutput->Add(fHistdPhidThetaPP);
464 fHistdPhidThetaSS = new TH2F("fHistdPhidThetaSS","",2000,-1.,1.,1000,-0.25,.25);
465 fOutput->Add(fHistdPhidThetaSS);
466 fHistdPhidThetaComb = new TH2F("fHistdPhidThetaComb","",2000,-1.,1.,1000,-0.25,.25);
467 fOutput->Add(fHistdPhidThetaComb);
39f606ae 468
469 fHistDeVtx = new TH2F("fHistDeVtx","",80,-20.,20.,5000,-0.5,0.5);
470 fHistDeVtx->GetXaxis()->SetTitle("z_{MCvtx} [cm]");
471 fHistDeVtx->GetYaxis()->SetTitle("z_{MCvtx}-z_{SPDvtx} [cm]");
472 fOutput->Add(fHistDeVtx);
4ab6aaf0 473
ea23c361 474 }
39f606ae 475 PostData(1, fOutput);
476// Printf("Output objects created...");
ea23c361 477}
478
479//________________________________________________________________________
39f606ae 480void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *)
ea23c361 481{
482 // Main loop
39f606ae 483
ea23c361 484 // Called for each event
39f606ae 485// Printf("User exec..........");
ea23c361 486
39f606ae 487 AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
56df5f44 488 if (!hand) { printf("No RP handler\n"); return; }
ea23c361 489
39f606ae 490 fmyESD = dynamic_cast<AliESDEvent*>(InputEvent());
491 if (!fmyESD) {
ea23c361 492 Printf("ERROR: fESD not available");
493 return;
494 }
86841765 495
496 AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
497 if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
ea23c361 498
ccfafc80 499
39f606ae 500 // Trigger selection
501 Bool_t eventTriggered = kTRUE;
502 static AliTriggerAnalysis* triggerAnalysis = 0;
503 if (eventTriggered) // applying an offline trigger
504 eventTriggered = triggerAnalysis->IsTriggerFired(fmyESD, fTrigger);
505 if (!eventTriggered)
506 Printf("No trigger");
ea23c361 507
39f606ae 508
509 // Centrality selection
510 Bool_t eventInCentralityBin = kFALSE;
ccfafc80 511/* AliESDCentrality *centrality = fmyESD->GetCentrality();
39f606ae 512 if (fCentrEst=="") eventInCentralityBin = kTRUE;
4ab6aaf0 513 else {
39f606ae 514 if(!centrality) {
515 AliError("Centrality object not available");
516 } else {
517 if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
518 }
4ab6aaf0 519 }
ccfafc80 520*/
521 if (fCentrEst) {
522 AliESDVZERO* esdV0 = fmyESD->GetVZEROData();
523 Float_t multV0A=esdV0->GetMTotV0A();
524 Float_t multV0C=esdV0->GetMTotV0C();
525 fV0Ampl->Fill(multV0A+multV0C);
526 if (multV0A+multV0C>=fMinMultV0) eventInCentralityBin = kTRUE;
527 } else if (!fCentrEst) {
528 eventInCentralityBin = kTRUE;
529 }
4ab6aaf0 530
ccfafc80 531 const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
4ab6aaf0 532
39f606ae 533 // ESD vertex
534 Bool_t eventWithVertex = kFALSE;
535 const AliESDVertex* vtxESD = fmyESD->GetVertex();
ccfafc80 536 const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
39f606ae 537 Double_t esdvtx[3];
538 vtxESD->GetXYZ(esdvtx);
86841765 539 Int_t nContrib = vtxESD->GetNContributors();
ccfafc80 540 Int_t nContribTPC = vtxTPC->GetNContributors();
541 if (nContrib>0&&nContribTPC>0) {
542 if (vtxESD->GetDispersion()<0.04)
543 if (vtxESD->GetZRes()<0.25)
544 if (nContribTPC>(-10.+0.25*multESD->GetNumberOfITSClusters(0))) {
545 fHistSPDvtx->Fill(esdvtx[2]);
546 if (TMath::Abs(esdvtx[2])<fVtxLim)
547 eventWithVertex = kTRUE;
548 }
86841765 549 }
ea23c361 550
39f606ae 551 // Reconstructing or loading tracklets...
39f606ae 552 Int_t multSPD = multESD->GetNumberOfTracklets();
553 Int_t nSingleCl1 = multESD->GetNumberOfSingleClusters();
554 Int_t multSPDcl1 = nSingleCl1 + multSPD;
555 Int_t multSPDEtacut = 0;
556 Int_t multSPDAna = 0;
ea23c361 557
86841765 558 Int_t multSPDcl2 = multESD->GetNumberOfITSClusters(1);
ea23c361 559
39f606ae 560 Float_t trackletCoord[multSPDcl1][5];
561 Float_t trackletLab[multSPDcl1][2];
ea23c361 562
ccfafc80 563 Bool_t eventSelected = kFALSE;
86841765 564 // Event selection: in centrality bin, triggered with vertex
ccfafc80 565 if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
566 eventSelected = kTRUE;
86841765 567 fHistSPDmultcl1->Fill(multSPDcl1);
568 fHistSPDmultcl2->Fill(multSPDcl2);
569 fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2);
570 fHistSPDmultvscl1->Fill(multSPDcl1,multSPD);
571 fHistSPDmultvscl2->Fill(multSPDcl2,multSPD);
39f606ae 572 if (fRecoTracklets) {
573 // Load ITS rec points tree
574 TTree *itsClusterTree = hand->GetTreeR("ITS");
575 if (!itsClusterTree) {
576 AliError(" Invalid ITS cluster tree !\n");
577 return;
578 }
579 float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()};
580
86841765 581 fMultReco->SetPhiWindow(fPhiWindow);
582 fMultReco->SetThetaWindow(fThetaWindow);
583 fMultReco->SetPhiShift(fPhiShift);
584 fMultReco->SetRemoveClustersFromOverlaps(fRemoveClustersFromOverlaps);
585 fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
586 fMultReco->SetZetaOverlapCut(fZetaOverlapCut);
587 fMultReco->SetPhiRotationAngle(fPhiRotationAngle);
588 fMultReco->SetHistOn(kFALSE);
589
39f606ae 590 fMultReco->Reconstruct(itsClusterTree,vtxf,vtxf);
591// Printf("cl 1 from alg %d",fMultReco->GetNClustersLayer1());
592// Printf("cl 2 from alg %d",fMultReco->GetNClustersLayer2());
593 multSPD = fMultReco->GetNTracklets();
594// Printf("tracklets found...%d",multSPD);
595 nSingleCl1 = fMultReco->GetNSingleClusters();
596// Printf("singl found...%d",nSingleCl1);
597 for (Int_t itr = 0; itr<multSPD;++itr) {
598 trackletCoord[itr][3] = fMultReco->GetTracklet(itr)[0]; //theta
599 trackletCoord[itr][2] = fMultReco->GetTracklet(itr)[1]; //phi
600 trackletCoord[itr][1] = fMultReco->GetTracklet(itr)[3]; //deTheta
601 trackletCoord[itr][0] = fMultReco->GetTracklet(itr)[2]; //dephi
86841765 602 trackletCoord[itr][4] = -TMath::Log(TMath::Tan(trackletCoord[itr][3]/2.)); //eta
603 trackletLab[itr][0] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[4]); //label lay1
604 trackletLab[itr][1] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[5]); //label lay2
605 }
606 for (Int_t icl1 = 0; icl1<multSPDcl1;++icl1) {
607 fHistSPDphicl1->Fill(fMultReco->GetClusterLayer1(icl1)[1]);
608 }
609 for (Int_t icl2 = 0; icl2<multSPDcl2; ++icl2) {
610 fHistSPDphicl2->Fill(fMultReco->GetClusterLayer2(icl2)[1]);
39f606ae 611 }
86841765 612
39f606ae 613 } else {
614 // set variables from ESD
615 for (Int_t itr = 0; itr<multSPD;++itr) {
616 trackletCoord[itr][3] = multESD->GetTheta(itr); //theta
617 trackletCoord[itr][2] = multESD->GetPhi(itr); //phi
618 trackletCoord[itr][1] = multESD->GetDeltaTheta(itr); //deTheta
619 trackletCoord[itr][0] = multESD->GetDeltaPhi(itr); //dePhi
86841765 620 trackletCoord[itr][4] = multESD->GetEta(itr); //eta
39f606ae 621 trackletLab[itr][0] = multESD->GetLabel(itr,0); //label lay1
622 trackletLab[itr][1] = multESD->GetLabel(itr,1); //label lay2
623 }
624 }
625// Printf("tracklets in ESD...%d",multESD->GetNumberOfTracklets());
ea23c361 626
627 for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
86841765 628 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) { // select tracklets with tighter cuts
39f606ae 629 fHistSPDRAWEtavsZ->Fill(trackletCoord[itracklet][4],esdvtx[2]);
630 fHistSPDphi->Fill(trackletCoord[itracklet][2]);
631 fHistSPDtheta->Fill(trackletCoord[itracklet][3]);
632 fHistSPDdePhi->Fill(trackletCoord[itracklet][0]);
633 fHistSPDdeTheta->Fill(trackletCoord[itracklet][1]);
4ab6aaf0 634
39f606ae 635 fHistSPDphivsSPDeta->Fill(trackletCoord[itracklet][4],trackletCoord[itracklet][2]);
636 multSPDAna++;
637 // Calculate multiplicity in etacut
638 if (TMath::Abs(trackletCoord[itracklet][4])<1.4) multSPDEtacut++;
639 }
640 // for combinatorial background normalization
641 fHistSPDdePhideTheta->Fill(trackletCoord[itracklet][0], trackletCoord[itracklet][1]);
ea23c361 642
ea23c361 643 }
39f606ae 644 if (multSPDAna!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
ea23c361 645 fHistSPDmultEtacut->Fill(multSPDEtacut);
39f606ae 646 fHistSPDmult->Fill(multSPDAna);
647
648 fHistSPDRAWMultvsZ->Fill(multSPDAna,esdvtx[2]);
ea23c361 649
f1fc7d61 650 } // End selected events
ea23c361 651
39f606ae 652 if (eventInCentralityBin) {
653 fHistSPDRAWMultvsZCentrEvts->Fill(multSPDAna,esdvtx[2]);
654 if (eventTriggered) fHistSPDRAWMultvsZTriggCentrEvts->Fill(multSPDAna,esdvtx[2]);
4ab6aaf0 655 }
ea23c361 656
39f606ae 657 if (fUseMC) {
ea23c361 658
39f606ae 659 if (!fMCEvent) {
660 AliError("No MC info found");
661 return;
662 }
ea23c361 663 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
664 if (!eventHandler) {
665 Printf("ERROR: Could not retrieve MC event handler");
666 return;
667 }
668
39f606ae 669 AliStack* stack = fMCEvent->Stack();
ea23c361 670 if (!stack) {
671 AliDebug(AliLog::kError, "Stack not available");
672 return;
673 }
674
ccfafc80 675 Float_t impactParameter = 0.;
676 Int_t npart = 0;
677
39f606ae 678 AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
ccfafc80 679 AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
680
681 if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
682 else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
683
39f606ae 684 Bool_t IsEventInMCCentralityBin = kFALSE;
685 switch (fMCCentralityBin) {
686
86841765 687 case 1: if (impactParameter<3) IsEventInMCCentralityBin = kTRUE;
688 break;
689 case 2: if (impactParameter>3&&impactParameter<6) IsEventInMCCentralityBin = kTRUE;
690 break;
691 case 3: if (impactParameter>6&&impactParameter<9) IsEventInMCCentralityBin = kTRUE;
692 break;
693 case 4: if (impactParameter>9&&impactParameter<12) IsEventInMCCentralityBin = kTRUE;
694 break;
695 case 5: if (impactParameter>12&&impactParameter<15) IsEventInMCCentralityBin = kTRUE;
696 break;
697 case 6: if (impactParameter>15) IsEventInMCCentralityBin = kTRUE;
698 break;
699 case 7: IsEventInMCCentralityBin = kTRUE;
700 break;
b8819a80 701
4ab6aaf0 702 }
39f606ae 703
704 if (IsEventInMCCentralityBin) {
86841765 705 // MC vertex
706 TArrayF vtxMC(3);
707 fMCEvent->GenEventHeader()->PrimaryVertex(vtxMC);
708 fHistMCvtxx->Fill(vtxMC[0]);
709 fHistMCvtxy->Fill(vtxMC[1]);
710 fHistMCvtxz->Fill(vtxMC[2]);
ccfafc80 711
86841765 712// Printf("Impact parameter gen: %f", impactParameter);
ccfafc80 713 if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
714 else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
715
86841765 716 //Rec centrality vs gen centrality
717 AliESDZDC *zdcRec = fmyESD->GetESDZDC();
718 Double_t impactParameterZDC = zdcRec->GetImpactParameter();
719// Printf("Impact parameter rec: %f", impactParameterZDC);
720
721 fHistRecvsGenImpactPar->Fill(impactParameter,impactParameterZDC);
722 fHistMCNpart->Fill(npart);
723
724 // Tracks from MC
725 Int_t multMCCharged = 0;
726 Int_t multMCChargedEtacut = 0;
ccfafc80 727// Int_t nMCPart = stack->GetNprimary();
728 Int_t nMCPart = stack->GetNtrack(); // decay products of D and B mesons are also primaries and produced in HIJING during transport
86841765 729 Float_t* etagen = new Float_t[nMCPart];
730 Int_t* stackIndexOfPrimaryParts = new Int_t[nMCPart];
731 Bool_t* reconstructedPrimaryPart = new Bool_t[nMCPart];
732 Bool_t* detectedPrimaryPartLay1 = new Bool_t[nMCPart];
733 Bool_t* detectedPrimaryPartLay2 = new Bool_t[nMCPart];
734 Bool_t* detectablePrimaryPart = new Bool_t[nMCPart];
86841765 735
736 TTree* tRefTree;
39f606ae 737 if (fTR) {
86841765 738 tRefTree = eventHandler->TreeTR();
739 fMCEvent->ConnectTreeTR(tRefTree);
39f606ae 740 }
ea23c361 741
86841765 742 // Loop over MC particles
743 for (Int_t imc=0; imc<nMCPart; imc++) {
744 AliMCParticle *mcpart = (AliMCParticle*)fMCEvent->GetTrack(imc);
ccfafc80 745
86841765 746 Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
ccfafc80 747 if (!isPrimary) continue;
86841765 748 if (mcpart->Charge() == 0) continue;
749 Float_t theta = mcpart->Theta();
750 if (theta==0 || theta==TMath::Pi()) continue;
751 Float_t eta = mcpart->Eta();
752// Float_t pt = mcpart->Pt();
753 etagen[multMCCharged] = eta;
754 stackIndexOfPrimaryParts[multMCCharged] = imc;
755
756 reconstructedPrimaryPart[multMCCharged]=kFALSE;
757 detectedPrimaryPartLay1[multMCCharged]=kFALSE;
758 detectedPrimaryPartLay2[multMCCharged]=kFALSE;
759 detectablePrimaryPart[multMCCharged]=kFALSE;
ea23c361 760
39f606ae 761 if (fTR) {
86841765 762 Int_t nref = mcpart->GetNumberOfTrackReferences();
763 if (nref==0) {
764 detectablePrimaryPart[multMCCharged]= kTRUE;
765 } else if (nref>0) {
766 // Check if the primary is detectable
767 detectablePrimaryPart[multMCCharged] = IsDetectablePrimary(nref,mcpart);
768 // Check if the primary is detected (if SPD 100% efficient)
769 if (detectablePrimaryPart[multMCCharged]) {
770 detectedPrimaryPartLay1[multMCCharged] = IsDetectedPrimary(nref,mcpart,0);
771 detectedPrimaryPartLay2[multMCCharged] = IsDetectedPrimary(nref,mcpart,1);
772 }
773 }
ea23c361 774 }
ccfafc80 775 if (eventSelected&&fPartSpecies) {
776 if (TMath::Abs(mcpart->PdgCode())==2212) fHistProtons->Fill(etagen[multMCCharged],vtxMC[2]);
777 else if (TMath::Abs(mcpart->PdgCode())==321) fHistKaons->Fill(etagen[multMCCharged],vtxMC[2]);
778 else if (TMath::Abs(mcpart->PdgCode())==211) fHistPions->Fill(etagen[multMCCharged],vtxMC[2]);
779 else fHistOthers->Fill(etagen[multMCCharged],vtxMC[2]); //includes leptons pdg->11,13
780 }
86841765 781 multMCCharged++;
782 if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
783 } // end of MC particle loop
ea23c361 784
86841765 785 fHistMCmultEtacut->Fill(multMCChargedEtacut);
4ab6aaf0 786
86841765 787 // Event selection: in centrality bin, triggered with vertex
ccfafc80 788 if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
86841765 789
790 fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);
791
792 for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
ccfafc80 793 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna)
86841765 794 fHistBkgCorrDen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
795
ccfafc80 796 if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) {
797 Bool_t trakletByPrim = kFALSE;
798 for (Int_t imc=0; imc<multMCCharged; imc++) {
799 if (trackletLab[itracklet][0]==stackIndexOfPrimaryParts[imc]) {
800 fHistdPhidThetaPP->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
801 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
86841765 802 fHistBkgCorrDenPrimGen->Fill(etagen[imc],vtxMC[2]);
ccfafc80 803 if (fPartSpecies) {
804 if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==2212)
805 fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
806 else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==321)
807 fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
808 else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==211)
809 fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
810 else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
811 }
812 }
813 trakletByPrim = kTRUE;
814 if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
815 reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers
816 }
817 break;
818 }
819 }
820 if (!trakletByPrim) {
821 fHistdPhidThetaSS->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
822 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
823
86841765 824 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
ccfafc80 825 if (fPartSpecies) {
826 Int_t motherlab = ((AliMCParticle*)fMCEvent->GetTrack((Int_t)trackletLab[itracklet][0]))->GetMother();
827 if (motherlab>-1) {
828 if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==2212) fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
829 else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==321) fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
830 else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==211) fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
831 else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
832 } else fHistReconstructedSec->Fill(trackletCoord[itracklet][4],esdvtx[2]);
833 }
86841765 834 }
ccfafc80 835 }
836 } else {
837 fHistdPhidThetaComb->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
838 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
86841765 839 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
ccfafc80 840 fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]);
86841765 841 }
ccfafc80 842 }
86841765 843 } // end loop tracklets
844
845 for (Int_t imc=0; imc<multMCCharged; imc++) {
846 if (fTR) {
847 if (reconstructedPrimaryPart[imc]) fHistBkgCorrNum->Fill(etagen[imc],vtxMC[2]); // only for separate corrections
848 if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) fHistAlgEffNum->Fill(etagen[imc],vtxMC[2]);
849 if (detectablePrimaryPart[imc]) fHistNonDetectableCorrDen->Fill(etagen[imc],vtxMC[2]);
850 }
851 fHistNonDetectableCorrNum->Fill(etagen[imc],vtxMC[2]);
852 }
4ab6aaf0 853
86841765 854 fHistSelEvts->Fill(multSPDAna,vtxMC[2]);
855 } // end of selected events
ea23c361 856
86841765 857 fHistMCmultEtacutvsSPDmultEtacut->Fill(multMCChargedEtacut,multSPDEtacut);
858 fHistMCmultEtacutvsSPDmultcl1->Fill(multMCChargedEtacut,multSPDcl1);
859 fHistMCmultEtacutvsSPDmultcl2->Fill(multMCChargedEtacut,multSPDcl2);
860 fHistAllEvts->Fill(multSPDAna,vtxMC[2]);
ea23c361 861
39f606ae 862 if (eventInCentralityBin) {
86841765 863 fHistCentrEvts->Fill(multSPDAna,vtxMC[2]);
864 if (eventTriggered) {
865 fHistTrigCentrEvts->Fill(multSPDAna,vtxMC[2]);
866 }
39f606ae 867 }
86841765 868
869 for (Int_t imc=0; imc<multMCCharged; imc++) {
870 fHistAllPrimaries->Fill(etagen[imc],vtxMC[2]);
871 if (eventInCentralityBin) {
872 fHistTrackCentrEvts->Fill(etagen[imc],vtxMC[2]);
873 if (eventTriggered) fHistTrackTrigCentrEvts->Fill(etagen[imc],vtxMC[2]);
874 }
875 }
39f606ae 876
86841765 877 delete[] etagen;
878 delete[] stackIndexOfPrimaryParts;
879 delete[] reconstructedPrimaryPart;
880 delete[] detectedPrimaryPartLay1;
881 delete[] detectedPrimaryPartLay2;
882 delete[] detectablePrimaryPart;
883 }
ea23c361 884 }
39f606ae 885 PostData(1, fOutput);
ea23c361 886}
ea23c361 887//________________________________________________________________________
86841765 888Bool_t AliAnalysisTaskSPDdNdEta::IsDetectedPrimary(Int_t nref, AliMCParticle* mcpart, Int_t Layer) {
39f606ae 889 Double_t rMinLay1 = 3.4; //min rad for track ref first SPD layer
890 Double_t rMaxLay1 = 4.4; //max rad for track ref first SPD layer
891 Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer
892 Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer
893 Double_t rmin = (Layer > 0 ? rMinLay2 : rMinLay1);
894 Double_t rmax = (Layer > 0 ? rMaxLay2 : rMaxLay1);
895 AliTrackReference *tref=0x0;
896
897 for (Int_t iref=0;iref<nref;iref++) {
898 tref = mcpart->GetTrackReference(iref);
899 if (tref) {
900 if (tref->R()>rmin&&tref->R()<rmax) {
901 if (tref->DetectorId()==0) {
902 return kTRUE;
903 }
904 }
905 }
ea23c361 906 }
39f606ae 907 return kFALSE;
908}
909//________________________________________________________________________
86841765 910Bool_t AliAnalysisTaskSPDdNdEta::IsDetectablePrimary(Int_t nref, AliMCParticle* mcpart) {
911
39f606ae 912 Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer
913 Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer
914
915 AliTrackReference *tref= mcpart->GetTrackReference(nref-1);
916 if (tref->DetectorId()!=-1) {
917 return kTRUE;
918 } else { //last is -1 -> particle disappeared. Where?
919 if (tref->R()>rMaxLay2) {
920 return kTRUE;
921 } else if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2) { // this last tr is in lay 2
922 for (Int_t iref=0;iref<nref;iref++) { // look for other tr in lay 2
923 tref = mcpart->GetTrackReference(iref);
924 if (tref)
925 if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2)
926 if (tref->DetectorId()==0) return kTRUE;
927 }
928 } // else particle disappeared before lay2
ea23c361 929 }
39f606ae 930 return kFALSE;
931}
932//________________________________________________________________________
933void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
934{
935 Printf("Terminating...");
936// AliAnalysisTaskSE::Terminate();
937
a083f34c 938}