Coverity fixes
[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliAnalysisTaskSPDdNdEta.cxx
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
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 //////////////////////////////////////////////////////////////////////////
23
24 #include "TChain.h"
25 #include "TTree.h"
26
27 #include "TH1F.h"
28 #include "TH2F.h" 
29 #include "TH3F.h" 
30 #include "TList.h"
31 #include "TGeoGlobalMagField.h"
32 #include "AliMagF.h"
33 #include "AliAnalysisManager.h"
34
35 #include "AliMultiplicity.h"
36 #include "AliESDEvent.h"  
37 #include "AliESDInputHandler.h"
38
39 #include "AliESDInputHandlerRP.h"
40 #include "AliESDVZERO.h"
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
48 #include "AliMCEventHandler.h"
49 #include "AliMCEvent.h"
50 #include "AliStack.h"
51 #include "AliTrackReference.h"
52
53 #include "AliGenHijingEventHeader.h" 
54 #include "AliGenDPMjetEventHeader.h"
55
56 #include "AliLog.h"
57
58 #include "AliTriggerAnalysis.h" 
59 #include "AliPhysicsSelection.h"
60 #include "AliTrackletAlg.h" 
61 #include "AliAnalysisTaskSPDdNdEta.h"
62
63
64 ClassImp(AliAnalysisTaskSPDdNdEta)
65 //________________________________________________________________________
66 AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name) 
67   : AliAnalysisTaskSE(name), 
68
69   fmyESD(0), 
70   fOutput(0), 
71   
72   fUseMC(kFALSE), 
73   fTrigger(AliTriggerAnalysis::kAcceptAll),
74   fTR(kFALSE),
75   fRecoTracklets(kFALSE),
76
77   fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall),
78   fCentrLowLim(0),
79   fCentrUpLim(0),
80   fCentrEst(kFALSE),
81   fMinClMultLay2(0),
82   fMaxClMultLay2(0),
83   fMinMultV0(0),
84   fVtxLim(0),
85   fPartSpecies(0),
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
98   fV0Ampl(0),
99   fHistSPDRAWMultvsZ(0),
100   fHistSPDRAWMultvsZTriggCentrEvts(0),
101   fHistSPDRAWMultvsZCentrEvts(0),
102   fHistSPDRAWEtavsZ(0),
103  
104   fHistSPDmultEtacut(0),
105   fHistSPDmult(0),
106   fHistSPDmultcl1(0),
107   fHistSPDmultcl2(0),
108   fHistSPDmultcl1vscl2(0),
109   fHistSPDmultvscl1(0),
110   fHistSPDmultvscl2(0),
111
112   fHistSPDeta(0),
113   fHistSPDphi(0),
114   fHistSPDtheta(0),
115   fHistSPDdePhi(0),
116   fHistSPDphivsSPDeta(0),  
117   fHistSPDdeTheta(0),
118   fHistSPDvtx(0),
119   fHistSPDvtxAnalysis(0),
120   fHistSPDdePhideTheta(0),
121
122   fHistSPDphicl1(0),
123   fHistSPDphicl2(0),
124
125   fHistBkgCorrDen(0),
126   fHistReconstructedProtons(0),
127   fHistReconstructedKaons(0),
128   fHistReconstructedPions(0),
129   fHistReconstructedOthers(0), 
130   fHistReconstructedSec(0),
131   fHistBkgCorrDenPrimGen(0),
132   fHistBkgCombLabels(0),
133   fHistBkgCorrNum(0),
134   fHistAlgEffNum(0),
135   fHistNonDetectableCorrDen(0),
136
137   fHistNonDetectableCorrNum(0),
138   fHistProtons(0),
139   fHistKaons(0),
140   fHistPions(0),
141   fHistOthers(0),
142   fHistAllPrimaries(0),
143   fHistTrackCentrEvts(0),
144   fHistTrackTrigCentrEvts(0),
145
146   fHistAllEvts(0),
147   fHistCentrEvts(0),
148   fHistTrigCentrEvts(0),
149   fHistSelEvts(0),
150
151   fHistMCmultEtacut(0),
152   fHistMCmultEtacutvsSPDmultEtacut(0),
153   fHistMCmultEtacutvsSPDmultcl1(0),
154   fHistMCmultEtacutvsSPDmultcl2(0),
155  
156   fHistMCvtxx(0),
157   fHistMCvtxy(0),
158   fHistMCvtxz(0),
159
160   fHistRecvsGenImpactPar(0),
161   fHistMCNpart(0), 
162
163   fHistdPhidThetaPP(0),
164   fHistdPhidThetaSS(0),
165   fHistdPhidThetaComb(0),
166
167   fHistDeVtx(0)
168
169 {
170
171   // Constructor
172
173   // Define input and output slots here
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());  
178
179 }
180 //________________________________________________________________________
181 AliAnalysisTaskSPDdNdEta::~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
188   delete fMultReco;
189
190   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
191     delete fOutput;
192     fOutput = 0;
193   }
194
195 }
196 //________________________________________________________________________
197 void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects() 
198 {
199   if (fRecoTracklets) {
200     AliCDBManager *man = AliCDBManager::Instance();
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");
204     man->SetRun(137161); 
205     AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
206     if(!obj) AliFatal("Unable to load geometry from CDB!");
207     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
208     AliGeomManager::GetNalignable("ITS");
209     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
210   }
211
212   fMultReco = new AliTrackletAlg();
213
214   // Create histograms
215   fOutput = new TList();
216   fOutput->SetOwner(); 
217
218   Int_t nBinVtx = 40;
219   Double_t MaxVtx = 20.;
220
221   Int_t nBinMultCorr = 200;
222   Float_t nMaxMult = 20000.;
223   Double_t binLimMultCorr[nBinMultCorr+1];
224   binLimMultCorr[0] = 0.;
225   binLimMultCorr[1] = 1.;
226   for (Int_t i = 2; i<=nBinMultCorr;++i) {
227     binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
228   }
229
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   } 
237
238   // Data to be corrected
239   // ...event level    
240   
241   fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
242   fOutput->Add(fV0Ampl);
243
244   fHistSPDRAWMultvsZ = new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
245   fHistSPDRAWMultvsZ->GetXaxis()->SetTitle("Tracklet multiplicity");
246   fHistSPDRAWMultvsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
247   fOutput->Add(fHistSPDRAWMultvsZ); 
248
249   fHistSPDRAWMultvsZTriggCentrEvts = new TH2F("fHistSPDRAWMultvsZTriggCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
250   fHistSPDRAWMultvsZTriggCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
251   fHistSPDRAWMultvsZTriggCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
252   fOutput->Add(fHistSPDRAWMultvsZTriggCentrEvts);
253
254   fHistSPDRAWMultvsZCentrEvts = new TH2F("fHistSPDRAWMultvsZCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
255   fHistSPDRAWMultvsZCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
256   fHistSPDRAWMultvsZCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
257   fOutput->Add(fHistSPDRAWMultvsZCentrEvts);
258
259   // ...track level
260   fHistSPDRAWEtavsZ = new TH2F("fHistSPDRAWEtavsZ", "Tracklet pseudorapidity distribution", nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
261   fHistSPDRAWEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
262   fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
263   fOutput->Add(fHistSPDRAWEtavsZ);
264   
265   fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution",nBinMultCorr,binLimMultCorr);
266   fHistSPDmultEtacut->GetXaxis()->SetTitle("Tracklet multiplicity (|#eta|<1.4)");
267   fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
268   fOutput->Add(fHistSPDmultEtacut);
269
270   fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", nBinMultCorr,binLimMultCorr);
271   fHistSPDmult->GetXaxis()->SetTitle("Tracklet multiplicity");
272   fHistSPDmult->GetYaxis()->SetTitle("Entries");
273   fOutput->Add(fHistSPDmult);
274
275   fHistSPDmultcl1 = new TH1F("fHistSPDmultcl1", "Inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
276   fHistSPDmultcl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
277   fHistSPDmultcl1->GetYaxis()->SetTitle("Entries");
278   fOutput->Add(fHistSPDmultcl1);
279
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
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
305   fHistSPDtheta = new TH1F("fHistSPDtheta", "Tracklet #theta  distribution", 180, 0.,TMath::Pi());
306   fHistSPDtheta->GetXaxis()->SetTitle("#theta [rad]");
307   fHistSPDtheta->GetYaxis()->SetTitle("Entries");
308   fOutput->Add(fHistSPDtheta);
309
310   fHistSPDdePhi= new TH1F("fHistSPDdePhi", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
311   fHistSPDdePhi->GetXaxis()->SetTitle("#Delta#varphi [rad]");
312   fHistSPDdePhi->GetYaxis()->SetTitle("Entries");
313   fOutput->Add(fHistSPDdePhi);
314
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
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
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);
331   fHistSPDvtxAnalysis->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
332   fHistSPDvtxAnalysis->GetYaxis()->SetTitle("Entries");
333   fOutput->Add(fHistSPDvtxAnalysis);
334
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);
339
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
350
351   if (fUseMC) {
352
353     // Track level correction histograms 
354     fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
355     fOutput->Add(fHistBkgCorrDen);
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  
366     fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
367     fOutput->Add(fHistBkgCorrDenPrimGen);
368
369     fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
370     fOutput->Add(fHistBkgCombLabels);
371
372     if (fTR) {
373       fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
374       fOutput->Add(fHistBkgCorrNum);
375  
376       fHistAlgEffNum = new TH2F("fHistAlgEffNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
377       fOutput->Add(fHistAlgEffNum);  
378
379       fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
380       fOutput->Add(fHistNonDetectableCorrDen);
381
382     }
383
384     fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
385     fOutput->Add(fHistNonDetectableCorrNum);
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  
399     fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
400     fOutput->Add(fHistAllPrimaries);
401
402     fHistTrackCentrEvts = new TH2F("fHistTrackCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
403     fOutput->Add(fHistTrackCentrEvts);
404
405     fHistTrackTrigCentrEvts = new TH2F("fHistTrackTrigCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
406     fOutput->Add(fHistTrackTrigCentrEvts);
407
408
409     // Event level correction histograms  
410     fHistAllEvts = new TH2F("fHistAllEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
411     fOutput->Add(fHistAllEvts);
412
413     fHistCentrEvts = new TH2F("fHistCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
414     fOutput->Add(fHistCentrEvts);
415
416     fHistTrigCentrEvts = new TH2F("fHistTrigCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
417     fOutput->Add(fHistTrigCentrEvts);
418
419     fHistSelEvts = new TH2F("fHistSelEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
420     fOutput->Add(fHistSelEvts);
421
422
423     // MC distributions
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
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
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);
461  
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);
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);
473
474   }
475   PostData(1, fOutput);
476 //  Printf("Output objects created...");
477 }
478
479 //________________________________________________________________________
480 void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 
481 {
482   // Main loop
483
484   // Called for each event
485 //  Printf("User exec..........");
486
487   AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
488   if (!hand) { printf("No RP handler\n"); return; }
489
490   fmyESD = dynamic_cast<AliESDEvent*>(InputEvent()); 
491   if (!fmyESD) {
492     Printf("ERROR: fESD not available");
493     return;
494   }
495     
496   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
497   if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
498
499
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");
507
508
509   // Centrality selection 
510   Bool_t eventInCentralityBin = kFALSE;
511 /*  AliESDCentrality *centrality = fmyESD->GetCentrality();
512   if (fCentrEst=="") eventInCentralityBin = kTRUE;
513   else {
514     if(!centrality) {
515       AliError("Centrality object not available"); 
516     }  else {
517       if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
518     }
519   }
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   }
530
531   const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
532
533   // ESD vertex
534   Bool_t eventWithVertex = kFALSE;
535   const AliESDVertex* vtxESD = fmyESD->GetVertex();
536   const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
537   Double_t esdvtx[3];
538   vtxESD->GetXYZ(esdvtx);
539   Int_t nContrib = vtxESD->GetNContributors();
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         }
549   } 
550
551   // Reconstructing or loading tracklets...
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;
557
558   Int_t multSPDcl2 = multESD->GetNumberOfITSClusters(1);
559
560   Float_t trackletCoord[multSPDcl1][5];
561   Float_t trackletLab[multSPDcl1][2];
562
563   Bool_t eventSelected = kFALSE;
564   // Event selection: in centrality bin, triggered with vertex
565   if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
566     eventSelected = kTRUE; 
567     fHistSPDmultcl1->Fill(multSPDcl1);
568     fHistSPDmultcl2->Fill(multSPDcl2);
569     fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2);
570     fHistSPDmultvscl1->Fill(multSPDcl1,multSPD);
571     fHistSPDmultvscl2->Fill(multSPDcl2,multSPD);
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
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
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
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]);          
611       }
612
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
620         trackletCoord[itr][4] = multESD->GetEta(itr);         //eta 
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());
626
627     for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
628       if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) { // select tracklets with tighter cuts
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]);
634  
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]);
642
643     }
644     if (multSPDAna!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
645     fHistSPDmultEtacut->Fill(multSPDEtacut);
646     fHistSPDmult->Fill(multSPDAna); 
647
648     fHistSPDRAWMultvsZ->Fill(multSPDAna,esdvtx[2]);
649
650   } // End selected events
651
652   if (eventInCentralityBin) {
653     fHistSPDRAWMultvsZCentrEvts->Fill(multSPDAna,esdvtx[2]);  
654     if (eventTriggered) fHistSPDRAWMultvsZTriggCentrEvts->Fill(multSPDAna,esdvtx[2]);
655   }
656
657   if (fUseMC) {
658
659     if (!fMCEvent) {
660       AliError("No MC info found");
661       return;
662     } 
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
669     AliStack* stack = fMCEvent->Stack();
670     if (!stack) {
671       AliDebug(AliLog::kError, "Stack not available");
672       return;
673     }
674
675     Float_t impactParameter = 0.;
676     Int_t npart = 0;
677
678     AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
679     AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
680
681     if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
682     else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
683
684     Bool_t IsEventInMCCentralityBin = kFALSE;
685     switch (fMCCentralityBin) {
686
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;
701
702     }
703
704     if (IsEventInMCCentralityBin) {
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]);
711
712 //      Printf("Impact parameter gen: %f", impactParameter);
713        if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
714        else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
715
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;
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
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];
735
736       TTree* tRefTree;  
737       if (fTR) {
738         tRefTree = eventHandler->TreeTR(); 
739         fMCEvent->ConnectTreeTR(tRefTree); 
740       }
741
742       // Loop over MC particles
743       for (Int_t imc=0; imc<nMCPart; imc++) {
744         AliMCParticle *mcpart  = (AliMCParticle*)fMCEvent->GetTrack(imc);
745
746         Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
747         if (!isPrimary)            continue;
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;
760
761         if (fTR) {
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           }  
774         }
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         }
781         multMCCharged++;
782         if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
783       } // end of MC particle loop
784
785       fHistMCmultEtacut->Fill(multMCChargedEtacut);
786
787       // Event selection: in centrality bin, triggered with vertex 
788       if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
789     
790         fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);      
791
792         for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
793           if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) 
794             fHistBkgCorrDen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
795
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) {
802                   fHistBkgCorrDenPrimGen->Fill(etagen[imc],vtxMC[2]);
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
824                 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
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                 }
834               }
835             }
836           } else {
837             fHistdPhidThetaComb->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
838             if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
839               fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
840               fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]); 
841             }
842           }
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         }
853
854         fHistSelEvts->Fill(multSPDAna,vtxMC[2]);
855       } // end of selected events
856
857       fHistMCmultEtacutvsSPDmultEtacut->Fill(multMCChargedEtacut,multSPDEtacut);
858       fHistMCmultEtacutvsSPDmultcl1->Fill(multMCChargedEtacut,multSPDcl1);
859       fHistMCmultEtacutvsSPDmultcl2->Fill(multMCChargedEtacut,multSPDcl2);
860       fHistAllEvts->Fill(multSPDAna,vtxMC[2]);
861
862       if (eventInCentralityBin) {
863         fHistCentrEvts->Fill(multSPDAna,vtxMC[2]);
864         if (eventTriggered) {
865           fHistTrigCentrEvts->Fill(multSPDAna,vtxMC[2]);
866         }
867       }
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       } 
876     
877       delete[] etagen;
878       delete[] stackIndexOfPrimaryParts;
879       delete[] reconstructedPrimaryPart;
880       delete[] detectedPrimaryPartLay1;
881       delete[] detectedPrimaryPartLay2;
882       delete[] detectablePrimaryPart;
883     }  
884   }
885   PostData(1, fOutput);
886 }      
887 //________________________________________________________________________
888 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectedPrimary(Int_t nref, AliMCParticle* mcpart, Int_t Layer) {
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      }
906   }
907   return kFALSE;
908 }
909 //________________________________________________________________________
910 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectablePrimary(Int_t nref, AliMCParticle* mcpart) {
911
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
929   }
930   return kFALSE; 
931 }
932 //________________________________________________________________________
933 void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *) 
934 {
935   Printf("Terminating...");
936 //  AliAnalysisTaskSE::Terminate();
937
938 }