AliESDCentrality replaced by AliCentrality
[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     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
207     AliGeomManager::GetNalignable("ITS");
208     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
209   }
210
211   fMultReco = new AliTrackletAlg();
212
213   // Create histograms
214   fOutput = new TList();
215   fOutput->SetOwner(); 
216
217   Int_t nBinVtx = 40;
218   Double_t MaxVtx = 20.;
219
220   Int_t nBinMultCorr = 200;
221   Float_t nMaxMult = 20000.;
222   Double_t binLimMultCorr[nBinMultCorr+1];
223   binLimMultCorr[0] = 0.;
224   binLimMultCorr[1] = 1.;
225   for (Int_t i = 2; i<=nBinMultCorr;++i) {
226     binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
227   }
228
229   Int_t nBinEtaCorr = 60; 
230   Float_t etaMax = 3.;
231   Float_t etaMin = -3.;
232   Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
233   for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
234     binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
235   } 
236
237   // Data to be corrected
238   // ...event level    
239   
240   fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
241   fOutput->Add(fV0Ampl);
242
243   fHistSPDRAWMultvsZ = new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
244   fHistSPDRAWMultvsZ->GetXaxis()->SetTitle("Tracklet multiplicity");
245   fHistSPDRAWMultvsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
246   fOutput->Add(fHistSPDRAWMultvsZ); 
247
248   fHistSPDRAWMultvsZTriggCentrEvts = new TH2F("fHistSPDRAWMultvsZTriggCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
249   fHistSPDRAWMultvsZTriggCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
250   fHistSPDRAWMultvsZTriggCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
251   fOutput->Add(fHistSPDRAWMultvsZTriggCentrEvts);
252
253   fHistSPDRAWMultvsZCentrEvts = new TH2F("fHistSPDRAWMultvsZCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
254   fHistSPDRAWMultvsZCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
255   fHistSPDRAWMultvsZCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
256   fOutput->Add(fHistSPDRAWMultvsZCentrEvts);
257
258   // ...track level
259   fHistSPDRAWEtavsZ = new TH2F("fHistSPDRAWEtavsZ", "Tracklet pseudorapidity distribution", nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
260   fHistSPDRAWEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
261   fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
262   fOutput->Add(fHistSPDRAWEtavsZ);
263   
264   fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution",nBinMultCorr,binLimMultCorr);
265   fHistSPDmultEtacut->GetXaxis()->SetTitle("Tracklet multiplicity (|#eta|<1.4)");
266   fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
267   fOutput->Add(fHistSPDmultEtacut);
268
269   fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", nBinMultCorr,binLimMultCorr);
270   fHistSPDmult->GetXaxis()->SetTitle("Tracklet multiplicity");
271   fHistSPDmult->GetYaxis()->SetTitle("Entries");
272   fOutput->Add(fHistSPDmult);
273
274   fHistSPDmultcl1 = new TH1F("fHistSPDmultcl1", "Inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
275   fHistSPDmultcl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
276   fHistSPDmultcl1->GetYaxis()->SetTitle("Entries");
277   fOutput->Add(fHistSPDmultcl1);
278
279   fHistSPDmultcl2 = new TH1F("fHistSPDmultcl2", "Outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
280   fHistSPDmultcl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
281   fHistSPDmultcl2->GetYaxis()->SetTitle("Entries");
282   fOutput->Add(fHistSPDmultcl2);
283
284   fHistSPDmultcl1vscl2 = new TH2F("fHistSPDmultcl1vscl2", "Inner layer cluster vs outer multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
285   fHistSPDmultcl1vscl2->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
286   fHistSPDmultcl1vscl2->GetYaxis()->SetTitle("Outer layer cluster multiplicity");
287   fOutput->Add(fHistSPDmultcl1vscl2);
288
289   fHistSPDmultvscl1 = new TH2F("fHistSPDmultvscl1", "Tracklet vs inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
290   fHistSPDmultvscl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
291   fHistSPDmultvscl1->GetYaxis()->SetTitle("Tracklet multiplicity");
292   fOutput->Add(fHistSPDmultvscl1);
293
294   fHistSPDmultvscl2 = new TH2F("fHistSPDmultvscl2", "Tracklet vs outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
295   fHistSPDmultvscl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
296   fHistSPDmultvscl2->GetYaxis()->SetTitle("Tracklet multiplicity");
297   fOutput->Add(fHistSPDmultvscl2);
298
299   fHistSPDphi = new TH1F("fHistSPDphi", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
300   fHistSPDphi->GetXaxis()->SetTitle("#varphi [rad]");
301   fHistSPDphi->GetYaxis()->SetTitle("Entries");
302   fOutput->Add(fHistSPDphi);
303
304   fHistSPDtheta = new TH1F("fHistSPDtheta", "Tracklet #theta  distribution", 180, 0.,TMath::Pi());
305   fHistSPDtheta->GetXaxis()->SetTitle("#theta [rad]");
306   fHistSPDtheta->GetYaxis()->SetTitle("Entries");
307   fOutput->Add(fHistSPDtheta);
308
309   fHistSPDdePhi= new TH1F("fHistSPDdePhi", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
310   fHistSPDdePhi->GetXaxis()->SetTitle("#Delta#varphi [rad]");
311   fHistSPDdePhi->GetYaxis()->SetTitle("Entries");
312   fOutput->Add(fHistSPDdePhi);
313
314   fHistSPDphivsSPDeta= new TH2F("fHistSPDphivsSPDeta", "Tracklets - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());
315   fHistSPDphivsSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");
316   fHistSPDphivsSPDeta->GetYaxis()->SetTitle("#varphi [rad]");
317   fOutput->Add(fHistSPDphivsSPDeta);
318
319   fHistSPDdeTheta= new TH1F("fHistSPDdeTheta", "Tracklet #Delta#theta distribution",100,-0.05,.05);
320   fHistSPDdeTheta->GetXaxis()->SetTitle("#Delta#theta [rad]");
321   fHistSPDdeTheta->GetYaxis()->SetTitle("Entries");
322   fOutput->Add(fHistSPDdeTheta);
323
324   fHistSPDvtx = new TH1F("fHistSPDvtx", "SPD vertex distribution",nBinVtx,-MaxVtx,MaxVtx);
325   fHistSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
326   fHistSPDvtx->GetYaxis()->SetTitle("Entries");
327   fOutput->Add(fHistSPDvtx);
328
329   fHistSPDvtxAnalysis = new TH1F("fHistSPDvtxAnalysis", "SPD vertex distribution: events selected for the analysis",nBinVtx,-MaxVtx,MaxVtx);
330   fHistSPDvtxAnalysis->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
331   fHistSPDvtxAnalysis->GetYaxis()->SetTitle("Entries");
332   fOutput->Add(fHistSPDvtxAnalysis);
333
334   fHistSPDdePhideTheta= new TH2F("fHistSPDdePhideTheta", "Tracklet #Delta#varphi  distribution",2000,-1.,1.,1000,-0.25,.25);
335   fHistSPDdePhideTheta->GetXaxis()->SetTitle("#Delta#varphi [rad]");
336   fHistSPDdePhideTheta->GetYaxis()->SetTitle("#Delta#theta [rad]");
337   fOutput->Add(fHistSPDdePhideTheta);
338
339   fHistSPDphicl1 = new TH1F("fHistSPDphicl1", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
340   fHistSPDphicl1->GetXaxis()->SetTitle("#varphi [rad]");
341   fHistSPDphicl1->GetYaxis()->SetTitle("Entries");
342   fOutput->Add(fHistSPDphicl1);
343
344   fHistSPDphicl2 = new TH1F("fHistSPDphicl2", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
345   fHistSPDphicl2->GetXaxis()->SetTitle("#varphi [rad]");
346   fHistSPDphicl2->GetYaxis()->SetTitle("Entries");
347   fOutput->Add(fHistSPDphicl2);
348
349
350   if (fUseMC) {
351
352     // Track level correction histograms 
353     fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
354     fOutput->Add(fHistBkgCorrDen);
355    
356     fHistReconstructedProtons = new TH2F("fHistReconstructedProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
357     fOutput->Add(fHistReconstructedProtons);
358     fHistReconstructedKaons = new TH2F("fHistReconstructedKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
359     fOutput->Add(fHistReconstructedKaons);
360     fHistReconstructedPions = new TH2F("fHistReconstructedPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
361     fOutput->Add(fHistReconstructedPions);
362     fHistReconstructedOthers = new TH2F("fHistReconstructedOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
363     fOutput->Add(fHistReconstructedOthers);
364  
365     fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
366     fOutput->Add(fHistBkgCorrDenPrimGen);
367
368     fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
369     fOutput->Add(fHistBkgCombLabels);
370
371     if (fTR) {
372       fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
373       fOutput->Add(fHistBkgCorrNum);
374  
375       fHistAlgEffNum = new TH2F("fHistAlgEffNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
376       fOutput->Add(fHistAlgEffNum);  
377
378       fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
379       fOutput->Add(fHistNonDetectableCorrDen);
380
381     }
382
383     fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
384     fOutput->Add(fHistNonDetectableCorrNum);
385
386     fHistProtons = new TH2F("fHistProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
387     fOutput->Add(fHistProtons);
388     fHistKaons = new TH2F("fHistKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
389     fOutput->Add(fHistKaons);
390     fHistPions = new TH2F("fHistPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
391     fOutput->Add(fHistPions);
392     fHistOthers = new TH2F("fHistOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
393     fOutput->Add(fHistOthers);
394     fHistReconstructedSec = new TH2F("fHistReconstructedSec","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
395     fOutput->Add(fHistReconstructedSec); 
396
397  
398     fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
399     fOutput->Add(fHistAllPrimaries);
400
401     fHistTrackCentrEvts = new TH2F("fHistTrackCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
402     fOutput->Add(fHistTrackCentrEvts);
403
404     fHistTrackTrigCentrEvts = new TH2F("fHistTrackTrigCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
405     fOutput->Add(fHistTrackTrigCentrEvts);
406
407
408     // Event level correction histograms  
409     fHistAllEvts = new TH2F("fHistAllEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
410     fOutput->Add(fHistAllEvts);
411
412     fHistCentrEvts = new TH2F("fHistCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
413     fOutput->Add(fHistCentrEvts);
414
415     fHistTrigCentrEvts = new TH2F("fHistTrigCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
416     fOutput->Add(fHistTrigCentrEvts);
417
418     fHistSelEvts = new TH2F("fHistSelEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
419     fOutput->Add(fHistSelEvts);
420
421
422     // MC distributions
423     fHistMCmultEtacut = new TH1F("fHistMCmultEtacut","Generated multiplicity",nBinMultCorr,binLimMultCorr);
424     fHistMCmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
425     fHistMCmultEtacut->GetYaxis()->SetTitle("Entries");
426     fOutput->Add(fHistMCmultEtacut);
427
428     fHistMCmultEtacutvsSPDmultEtacut = new TH2F("fHistMCmultEtacutvsSPDmultEtacut","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
429     fHistMCmultEtacutvsSPDmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
430     fHistMCmultEtacutvsSPDmultEtacut->GetYaxis()->SetTitle("Tracklet multiplicity |#eta|<1.4");
431     fOutput->Add(fHistMCmultEtacutvsSPDmultEtacut);
432
433     fHistMCmultEtacutvsSPDmultcl1 = new TH2F("fHistMCmultEtacutvsSPDmultcl1","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
434     fHistMCmultEtacutvsSPDmultcl1->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
435     fHistMCmultEtacutvsSPDmultcl1->GetYaxis()->SetTitle("Cluster inner layer multiplicity");
436     fOutput->Add(fHistMCmultEtacutvsSPDmultcl1);
437
438     fHistMCmultEtacutvsSPDmultcl2 = new TH2F("fHistMCmultEtacutvsSPDmultcl2","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
439     fHistMCmultEtacutvsSPDmultcl2->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
440     fHistMCmultEtacutvsSPDmultcl2->GetYaxis()->SetTitle("Cluster outer layer multiplicity");
441     fOutput->Add(fHistMCmultEtacutvsSPDmultcl2);
442
443
444     fHistMCvtxx = new TH1F("fHistMCvtxx", "MC vertex distribution - x",100,-.5,.5);
445     fOutput->Add(fHistMCvtxx);
446     fHistMCvtxy = new TH1F("fHistMCvtxy", "MC vertex distribution - y",100,-.5,.5);
447     fOutput->Add(fHistMCvtxy);
448     fHistMCvtxz = new TH1F("fHistMCvtxz", "MC vertex distribution - z",500,-50.,50.);
449     fOutput->Add(fHistMCvtxz);
450
451     fHistRecvsGenImpactPar= new TH2F("fHistRecvsGenImpactPar","",20,0.,20.,20,0.,20.);
452     fHistRecvsGenImpactPar->GetXaxis()->SetTitle("b_{gen} [fm]");
453     fHistRecvsGenImpactPar->GetYaxis()->SetTitle("b_{rec} [fm]");
454     fOutput->Add(fHistRecvsGenImpactPar);
455
456     fHistMCNpart = new TH1F("fHistMCNpart","Number of participants",450,0,450);
457     fHistMCNpart->GetXaxis()->SetTitle("N_{part}");
458     fHistMCNpart->GetYaxis()->SetTitle("Entries"); 
459     fOutput->Add(fHistMCNpart);
460  
461     fHistdPhidThetaPP = new TH2F("fHistdPhidThetaPP","",2000,-1.,1.,1000,-0.25,.25);
462     fOutput->Add(fHistdPhidThetaPP);
463     fHistdPhidThetaSS = new TH2F("fHistdPhidThetaSS","",2000,-1.,1.,1000,-0.25,.25);
464     fOutput->Add(fHistdPhidThetaSS);
465     fHistdPhidThetaComb = new TH2F("fHistdPhidThetaComb","",2000,-1.,1.,1000,-0.25,.25);
466     fOutput->Add(fHistdPhidThetaComb);
467
468     fHistDeVtx = new TH2F("fHistDeVtx","",80,-20.,20.,5000,-0.5,0.5);
469     fHistDeVtx->GetXaxis()->SetTitle("z_{MCvtx} [cm]");
470     fHistDeVtx->GetYaxis()->SetTitle("z_{MCvtx}-z_{SPDvtx} [cm]");
471     fOutput->Add(fHistDeVtx);
472
473   }
474   PostData(1, fOutput);
475 //  Printf("Output objects created...");
476 }
477
478 //________________________________________________________________________
479 void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *) 
480 {
481   // Main loop
482
483   // Called for each event
484 //  Printf("User exec..........");
485
486   AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
487
488   fmyESD = dynamic_cast<AliESDEvent*>(InputEvent()); 
489   if (!fmyESD) {
490     Printf("ERROR: fESD not available");
491     return;
492   }
493
494     
495   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
496   if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
497
498
499   // Trigger selection
500   Bool_t eventTriggered = kTRUE;
501   static AliTriggerAnalysis* triggerAnalysis = 0; 
502   if (eventTriggered) // applying an offline trigger
503     eventTriggered = triggerAnalysis->IsTriggerFired(fmyESD, fTrigger);
504   if (!eventTriggered)
505     Printf("No trigger");
506
507
508   // Centrality selection 
509   Bool_t eventInCentralityBin = kFALSE;
510 /*  AliESDCentrality *centrality = fmyESD->GetCentrality();
511   if (fCentrEst=="") eventInCentralityBin = kTRUE;
512   else {
513     if(!centrality) {
514       AliError("Centrality object not available"); 
515     }  else {
516       if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
517     }
518   }
519 */
520   if (fCentrEst) {
521     AliESDVZERO* esdV0 = fmyESD->GetVZEROData();
522     Float_t multV0A=esdV0->GetMTotV0A();
523     Float_t multV0C=esdV0->GetMTotV0C();
524     fV0Ampl->Fill(multV0A+multV0C);
525     if (multV0A+multV0C>=fMinMultV0) eventInCentralityBin = kTRUE; 
526   } else if (!fCentrEst) {
527     eventInCentralityBin = kTRUE;
528   }
529
530   const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
531
532   // ESD vertex
533   Bool_t eventWithVertex = kFALSE;
534   const AliESDVertex* vtxESD = fmyESD->GetVertex();
535   const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
536   Double_t esdvtx[3];
537   vtxESD->GetXYZ(esdvtx);
538   Int_t nContrib = vtxESD->GetNContributors();
539   Int_t nContribTPC = vtxTPC->GetNContributors(); 
540   if (nContrib>0&&nContribTPC>0) {
541     if (vtxESD->GetDispersion()<0.04) 
542       if (vtxESD->GetZRes()<0.25) 
543         if (nContribTPC>(-10.+0.25*multESD->GetNumberOfITSClusters(0))) {
544           fHistSPDvtx->Fill(esdvtx[2]);
545           if (TMath::Abs(esdvtx[2])<fVtxLim)  
546             eventWithVertex = kTRUE;   
547         }
548   } 
549
550   // Reconstructing or loading tracklets...
551   Int_t multSPD = multESD->GetNumberOfTracklets();
552   Int_t nSingleCl1 = multESD->GetNumberOfSingleClusters();
553   Int_t multSPDcl1 = nSingleCl1 + multSPD;  
554   Int_t multSPDEtacut = 0;
555   Int_t multSPDAna = 0;
556
557   Int_t multSPDcl2 = multESD->GetNumberOfITSClusters(1);
558
559   Float_t trackletCoord[multSPDcl1][5];
560   Float_t trackletLab[multSPDcl1][2];
561
562   Bool_t eventSelected = kFALSE;
563   // Event selection: in centrality bin, triggered with vertex
564   if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
565     eventSelected = kTRUE; 
566     fHistSPDmultcl1->Fill(multSPDcl1);
567     fHistSPDmultcl2->Fill(multSPDcl2);
568     fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2);
569     fHistSPDmultvscl1->Fill(multSPDcl1,multSPD);
570     fHistSPDmultvscl2->Fill(multSPDcl2,multSPD);
571     if (fRecoTracklets) {
572       // Load ITS rec points tree
573       TTree *itsClusterTree = hand->GetTreeR("ITS");
574       if (!itsClusterTree) {
575         AliError(" Invalid ITS cluster tree !\n");
576         return;
577       }
578       float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()};
579
580       fMultReco->SetPhiWindow(fPhiWindow);
581       fMultReco->SetThetaWindow(fThetaWindow); 
582       fMultReco->SetPhiShift(fPhiShift); 
583       fMultReco->SetRemoveClustersFromOverlaps(fRemoveClustersFromOverlaps);
584       fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
585       fMultReco->SetZetaOverlapCut(fZetaOverlapCut); 
586       fMultReco->SetPhiRotationAngle(fPhiRotationAngle);
587       fMultReco->SetHistOn(kFALSE); 
588
589       fMultReco->Reconstruct(itsClusterTree,vtxf,vtxf);
590 //      Printf("cl 1 from alg %d",fMultReco->GetNClustersLayer1());
591 //      Printf("cl 2 from alg %d",fMultReco->GetNClustersLayer2());
592       multSPD = fMultReco->GetNTracklets(); 
593 //      Printf("tracklets found...%d",multSPD);
594       nSingleCl1 = fMultReco->GetNSingleClusters();
595 //      Printf("singl found...%d",nSingleCl1);
596       for (Int_t itr = 0; itr<multSPD;++itr) {
597         trackletCoord[itr][3] = fMultReco->GetTracklet(itr)[0]; //theta
598         trackletCoord[itr][2] = fMultReco->GetTracklet(itr)[1]; //phi
599         trackletCoord[itr][1] = fMultReco->GetTracklet(itr)[3]; //deTheta
600         trackletCoord[itr][0] = fMultReco->GetTracklet(itr)[2]; //dephi
601         trackletCoord[itr][4] = -TMath::Log(TMath::Tan(trackletCoord[itr][3]/2.)); //eta
602         trackletLab[itr][0] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[4]);  //label lay1
603         trackletLab[itr][1] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[5]);  //label lay2
604       }
605       for (Int_t icl1 = 0; icl1<multSPDcl1;++icl1) {
606         fHistSPDphicl1->Fill(fMultReco->GetClusterLayer1(icl1)[1]); 
607       }
608       for (Int_t icl2 = 0; icl2<multSPDcl2; ++icl2) {
609         fHistSPDphicl2->Fill(fMultReco->GetClusterLayer2(icl2)[1]);          
610       }
611
612     } else {
613       // set variables from ESD
614       for (Int_t itr = 0; itr<multSPD;++itr) {
615         trackletCoord[itr][3] = multESD->GetTheta(itr);       //theta
616         trackletCoord[itr][2] = multESD->GetPhi(itr);         //phi
617         trackletCoord[itr][1] = multESD->GetDeltaTheta(itr);  //deTheta
618         trackletCoord[itr][0] = multESD->GetDeltaPhi(itr);    //dePhi
619         trackletCoord[itr][4] = multESD->GetEta(itr);         //eta 
620         trackletLab[itr][0] = multESD->GetLabel(itr,0);       //label lay1
621         trackletLab[itr][1] = multESD->GetLabel(itr,1);       //label lay2
622       }
623     }
624 //    Printf("tracklets in ESD...%d",multESD->GetNumberOfTracklets());
625
626     for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
627       if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) { // select tracklets with tighter cuts
628         fHistSPDRAWEtavsZ->Fill(trackletCoord[itracklet][4],esdvtx[2]);
629         fHistSPDphi->Fill(trackletCoord[itracklet][2]);
630         fHistSPDtheta->Fill(trackletCoord[itracklet][3]);
631         fHistSPDdePhi->Fill(trackletCoord[itracklet][0]);
632         fHistSPDdeTheta->Fill(trackletCoord[itracklet][1]);
633  
634         fHistSPDphivsSPDeta->Fill(trackletCoord[itracklet][4],trackletCoord[itracklet][2]);
635         multSPDAna++;
636         // Calculate multiplicity in etacut
637         if (TMath::Abs(trackletCoord[itracklet][4])<1.4) multSPDEtacut++;
638       }
639       // for combinatorial background normalization
640       fHistSPDdePhideTheta->Fill(trackletCoord[itracklet][0], trackletCoord[itracklet][1]);
641
642     }
643     if (multSPDAna!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
644     fHistSPDmultEtacut->Fill(multSPDEtacut);
645     fHistSPDmult->Fill(multSPDAna); 
646
647     fHistSPDRAWMultvsZ->Fill(multSPDAna,esdvtx[2]);
648
649   } // End selected events
650
651   if (eventInCentralityBin) {
652     fHistSPDRAWMultvsZCentrEvts->Fill(multSPDAna,esdvtx[2]);  
653     if (eventTriggered) fHistSPDRAWMultvsZTriggCentrEvts->Fill(multSPDAna,esdvtx[2]);
654   }
655
656   if (fUseMC) {
657
658     if (!fMCEvent) {
659       AliError("No MC info found");
660       return;
661     } 
662     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
663     if (!eventHandler) {
664       Printf("ERROR: Could not retrieve MC event handler");
665       return;
666     }
667
668     AliStack* stack = fMCEvent->Stack();
669     if (!stack) {
670       AliDebug(AliLog::kError, "Stack not available");
671       return;
672     }
673
674     Float_t impactParameter = 0.;
675     Int_t npart = 0;
676
677     AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
678     AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
679
680     if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
681     else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
682
683     Bool_t IsEventInMCCentralityBin = kFALSE;
684     switch (fMCCentralityBin) {
685
686       case 1: if (impactParameter<3) IsEventInMCCentralityBin = kTRUE;
687       break;
688       case 2: if (impactParameter>3&&impactParameter<6) IsEventInMCCentralityBin = kTRUE;
689       break;
690       case 3: if (impactParameter>6&&impactParameter<9) IsEventInMCCentralityBin = kTRUE;
691       break;
692       case 4: if (impactParameter>9&&impactParameter<12) IsEventInMCCentralityBin = kTRUE;
693       break;
694       case 5: if (impactParameter>12&&impactParameter<15) IsEventInMCCentralityBin = kTRUE;
695       break;
696       case 6: if (impactParameter>15) IsEventInMCCentralityBin = kTRUE;
697       break;
698       case 7:  IsEventInMCCentralityBin = kTRUE;
699       break;
700
701     }
702
703     if (IsEventInMCCentralityBin) {
704       // MC vertex
705       TArrayF vtxMC(3);
706       fMCEvent->GenEventHeader()->PrimaryVertex(vtxMC);
707       fHistMCvtxx->Fill(vtxMC[0]);   
708       fHistMCvtxy->Fill(vtxMC[1]);
709       fHistMCvtxz->Fill(vtxMC[2]);
710
711 //      Printf("Impact parameter gen: %f", impactParameter);
712        if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
713        else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
714
715       //Rec centrality vs gen centrality
716       AliESDZDC *zdcRec = fmyESD->GetESDZDC();
717       Double_t impactParameterZDC = zdcRec->GetImpactParameter();
718 //      Printf("Impact parameter rec: %f", impactParameterZDC);
719
720       fHistRecvsGenImpactPar->Fill(impactParameter,impactParameterZDC);
721       fHistMCNpart->Fill(npart);
722   
723       // Tracks from MC
724       Int_t  multMCCharged = 0;
725       Int_t  multMCChargedEtacut = 0;
726 //      Int_t  nMCPart = stack->GetNprimary();
727       Int_t  nMCPart = stack->GetNtrack();  // decay products of D and B mesons are also primaries and produced in HIJING during transport
728       Float_t* etagen = new Float_t[nMCPart];  
729       Int_t* stackIndexOfPrimaryParts = new Int_t[nMCPart];
730       Bool_t* reconstructedPrimaryPart = new Bool_t[nMCPart];
731       Bool_t* detectedPrimaryPartLay1 = new Bool_t[nMCPart];
732       Bool_t* detectedPrimaryPartLay2 = new Bool_t[nMCPart];
733       Bool_t* detectablePrimaryPart = new Bool_t[nMCPart];
734
735       TTree* tRefTree;  
736       if (fTR) {
737         tRefTree = eventHandler->TreeTR(); 
738         fMCEvent->ConnectTreeTR(tRefTree); 
739       }
740
741       // Loop over MC particles
742       for (Int_t imc=0; imc<nMCPart; imc++) {
743         AliMCParticle *mcpart  = (AliMCParticle*)fMCEvent->GetTrack(imc);
744
745         Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
746         if (!isPrimary)            continue;
747         if (mcpart->Charge() == 0) continue;
748         Float_t theta = mcpart->Theta();
749         if (theta==0 || theta==TMath::Pi())    continue;
750         Float_t eta = mcpart->Eta();
751 //        Float_t pt = mcpart->Pt();
752         etagen[multMCCharged] = eta; 
753         stackIndexOfPrimaryParts[multMCCharged] = imc;
754
755         reconstructedPrimaryPart[multMCCharged]=kFALSE;
756         detectedPrimaryPartLay1[multMCCharged]=kFALSE;
757         detectedPrimaryPartLay2[multMCCharged]=kFALSE;
758         detectablePrimaryPart[multMCCharged]=kFALSE;
759
760         if (fTR) {
761           Int_t nref = mcpart->GetNumberOfTrackReferences();
762           if (nref==0) {
763             detectablePrimaryPart[multMCCharged]= kTRUE;
764           } else if (nref>0) {
765             // Check if the primary is detectable 
766             detectablePrimaryPart[multMCCharged] = IsDetectablePrimary(nref,mcpart);
767             // Check if the primary is detected (if SPD 100% efficient) 
768             if (detectablePrimaryPart[multMCCharged]) {
769               detectedPrimaryPartLay1[multMCCharged] = IsDetectedPrimary(nref,mcpart,0);
770               detectedPrimaryPartLay2[multMCCharged] = IsDetectedPrimary(nref,mcpart,1);
771             }
772           }  
773         }
774         if (eventSelected&&fPartSpecies) {
775           if (TMath::Abs(mcpart->PdgCode())==2212) fHistProtons->Fill(etagen[multMCCharged],vtxMC[2]);
776           else if (TMath::Abs(mcpart->PdgCode())==321) fHistKaons->Fill(etagen[multMCCharged],vtxMC[2]);
777           else if (TMath::Abs(mcpart->PdgCode())==211) fHistPions->Fill(etagen[multMCCharged],vtxMC[2]);
778           else fHistOthers->Fill(etagen[multMCCharged],vtxMC[2]);  //includes leptons pdg->11,13
779         }
780         multMCCharged++;
781         if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
782       } // end of MC particle loop
783
784       fHistMCmultEtacut->Fill(multMCChargedEtacut);
785
786       // Event selection: in centrality bin, triggered with vertex 
787       if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
788     
789         fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);      
790
791         for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
792           if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) 
793             fHistBkgCorrDen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
794
795           if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) {
796             Bool_t trakletByPrim = kFALSE; 
797             for (Int_t imc=0; imc<multMCCharged; imc++) {
798               if (trackletLab[itracklet][0]==stackIndexOfPrimaryParts[imc]) {
799                 fHistdPhidThetaPP->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
800                 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
801                   fHistBkgCorrDenPrimGen->Fill(etagen[imc],vtxMC[2]);
802                   if (fPartSpecies) {
803                     if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==2212) 
804                       fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
805                     else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==321) 
806                       fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
807                     else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==211) 
808                       fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
809                     else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
810                   }     
811                 }
812                 trakletByPrim = kTRUE; 
813                 if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
814                   reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers 
815                 }
816                 break;
817               }  
818             }
819             if (!trakletByPrim) {
820               fHistdPhidThetaSS->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
821               if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
822
823                 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
824                 if (fPartSpecies) {
825                   Int_t motherlab = ((AliMCParticle*)fMCEvent->GetTrack((Int_t)trackletLab[itracklet][0]))->GetMother(); 
826                   if (motherlab>-1) { 
827                     if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==2212) fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
828                     else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==321) fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
829                     else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==211) fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
830                     else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
831                   } else fHistReconstructedSec->Fill(trackletCoord[itracklet][4],esdvtx[2]);
832                 }
833               }
834             }
835           } else {
836             fHistdPhidThetaComb->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
837             if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
838               fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
839               fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]); 
840             }
841           }
842         } // end loop tracklets
843
844         for (Int_t imc=0; imc<multMCCharged; imc++) {
845           if (fTR) {
846             if (reconstructedPrimaryPart[imc]) fHistBkgCorrNum->Fill(etagen[imc],vtxMC[2]); // only for separate corrections
847             if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) fHistAlgEffNum->Fill(etagen[imc],vtxMC[2]);
848             if (detectablePrimaryPart[imc]) fHistNonDetectableCorrDen->Fill(etagen[imc],vtxMC[2]); 
849           }
850           fHistNonDetectableCorrNum->Fill(etagen[imc],vtxMC[2]);
851         }
852
853         fHistSelEvts->Fill(multSPDAna,vtxMC[2]);
854       } // end of selected events
855
856       fHistMCmultEtacutvsSPDmultEtacut->Fill(multMCChargedEtacut,multSPDEtacut);
857       fHistMCmultEtacutvsSPDmultcl1->Fill(multMCChargedEtacut,multSPDcl1);
858       fHistMCmultEtacutvsSPDmultcl2->Fill(multMCChargedEtacut,multSPDcl2);
859       fHistAllEvts->Fill(multSPDAna,vtxMC[2]);
860
861       if (eventInCentralityBin) {
862         fHistCentrEvts->Fill(multSPDAna,vtxMC[2]);
863         if (eventTriggered) {
864           fHistTrigCentrEvts->Fill(multSPDAna,vtxMC[2]);
865         }
866       }
867
868       for (Int_t imc=0; imc<multMCCharged; imc++) { 
869         fHistAllPrimaries->Fill(etagen[imc],vtxMC[2]);
870         if (eventInCentralityBin) {
871           fHistTrackCentrEvts->Fill(etagen[imc],vtxMC[2]);
872           if (eventTriggered) fHistTrackTrigCentrEvts->Fill(etagen[imc],vtxMC[2]); 
873         }
874       } 
875     
876       delete[] etagen;
877       delete[] stackIndexOfPrimaryParts;
878       delete[] reconstructedPrimaryPart;
879       delete[] detectedPrimaryPartLay1;
880       delete[] detectedPrimaryPartLay2;
881       delete[] detectablePrimaryPart;
882     }  
883   }
884   PostData(1, fOutput);
885 }      
886 //________________________________________________________________________
887 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectedPrimary(Int_t nref, AliMCParticle* mcpart, Int_t Layer) {
888   Double_t rMinLay1 = 3.4; //min rad for track ref first SPD layer 
889   Double_t rMaxLay1 = 4.4; //max rad for track ref first SPD layer 
890   Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer  
891   Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer 
892   Double_t rmin = (Layer > 0 ? rMinLay2 : rMinLay1);
893   Double_t rmax = (Layer > 0 ? rMaxLay2 : rMaxLay1);
894   AliTrackReference *tref=0x0;
895
896   for (Int_t iref=0;iref<nref;iref++) {
897      tref = mcpart->GetTrackReference(iref);
898    if (tref) {
899        if (tref->R()>rmin&&tref->R()<rmax) {
900          if (tref->DetectorId()==0) {
901            return kTRUE;
902          }
903        }
904      }
905   }
906   return kFALSE;
907 }
908 //________________________________________________________________________
909 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectablePrimary(Int_t nref, AliMCParticle* mcpart) {
910
911   Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer  
912   Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer 
913
914   AliTrackReference *tref= mcpart->GetTrackReference(nref-1);
915   if (tref->DetectorId()!=-1) {
916     return kTRUE;
917   } else { //last is -1 -> particle disappeared. Where?
918     if (tref->R()>rMaxLay2) {
919       return kTRUE;  
920     } else if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2) { // this last tr is in lay 2
921       for (Int_t iref=0;iref<nref;iref++) {  // look for other tr in lay 2
922         tref = mcpart->GetTrackReference(iref);
923         if (tref) 
924           if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2) 
925             if (tref->DetectorId()==0) return kTRUE; 
926       } 
927     } // else particle disappeared before lay2
928   }
929   return kFALSE; 
930 }
931 //________________________________________________________________________
932 void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *) 
933 {
934   Printf("Terminating...");
935 //  AliAnalysisTaskSE::Terminate();
936
937 }