1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //////////////////////////////////////////////////////////////////////////
17 // Class AliAnalysisTaskSPDdNdEta //
18 // Analysis task for dN/dEta reconstruction with the SPD //
20 // Author: M. Nicassio (INFN Bari) //
21 // Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it //
22 //////////////////////////////////////////////////////////////////////////
31 #include "TGeoGlobalMagField.h"
33 #include "AliAnalysisManager.h"
35 #include "AliMultiplicity.h"
36 #include "AliESDEvent.h"
37 #include "AliESDInputHandler.h"
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"
48 #include "AliMCEventHandler.h"
49 #include "AliMCEvent.h"
51 #include "AliTrackReference.h"
53 #include "AliGenHijingEventHeader.h"
54 #include "AliGenDPMjetEventHeader.h"
58 #include "AliTriggerAnalysis.h"
59 #include "AliPhysicsSelection.h"
60 #include "AliTrackletAlg.h"
61 #include "AliAnalysisTaskSPDdNdEta.h"
64 ClassImp(AliAnalysisTaskSPDdNdEta)
65 //________________________________________________________________________
66 AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
67 : AliAnalysisTaskSE(name),
73 fTrigger(AliTriggerAnalysis::kAcceptAll),
75 fRecoTracklets(kFALSE),
77 fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall),
90 fRemoveClustersFromOverlaps(0),
99 fHistSPDRAWMultvsZ(0),
100 fHistSPDRAWMultvsZTriggCentrEvts(0),
101 fHistSPDRAWMultvsZCentrEvts(0),
102 fHistSPDRAWEtavsZ(0),
104 fHistSPDmultEtacut(0),
108 fHistSPDmultcl1vscl2(0),
109 fHistSPDmultvscl1(0),
110 fHistSPDmultvscl2(0),
116 fHistSPDphivsSPDeta(0),
119 fHistSPDvtxAnalysis(0),
120 fHistSPDdePhideTheta(0),
126 fHistReconstructedProtons(0),
127 fHistReconstructedKaons(0),
128 fHistReconstructedPions(0),
129 fHistReconstructedOthers(0),
130 fHistReconstructedSec(0),
131 fHistBkgCorrDenPrimGen(0),
132 fHistBkgCombLabels(0),
135 fHistNonDetectableCorrDen(0),
137 fHistNonDetectableCorrNum(0),
142 fHistAllPrimaries(0),
143 fHistTrackCentrEvts(0),
144 fHistTrackTrigCentrEvts(0),
148 fHistTrigCentrEvts(0),
151 fHistMCmultEtacut(0),
152 fHistMCmultEtacutvsSPDmultEtacut(0),
153 fHistMCmultEtacutvsSPDmultcl1(0),
154 fHistMCmultEtacutvsSPDmultcl2(0),
160 fHistRecvsGenImpactPar(0),
163 fHistdPhidThetaPP(0),
164 fHistdPhidThetaSS(0),
165 fHistdPhidThetaComb(0),
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());
180 //________________________________________________________________________
181 AliAnalysisTaskSPDdNdEta::~AliAnalysisTaskSPDdNdEta()
186 // histograms are in the output list and deleted when the output
187 // list is deleted by the TSelector dtor
190 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
196 //________________________________________________________________________
197 void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects()
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");
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");
212 fMultReco = new AliTrackletAlg();
215 fOutput = new TList();
219 Double_t MaxVtx = 20.;
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;
230 Int_t nBinEtaCorr = 60;
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;
238 // Data to be corrected
241 fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
242 fOutput->Add(fV0Ampl);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
353 // Track level correction histograms
354 fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
355 fOutput->Add(fHistBkgCorrDen);
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);
366 fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
367 fOutput->Add(fHistBkgCorrDenPrimGen);
369 fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
370 fOutput->Add(fHistBkgCombLabels);
373 fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
374 fOutput->Add(fHistBkgCorrNum);
376 fHistAlgEffNum = new TH2F("fHistAlgEffNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
377 fOutput->Add(fHistAlgEffNum);
379 fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
380 fOutput->Add(fHistNonDetectableCorrDen);
384 fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
385 fOutput->Add(fHistNonDetectableCorrNum);
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);
399 fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
400 fOutput->Add(fHistAllPrimaries);
402 fHistTrackCentrEvts = new TH2F("fHistTrackCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
403 fOutput->Add(fHistTrackCentrEvts);
405 fHistTrackTrigCentrEvts = new TH2F("fHistTrackTrigCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
406 fOutput->Add(fHistTrackTrigCentrEvts);
409 // Event level correction histograms
410 fHistAllEvts = new TH2F("fHistAllEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
411 fOutput->Add(fHistAllEvts);
413 fHistCentrEvts = new TH2F("fHistCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
414 fOutput->Add(fHistCentrEvts);
416 fHistTrigCentrEvts = new TH2F("fHistTrigCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
417 fOutput->Add(fHistTrigCentrEvts);
419 fHistSelEvts = new TH2F("fHistSelEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
420 fOutput->Add(fHistSelEvts);
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);
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);
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);
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);
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);
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);
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);
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);
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);
475 PostData(1, fOutput);
476 // Printf("Output objects created...");
479 //________________________________________________________________________
480 void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *)
484 // Called for each event
485 // Printf("User exec..........");
487 AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
488 if (!hand) { printf("No RP handler\n"); return; }
490 fmyESD = dynamic_cast<AliESDEvent*>(InputEvent());
492 Printf("ERROR: fESD not available");
496 AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
497 if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
501 Bool_t eventTriggered = kTRUE;
502 static AliTriggerAnalysis* triggerAnalysis = 0;
503 if (eventTriggered) // applying an offline trigger
504 eventTriggered = triggerAnalysis->IsTriggerFired(fmyESD, fTrigger);
506 Printf("No trigger");
509 // Centrality selection
510 Bool_t eventInCentralityBin = kFALSE;
511 /* AliESDCentrality *centrality = fmyESD->GetCentrality();
512 if (fCentrEst=="") eventInCentralityBin = kTRUE;
515 AliError("Centrality object not available");
517 if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
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;
531 const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
534 Bool_t eventWithVertex = kFALSE;
535 const AliESDVertex* vtxESD = fmyESD->GetVertex();
536 const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
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;
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;
558 Int_t multSPDcl2 = multESD->GetNumberOfITSClusters(1);
560 Float_t trackletCoord[multSPDcl1][5];
561 Float_t trackletLab[multSPDcl1][2];
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");
579 float vtxf[3] = {static_cast<float>(vtxESD->GetX()),static_cast<float>(vtxESD->GetY()),static_cast<float>(vtxESD->GetZ())};
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);
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
606 for (Int_t icl1 = 0; icl1<multSPDcl1;++icl1) {
607 fHistSPDphicl1->Fill(fMultReco->GetClusterLayer1(icl1)[1]);
609 for (Int_t icl2 = 0; icl2<multSPDcl2; ++icl2) {
610 fHistSPDphicl2->Fill(fMultReco->GetClusterLayer2(icl2)[1]);
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
625 // Printf("tracklets in ESD...%d",multESD->GetNumberOfTracklets());
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]);
635 fHistSPDphivsSPDeta->Fill(trackletCoord[itracklet][4],trackletCoord[itracklet][2]);
637 // Calculate multiplicity in etacut
638 if (TMath::Abs(trackletCoord[itracklet][4])<1.4) multSPDEtacut++;
640 // for combinatorial background normalization
641 fHistSPDdePhideTheta->Fill(trackletCoord[itracklet][0], trackletCoord[itracklet][1]);
644 if (multSPDAna!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
645 fHistSPDmultEtacut->Fill(multSPDEtacut);
646 fHistSPDmult->Fill(multSPDAna);
648 fHistSPDRAWMultvsZ->Fill(multSPDAna,esdvtx[2]);
650 } // End selected events
652 if (eventInCentralityBin) {
653 fHistSPDRAWMultvsZCentrEvts->Fill(multSPDAna,esdvtx[2]);
654 if (eventTriggered) fHistSPDRAWMultvsZTriggCentrEvts->Fill(multSPDAna,esdvtx[2]);
660 AliError("No MC info found");
663 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
665 Printf("ERROR: Could not retrieve MC event handler");
669 AliStack* stack = fMCEvent->Stack();
671 AliDebug(AliLog::kError, "Stack not available");
675 Float_t impactParameter = 0.;
678 AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
679 AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
681 if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
682 else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
684 Bool_t IsEventInMCCentralityBin = kFALSE;
685 switch (fMCCentralityBin) {
687 case 1: if (impactParameter<3) IsEventInMCCentralityBin = kTRUE;
689 case 2: if (impactParameter>3&&impactParameter<6) IsEventInMCCentralityBin = kTRUE;
691 case 3: if (impactParameter>6&&impactParameter<9) IsEventInMCCentralityBin = kTRUE;
693 case 4: if (impactParameter>9&&impactParameter<12) IsEventInMCCentralityBin = kTRUE;
695 case 5: if (impactParameter>12&&impactParameter<15) IsEventInMCCentralityBin = kTRUE;
697 case 6: if (impactParameter>15) IsEventInMCCentralityBin = kTRUE;
699 case 7: IsEventInMCCentralityBin = kTRUE;
704 if (IsEventInMCCentralityBin) {
707 fMCEvent->GenEventHeader()->PrimaryVertex(vtxMC);
708 fHistMCvtxx->Fill(vtxMC[0]);
709 fHistMCvtxy->Fill(vtxMC[1]);
710 fHistMCvtxz->Fill(vtxMC[2]);
712 // Printf("Impact parameter gen: %f", impactParameter);
713 if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
714 else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
716 //Rec centrality vs gen centrality
717 AliESDZDC *zdcRec = fmyESD->GetESDZDC();
718 Double_t impactParameterZDC = zdcRec->GetImpactParameter();
719 // Printf("Impact parameter rec: %f", impactParameterZDC);
721 fHistRecvsGenImpactPar->Fill(impactParameter,impactParameterZDC);
722 fHistMCNpart->Fill(npart);
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];
738 tRefTree = eventHandler->TreeTR();
739 fMCEvent->ConnectTreeTR(tRefTree);
742 // Loop over MC particles
743 for (Int_t imc=0; imc<nMCPart; imc++) {
744 AliMCParticle *mcpart = (AliMCParticle*)fMCEvent->GetTrack(imc);
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;
756 reconstructedPrimaryPart[multMCCharged]=kFALSE;
757 detectedPrimaryPartLay1[multMCCharged]=kFALSE;
758 detectedPrimaryPartLay2[multMCCharged]=kFALSE;
759 detectablePrimaryPart[multMCCharged]=kFALSE;
762 Int_t nref = mcpart->GetNumberOfTrackReferences();
764 detectablePrimaryPart[multMCCharged]= kTRUE;
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);
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
782 if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
783 } // end of MC particle loop
785 fHistMCmultEtacut->Fill(multMCChargedEtacut);
787 // Event selection: in centrality bin, triggered with vertex
788 if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
790 fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);
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]);
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]);
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]);
813 trakletByPrim = kTRUE;
814 if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
815 reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers
820 if (!trakletByPrim) {
821 fHistdPhidThetaSS->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
822 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
824 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
826 Int_t motherlab = ((AliMCParticle*)fMCEvent->GetTrack((Int_t)trackletLab[itracklet][0]))->GetMother();
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]);
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]);
843 } // end loop tracklets
845 for (Int_t imc=0; imc<multMCCharged; imc++) {
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]);
851 fHistNonDetectableCorrNum->Fill(etagen[imc],vtxMC[2]);
854 fHistSelEvts->Fill(multSPDAna,vtxMC[2]);
855 } // end of selected events
857 fHistMCmultEtacutvsSPDmultEtacut->Fill(multMCChargedEtacut,multSPDEtacut);
858 fHistMCmultEtacutvsSPDmultcl1->Fill(multMCChargedEtacut,multSPDcl1);
859 fHistMCmultEtacutvsSPDmultcl2->Fill(multMCChargedEtacut,multSPDcl2);
860 fHistAllEvts->Fill(multSPDAna,vtxMC[2]);
862 if (eventInCentralityBin) {
863 fHistCentrEvts->Fill(multSPDAna,vtxMC[2]);
864 if (eventTriggered) {
865 fHistTrigCentrEvts->Fill(multSPDAna,vtxMC[2]);
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]);
878 delete[] stackIndexOfPrimaryParts;
879 delete[] reconstructedPrimaryPart;
880 delete[] detectedPrimaryPartLay1;
881 delete[] detectedPrimaryPartLay2;
882 delete[] detectablePrimaryPart;
885 PostData(1, fOutput);
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;
897 for (Int_t iref=0;iref<nref;iref++) {
898 tref = mcpart->GetTrackReference(iref);
900 if (tref->R()>rmin&&tref->R()<rmax) {
901 if (tref->DetectorId()==0) {
909 //________________________________________________________________________
910 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectablePrimary(Int_t nref, AliMCParticle* mcpart) {
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
915 AliTrackReference *tref= mcpart->GetTrackReference(nref-1);
916 if (tref->DetectorId()!=-1) {
918 } else { //last is -1 -> particle disappeared. Where?
919 if (tref->R()>rMaxLay2) {
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);
925 if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2)
926 if (tref->DetectorId()==0) return kTRUE;
928 } // else particle disappeared before lay2
932 //________________________________________________________________________
933 void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
935 Printf("Terminating...");
936 // AliAnalysisTaskSE::Terminate();