]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/correlationHF/AliAnalysisTaskSEmcCorr.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSEmcCorr.cxx
CommitLineData
d33db6d9
AR		 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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//
bb236562 18// Class AliAnalysisTaskSEmcCorr
19// AliAnalysisTaskSE for studying HF-(hadron,electrons) and hadron-hadron correlations
20// at MC level
21// Author: Andrea Rossi, andrea.rossi@cern.ch
d33db6d9
AR		 22/////////////////////////////////////////////////////////////
23
24
25#include <TH1F.h>
26#include <TH2F.h>
27#include <TAxis.h>
28#include <THnSparse.h>
29#include <TDatabasePDG.h>
30#include <TMath.h>
31#include <TROOT.h>
32#include <TArrayI.h>
3333be85 33#include <TString.h>
d33db6d9
AR		 34#include "AliAODEvent.h"
35#include "AliAODRecoDecayHF2Prong.h"
36#include "AliAODRecoDecayHF.h"
37#include "AliAODRecoDecay.h"
38#include "AliAnalysisDataSlot.h"
39#include "AliAnalysisDataContainer.h"
40#include "AliAODTrack.h"
41#include "AliAODHandler.h"
42#include "AliESDtrack.h"
43#include "AliAODVertex.h"
44#include "AliESDVertex.h"
45#include "AliVertexerTracks.h"
46#include "AliAODMCParticle.h"
47#include "AliAODPid.h"
48#include "AliTPCPIDResponse.h"
49#include "AliAODMCHeader.h"
50#include "AliAnalysisVertexingHF.h"
51#include "AliRDHFCutsD0toKpi.h"
52#include "AliAODInputHandler.h"
53#include "AliAnalysisManager.h"
54#include "AliNormalizationCounter.h"
55//#include "/Users/administrator/ALICE/CHARM/HFCJ/TestsForProduction/AliAnalysisTaskSEmcCorr.h"
56#include "AliAnalysisTaskSEmcCorr.h"
3333be85 57#include "AliGenEventHeader.h"
d33db6d9
AR		 58
59class TCanvas;
60class TTree;
61class TChain;
62class AliAnalysisTaskSE;
63
64
65ClassImp(AliAnalysisTaskSEmcCorr)
66
67AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr()
68: AliAnalysisTaskSE(),
69 fReadMC(kTRUE),
3333be85 70 fFastSimul(kFALSE),
71 fCheckDecay(kTRUE),
72 fDoHFCorrelations(kTRUE),
d33db6d9
AR		 73 fOnlyKine(kTRUE),
74 fDoHadronHadron(kFALSE),
75 fNentries(0),
76 fhNprongsD0(0),
77 fhNprongsD0chargedOnly(0),
78 fhNprongsD0chargedRef(0),
79 fYrangeTrig(0.5),
80 fEtarangeEleTrig(0.8),
81 fminDpt(1.),
82 fArrayDaugh(0x0),
83 flastdaugh(0),
84 fArrayAssoc(0x0),
85 fLastAss(0),
86 fminPtass(0.3),
87 fmaxEtaAss(99.),
88 fhMCcorrel(0x0),
89 fhMCtrigPart(0x0),
90 fhMChadroncorrel(0x0),
3333be85 91 fhMChadrontrigPart(0x0),
92 fhDzeroDecay(0x0),
93 fhDplusDecay(0x0),
94 fhLambdaCDecay(0x0),
95 fhAllBDecay(0x0),
96 fnpionPlus(0),
97 fnpionMinus(0),
98 fnpionZero(0),
99 fnkaonPlus(0),
100 fnkaonMinus(0),
101 fnkaonZeroShort(0),
102 fnProton(0),
103 fnAntiProton(0),
104 fnElePlus(0),
105 fnEleMinus(0),
106 fnMuonPlus(0),
107 fnMuonMinus(0),
108 fnNeutrino(0),
109 fnJpsi(0),
110 fnPhysPrim(0),
111 fpxtotdaugh(0),
112 fpytotdaugh(0),
113 fpztotdaugh(0),
114 fGeneratorString(0)
d33db6d9
AR		 115{// standard constructor
116
117}
118
119
120
121//________________________________________________________________________
122AliAnalysisTaskSEmcCorr::AliAnalysisTaskSEmcCorr(const char *name)
123 : AliAnalysisTaskSE(name),
124 fReadMC(kTRUE),
3333be85 125 fFastSimul(kFALSE),
126 fCheckDecay(kTRUE),
127 fDoHFCorrelations(kTRUE),
d33db6d9
AR		 128 fOnlyKine(kTRUE),
129 fDoHadronHadron(kFALSE),
130 fNentries(0),
131 fhNprongsD0(0),
132 fhNprongsD0chargedOnly(0),
133 fhNprongsD0chargedRef(0),
134 fYrangeTrig(0.5),
135 fEtarangeEleTrig(0.8),
136 fminDpt(1.),
137 fArrayDaugh(0x0),
138 flastdaugh(0),
139 fArrayAssoc(0x0),
140 fLastAss(0),
141 fminPtass(0.3),
142 fmaxEtaAss(99.),
143 fhMCcorrel(0x0),
144 fhMCtrigPart(0x0),
145 fhMChadroncorrel(0x0),
3333be85 146 fhMChadrontrigPart(0x0),
147 fhDzeroDecay(0x0),
148 fhDplusDecay(0x0),
149 fhLambdaCDecay(0x0),
150 fhAllBDecay(0x0),
151 fnpionPlus(0),
152 fnpionMinus(0),
153 fnpionZero(0),
154 fnkaonPlus(0),
155 fnkaonMinus(0),
156 fnkaonZeroShort(0),
157 fnProton(0),
158 fnAntiProton(0),
159 fnElePlus(0),
160 fnEleMinus(0),
161 fnMuonPlus(0),
162 fnMuonMinus(0),
163 fnNeutrino(0),
164 fnJpsi(0),
165 fnPhysPrim(0),
166 fpxtotdaugh(0),
167 fpytotdaugh(0),
168 fpztotdaugh(0),
169 fGeneratorString(0)
d33db6d9
AR		 170{ // default constructor
171
172 DefineOutput(1, TH1F::Class());
173 DefineOutput(2, TH1F::Class());
174 DefineOutput(3, TH1F::Class());
175 DefineOutput(4, TH1F::Class());
176 DefineOutput(5, THnSparseF::Class());
177 DefineOutput(6, THnSparseF::Class());
178 DefineOutput(7, THnSparseF::Class());
179 DefineOutput(8, THnSparseF::Class());
3333be85 180 DefineOutput(9,TH1F::Class());
181 DefineOutput(10,TH1F::Class());
182 DefineOutput(11,TH1F::Class());
183 DefineOutput(12,TH1F::Class());
d33db6d9
AR		 184}
185
186//________________________________________________________________________
187AliAnalysisTaskSEmcCorr::~AliAnalysisTaskSEmcCorr(){
188 // destructor
189 delete fNentries;
190 delete fhNprongsD0;
191 delete fhNprongsD0chargedOnly;
192 delete fhNprongsD0chargedRef;
193 delete fhMCcorrel;
194 delete fhMCtrigPart;
195 delete fhMChadroncorrel;
196 delete fhMChadrontrigPart;
197 delete fArrayDaugh;
198 delete fArrayAssoc;
3333be85 199 delete fhDzeroDecay;
200 delete fhDplusDecay;
201 delete fhLambdaCDecay;
202 delete fhAllBDecay;
d33db6d9
AR		 203}
204
205//________________________________________________________________________
206void AliAnalysisTaskSEmcCorr::Init()
207{
208 // Initialization
209 fArrayDaugh=new TArrayI(20);
210
211 fArrayAssoc=new TArrayI(200);
212
213}
3333be85 214//_______________________________________________
215void AliAnalysisTaskSEmcCorr::GetFinalStateParticles(AliAODMCParticle *part,TClonesArray *clarr){
216
217 Int_t fdaugh=part->GetDaughter(0);
218 Int_t ldaugh=part->GetDaughter(1);
219 if(fdaugh<0||ldaugh<0||(fdaugh>ldaugh))return;
220 for(Int_t j=fdaugh;j<=ldaugh;j++){
221 AliAODMCParticle *pd=(AliAODMCParticle*)clarr->At(j);
222 Int_t pdg=pd->GetPdgCode();
223 fnPhysPrim++;// increment it now and decrement it after all if, in case the daughter is not a phys prim
224 fpxtotdaugh+=pd->Px(); // increment it now and decrement it after all if, in case the daughter is not a phys prim
225 fpytotdaugh+=pd->Py();
226 fpztotdaugh+=pd->Pz();
227
228 if(pdg==211)fnpionPlus++;
229 else if(pdg==-211)fnpionMinus++;
230 else if(pdg==111)fnpionZero++;
231 else if(pdg==321)fnkaonPlus++;
232 else if(pdg==-321)fnkaonMinus++;
233 else if(pdg==310)fnkaonZeroShort++;
234 else if(pdg==2212)fnProton++;
235 else if(pdg==-2212)fnAntiProton++;
236 else if(pdg==-11)fnElePlus++;
237 else if(pdg==11)fnEleMinus++;
238 else if(pdg==-13)fnMuonPlus++;
239 else if(pdg==13)fnMuonMinus++;
240 else if(pdg==12||pdg==14||pdg==-12||pdg==-14)fnNeutrino++;
241 else if(pdg==443)fnJpsi++;
242 else if(pdg!=130&&pdg!=22){
243 fnPhysPrim--;
244 fpxtotdaugh-=pd->Px();
245 fpytotdaugh-=pd->Py();
246 fpztotdaugh-=pd->Pz();
247 GetFinalStateParticles(pd,clarr);
248 }
249 }
250}
251
252//_______________________________________________
253void AliAnalysisTaskSEmcCorr::CheckDzeroChannels(AliAODMCParticle *part,TClonesArray *clarr){
254 Int_t pdg=part->GetPdgCode();
255 if(!(TMath::Abs(pdg)==421))return;
256 fhDzeroDecay->Fill(0);
257 fnpionPlus=0;
258 fnpionMinus=0;
259 fnpionZero=0;
260 fnkaonPlus=0;
261 fnkaonMinus=0;
262 fnkaonZeroShort=0;
263 fnProton=0;
264 fnAntiProton=0;
265 fnElePlus=0;
266 fnEleMinus=0;
267 fnMuonPlus=0;
268 fnMuonMinus=0;
269 fnNeutrino=0;
270 fnJpsi=0;
271 fnPhysPrim=0;
272 fpxtotdaugh=0.;
273 fpytotdaugh=0.;
274 fpztotdaugh=0.;
275
276 GetFinalStateParticles(part,clarr);
277
278 Double_t px=part->Px();
279 Double_t py=part->Py();
280 Double_t pz=part->Pz();
281 // check exclusive channels in the list
282 if(TMath::Abs(px-fpxtotdaugh)<1.e-6&&TMath::Abs(py-fpytotdaugh)<1.e-6&&TMath::Abs(pz-fpztotdaugh)<1.e-6){
283 if(fnPhysPrim==2&&((pdg==421&&fnpionPlus==1&&fnkaonMinus==1)||(pdg==-421&&fnpionMinus==1&&fnkaonPlus==1))){
284 fhDzeroDecay->Fill(6);
285 }
286 if(fnPhysPrim==3&&((pdg==421&&fnpionPlus==1&&fnkaonMinus==1&&fnpionZero==1)||(pdg==-421&&fnpionMinus==1&&fnkaonPlus==1&&fnpionZero==1))){
287 fhDzeroDecay->Fill(7);
288 }
289 if(fnPhysPrim==4&&((pdg==421&&fnpionPlus==2&&fnkaonMinus==1&&fnpionMinus==1)||(pdg==-421&&fnpionMinus==2&&fnkaonPlus==1&&fnpionPlus==1))){
290 fhDzeroDecay->Fill(8);
291 }
292 if(fnPhysPrim==3&&((pdg==421&&fnElePlus==1&&fnkaonMinus==1&&fnNeutrino==1)||(pdg==-421&&fnEleMinus==1&&fnkaonPlus==1&&fnNeutrino==1))){
293 fhDzeroDecay->Fill(9);
294 }
295 if(fnPhysPrim==3&&((pdg==421&&fnMuonPlus==1&&fnkaonMinus==1&&fnNeutrino==1)||(pdg==-421&&fnMuonMinus==1&&fnkaonPlus==1&&fnNeutrino==1))){
296 fhDzeroDecay->Fill(10);
297 }
298 }
299
300 //now check decay to electron + X
301 if(fnElePlus>=1||fnEleMinus>=1){
302 fhDzeroDecay->Fill(11);
303 }
304 // no check number of prong
305 Int_t nprong=fnPhysPrim;
306 nprong-=fnpionZero;
307 nprong-=fnNeutrino;
308 nprong-=fnkaonZeroShort;
309 if(nprong==0){
310 fhDzeroDecay->Fill(1);
311 }
312 else if(nprong==1){// should never happen for D0
313 fhDzeroDecay->Fill(2);
314 }
315 else if(nprong==2){// should never happen for D0
316 fhDzeroDecay->Fill(3);
317 }
318 else if(nprong==3){
319 fhDzeroDecay->Fill(4);
320 }
321 else if(nprong==4){
322 fhDzeroDecay->Fill(5);
323 }
324
325 return;
326}
327
328
329//_______________________________________________
330void AliAnalysisTaskSEmcCorr::CheckDplusChannels(AliAODMCParticle *part,TClonesArray *clarr){
331 Int_t pdg=part->GetPdgCode();
332 if(!(TMath::Abs(pdg)==411))return;
333 fhDplusDecay->Fill(0);
334 fnpionPlus=0;
335 fnpionMinus=0;
336 fnpionZero=0;
337 fnkaonPlus=0;
338 fnkaonMinus=0;
339 fnkaonZeroShort=0;
340 fnProton=0;
341 fnAntiProton=0;
342 fnElePlus=0;
343 fnEleMinus=0;
344 fnMuonPlus=0;
345 fnMuonMinus=0;
346 fnNeutrino=0;
347 fnJpsi=0;
348 fnPhysPrim=0;
349 fpxtotdaugh=0.;
350 fpytotdaugh=0.;
351 fpztotdaugh=0.;
352 GetFinalStateParticles(part,clarr);
353
354 Double_t px=part->Px();
355 Double_t py=part->Py();
356 Double_t pz=part->Pz();
357 // check exclusive channels in the list
358 if(TMath::Abs(px-fpxtotdaugh)<1.e-6&&TMath::Abs(py-fpytotdaugh)<1.e-6&&TMath::Abs(pz-fpztotdaugh)<1.e-6){
359 if(fnPhysPrim==3&&((pdg==411&&fnpionPlus==2&&fnkaonMinus==1)||(pdg==-411&&fnpionMinus==2&&fnkaonPlus==1))){
360 fhDplusDecay->Fill(6);
361 }
362 if(fnPhysPrim==3&&((pdg==411&&fnkaonZeroShort==1&&fnpionPlus==1&&fnpionZero==1)||(pdg==-411&&fnkaonZeroShort==1&&fnpionMinus==1&&fnpionZero==1))){
363 fhDplusDecay->Fill(7);
364 }
365 if(fnPhysPrim==4&&((pdg==411&&fnpionPlus==2&&fnkaonMinus==1&&fnpionZero==1)||(pdg==-411&&fnpionMinus==2&&fnkaonPlus==1&&fnpionZero==1))){
366 fhDzeroDecay->Fill(8);
367 }
368 if(fnPhysPrim==3&&((pdg==411&&fnElePlus==1&&fnkaonZeroShort==1&&fnNeutrino==1)||(pdg==-411&&fnEleMinus==1&&fnkaonZeroShort==1&&fnNeutrino==1))){
369 fhDplusDecay->Fill(9);
370 }
371 if(fnPhysPrim==3&&((pdg==411&&fnMuonPlus==1&&fnkaonZeroShort==1&&fnNeutrino==1)||(pdg==-411&&fnMuonMinus==1&&fnkaonZeroShort==1&&fnNeutrino==1))){
372 fhDplusDecay->Fill(10);
373 }
374 }
375
376 //now check decay to electron + X
377 if(fnElePlus>=1||fnEleMinus>=1){
378 fhDplusDecay->Fill(11);
379 }
380 // no check number of prong
381 Int_t nprong=fnPhysPrim;
382 nprong-=fnpionZero;
383 nprong-=fnNeutrino;
384 nprong-=fnkaonZeroShort;
385 if(nprong==0){
386 fhDplusDecay->Fill(1);
387 }
388 else if(nprong==1){// should never happen for D0
389 fhDplusDecay->Fill(2);
390 }
391 else if(nprong==2){// should never happen for D0
392 fhDplusDecay->Fill(3);
393 }
394 else if(nprong==3){
395 fhDplusDecay->Fill(4);
396 }
397 else if(nprong==4){
398 fhDplusDecay->Fill(5);
399 }
400
401 return;
402}
403//________________________________________________________________________
404// void AliAnalysisTaskSEmcCorr::CheckLambdaChannels(AliAODMCParticle *part,TClonesArray *clarr){
405// // not implemented yet
406// }
407
408//________________________________________________________________________
409// void AliAnalysisTaskSEmcCorr::CheckBallChannels(AliAODMCParticle *part,TClonesArray *clarr){
410// // not implemented yet
411// }
d33db6d9
AR		 412
413
414//________________________________________________________________________
415void AliAnalysisTaskSEmcCorr::UserCreateOutputObjects(){
416
417 if(!fArrayDaugh) fArrayDaugh=new TArrayI(20);
418
419 if(!fArrayAssoc)fArrayAssoc=new TArrayI(200);
3333be85 420
d33db6d9
AR		 421
422 fNentries=new TH1F("nentriesChFr", "Analyzed sample properties", 21,0.5,21.5);
423 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
424 fNentries->GetXaxis()->SetBinLabel(2,"nEvSel");
425 fNentries->GetXaxis()->SetBinLabel(3,"nEvPile-up Rej");
426 fNentries->GetXaxis()->SetBinLabel(4,"nEvGoodVtxS");
427 fNentries->GetXaxis()->SetBinLabel(5,"nEvRejVtxZ");
428 fNentries->GetXaxis()->SetBinLabel(6,"nD0");
3333be85 429 fNentries->GetXaxis()->SetBinLabel(7,"nEventsSelGenName");
d33db6d9
AR		 430
431 fhNprongsD0=new TH1F("fhNprongsD0","fhNprongsD0",20,-0.5,19.5);
432 fhNprongsD0chargedOnly=new TH1F("fhNprongsD0chargedOnly","fhNprongsD0chargedOnly",20,-0.5,19.5);
433 fhNprongsD0chargedRef=new TH1F("fhNprongsD0chargedRef","fhNprongsD0chargedRef",20,-0.5,19.5);
434
435
2942f542 436 Int_t nbinsCorrMC[8]={15,20,20,20,50,63,20,11};
d33db6d9
AR		 437 Double_t binlowCorrMC[8]={-7.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
438 Double_t binupCorrMC[8]={7.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,9.5};
439 fhMCcorrel=new THnSparseF("fhMCcorrel","fhMCcorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMC,binlowCorrMC,binupCorrMC);
440
2942f542 441 Int_t nbinsTrigMC[3]={15,20,20};
d33db6d9
AR		 442 Double_t binlowTrigMC[3]={-7.5,0.,-1.};
443 Double_t binupTrigMC[3]={7.5,20.,1.};
444 fhMCtrigPart=new THnSparseF("fhMCtrigPart","fhMCcorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMC,binlowTrigMC,binupTrigMC);
445
3333be85 446 fhDzeroDecay=new TH1F("fhDzeroDecay","Dzero Decay channels",14,-0.5,13.5);
447 fhDzeroDecay->GetXaxis()->SetBinLabel(1,"All D0");
448 fhDzeroDecay->GetXaxis()->SetBinLabel(2,"0 prong");
449 fhDzeroDecay->GetXaxis()->SetBinLabel(3,"1 prong");
450 fhDzeroDecay->GetXaxis()->SetBinLabel(4,"2 prong");
451 fhDzeroDecay->GetXaxis()->SetBinLabel(5,"3 prong");
452 fhDzeroDecay->GetXaxis()->SetBinLabel(6,"4 prong");
453 fhDzeroDecay->GetXaxis()->SetBinLabel(7,"K pi (3.88)");
454 fhDzeroDecay->GetXaxis()->SetBinLabel(8,"K pi pi0 (13.9)");
455 fhDzeroDecay->GetXaxis()->SetBinLabel(9,"K pi pi pi (8.09)");
456 fhDzeroDecay->GetXaxis()->SetBinLabel(10,"e K nu (3.55)");
457 fhDzeroDecay->GetXaxis()->SetBinLabel(11,"mu K nu (3.31)");
458 fhDzeroDecay->GetXaxis()->SetBinLabel(12,"e X (6.49)");
459
460
461 fhDplusDecay=new TH1F("fhDplusDecay","Dplus Decay channels",14,-0.5,13.5);
462 fhDplusDecay->GetXaxis()->SetBinLabel(1,"All D+");
463 fhDplusDecay->GetXaxis()->SetBinLabel(2,"0 prong");
464 fhDplusDecay->GetXaxis()->SetBinLabel(3,"1 prong");
465 fhDplusDecay->GetXaxis()->SetBinLabel(4,"2 prong");
466 fhDplusDecay->GetXaxis()->SetBinLabel(5,"3 prong");
467 fhDplusDecay->GetXaxis()->SetBinLabel(6,"4 prong");
468 fhDplusDecay->GetXaxis()->SetBinLabel(7,"K pi pi (9.4)");
469 fhDplusDecay->GetXaxis()->SetBinLabel(8,"K0s pi pi0 (6.9)");
470 fhDplusDecay->GetXaxis()->SetBinLabel(9,"K pi pi pi0 (6.08)");
471 fhDplusDecay->GetXaxis()->SetBinLabel(10,"e K0bar nu (8.83)");
472 fhDplusDecay->GetXaxis()->SetBinLabel(11,"mu K0bar nu (9.4)");
473 fhDplusDecay->GetXaxis()->SetBinLabel(12,"e X (16.07");
474
475 fhLambdaCDecay=new TH1F("fhLambdaCDecay","LambdaC Decay channels",10,-0.5,9.5);
476 fhLambdaCDecay->GetXaxis()->SetBinLabel(1,"All LambdaC");
477 fhLambdaCDecay->GetXaxis()->SetBinLabel(2,"p K0s");
478 fhLambdaCDecay->GetXaxis()->SetBinLabel(3,"p K pi");
479
480 fhAllBDecay=new TH1F("fhAllBDecay","All bhadron (also LambdaB) Decay channels",10,-0.5,9.5);
481 fhAllBDecay->GetXaxis()->SetBinLabel(1,"All b-hadrons");
482 fhAllBDecay->GetXaxis()->SetBinLabel(2,"D0bar X (59.6)");
483 fhAllBDecay->GetXaxis()->SetBinLabel(3,"D0bar D0 X (5.1) ");
484 fhAllBDecay->GetXaxis()->SetBinLabel(4,"D+- D0 X (2.7)");
485 fhAllBDecay->GetXaxis()->SetBinLabel(5,"D- X (22.7)");
486 fhAllBDecay->GetXaxis()->SetBinLabel(6,"J/Psi X (1.16)");
d33db6d9 487
3333be85 488
489
2942f542 490 Int_t nbinsCorrMChadron[8]={6,20,20,20,50,63,20,7};
d33db6d9
AR		 491 Double_t binlowCorrMChadron[8]={-1.5,0.,-1.,0.,-5.,-0.5*TMath::Pi(),-2,-1.5};
492 Double_t binupCorrMChadron[8]={4.5,20.,1.,2.,5.,1.5*TMath::Pi(),2.,5.5};
493 fhMChadroncorrel=new THnSparseF("fhMChadroncorrel","fhMChadroncorrel;pdg;ptTrig;etaTrig;ptAss;etaAss;deltaPhi;deltaEta;pdgAss;",8,nbinsCorrMChadron,binlowCorrMChadron,binupCorrMChadron);
494
2942f542 495 Int_t nbinsTrigMChadron[3]={6,20,20};
d33db6d9
AR		 496 Double_t binlowTrigMChadron[3]={-1.5,0.,-1.};
497 Double_t binupTrigMChadron[3]={4.5,20.,1.};
498 fhMChadrontrigPart=new THnSparseF("fhMChadrontrigPart","fhMChadroncorrel;pdg;ptTrig;etaTrig;",3,nbinsTrigMChadron,binlowTrigMChadron,binupTrigMChadron);
499
500
501 PostData(1,fNentries);
502 PostData(2,fhNprongsD0);
503 PostData(3,fhNprongsD0chargedOnly);
504 PostData(4,fhNprongsD0chargedRef);
505 PostData(5,fhMCcorrel);
506 PostData(6,fhMCtrigPart);
507 PostData(7,fhMChadroncorrel);
508 PostData(8,fhMChadrontrigPart);
3333be85 509 PostData(9,fhDzeroDecay);
510 PostData(10,fhDplusDecay);
511 PostData(11,fhLambdaCDecay);
512 PostData(12,fhAllBDecay);
513
514
d33db6d9
AR		 515}
516
517
518
519//________________________________________________________________________
520void AliAnalysisTaskSEmcCorr::UserExec(Option_t */*option*/){
521
522
523
524 // Execute analysis for current event:
525 // heavy flavor candidates association to MC truth
526
527 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
528 if (!aod) {
529 Printf("ERROR: aod not available");
530 return;
531 }
532
533
d33db6d9
AR		 534
535 fNentries->Fill(1);
536
537 // AOD primary vertex
538
539
540 if(!fOnlyKine){
541 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
542 TString primTitle = vtx1->GetTitle();
543 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
544
545 fNentries->Fill(3);
546
547 }
548 else {
549 PostData(1,fNentries);
550 return;
551
552 }
553 }
3333be85 554
d33db6d9
AR		 555
556 TClonesArray *arrayMC=0x0;
557 AliAODMCHeader *aodmcHeader=0x0;
558 Double_t vtxTrue[3];
559
560
561 if(fReadMC){
562 // aod->GetList()->Print();
563 // load MC particles
3333be85 564 if(fFastSimul){
565 arrayMC =
566 (TClonesArray*)aod->GetList()->At(23);//FindObject("mcparticles");
567 }
568 else{
569 arrayMC =
570 (TClonesArray*)aod->GetList()->FindObject("mcparticles");
571 }
d33db6d9
AR		 572 if(!arrayMC) {
573 Printf("AliAnalysisTaskSEmcCorr::UserExec: MC particles branch not found!\n");
574 fNentries->Fill(2);
575 PostData(1,fNentries);
576 return;
577 }
578 // load MC header
579 aodmcHeader =
580 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
581 if(!aodmcHeader) {
582 Printf("AliAnalysisTaskSEmcCorr::UserExec: MC header branch not found!\n");
3333be85 583 fNentries->Fill(3);
d33db6d9
AR		 584 PostData(1,fNentries);
585 return;
586 }
3333be85 587
588 if(!(fGeneratorString.IsNull())){
589 TString strGenTitle=((AliGenEventHeader*)aodmcHeader->GetCocktailHeader(0))->GetTitle();
590 // Printf("Gen Title: %s",strGenTitle.Data());
591 if(!fGeneratorString.Contains(strGenTitle.Data())){
592 // Printf("Event rejected from generator title");
593 return;
594 }
595 }
d33db6d9
AR		 596 fNentries->Fill(4);
597 // printf("N MC entries: %d \n",arrayMC->GetEntriesFast());
598 // MC primary vertex
599 aodmcHeader->GetVertex(vtxTrue);
3333be85 600
601
d33db6d9 602 // FILL HISTOS FOR D0 mesons, c quarks and D0 from B properties
3333be85 603 if(fDoHFCorrelations) SelectAssociatedParticles(arrayMC);
d33db6d9 604 for(Int_t jp=0;jp<arrayMC->GetEntriesFast();jp++){
d33db6d9
AR		 605 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(jp);
606 Int_t pdg=TMath::Abs(part->GetPdgCode());
3333be85 607 if(fDoHFCorrelations){
608 if(pdg==421||pdg==411||pdg==413){
609
610 if(TMath::Abs(part->Y())>fYrangeTrig)continue;
611 if(TMath::Abs(part->Pt())<fminDpt)continue;
612 FillSkipParticleArray(part,arrayMC,jp);
613 FillCorrelationPlots(part,arrayMC);
614 }
615 else if(pdg==11){
616 if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
617 if(TMath::Abs(part->Pt())<fminDpt)continue;
618 FillSkipParticleArray(part,arrayMC,jp);
619 FillCorrelationPlots(part,arrayMC);
620 }
d33db6d9 621 }
3333be85 622 if(fCheckDecay){
623 CheckDzeroChannels(part,arrayMC);
624 CheckDplusChannels(part,arrayMC);
625 // CheckLambdaChannels(part,arrayMC);
626 // CheckBallChannels(part,arrayMC);
627
628 // NOW STUDY OF N-PRONGS (OLD PART OF THE TASK)
629 if(pdg==421){
630 Bool_t semilept=kFALSE;
631 Int_t nprongsD0=0;
632 Int_t nprongsD0charged=0;
633
634 Int_t lab1=part->GetDaughter(0);
635 Int_t lab2=part->GetDaughter(1);
636 if(lab1>=0&&lab2>=0){
637 for(Int_t jd=lab1;jd<=lab2;jd++){
638 AliAODMCParticle *d=(AliAODMCParticle*)arrayMC->At(jd);
639 Int_t pdgd=TMath::Abs(d->GetPdgCode());
640 if(pdgd>=11&&pdgd<=16){
641 semilept=kTRUE;
642 break;
643 }
644 if(pdgd==211||pdgd==111||pdgd==321||pdgd==311||pdgd==310||pdgd==130){
645 nprongsD0++;
646 }
647 if(pdgd==211||pdgd==321){
648 nprongsD0charged++;
649 }
650 else{
651 Int_t lab1d=d->GetDaughter(0);
652 Int_t lab2d=d->GetDaughter(1);
653 if(lab1d>=0&&lab2d>=0){
654 for(Int_t jdd=lab1d;jdd<=lab2d;jdd++){
655 AliAODMCParticle *dd=(AliAODMCParticle*)arrayMC->At(jdd);
656 Int_t pdgdd=TMath::Abs(dd->GetPdgCode());
657 if(pdgdd==211||pdgdd==111||pdgdd==321||pdgdd==311||pdgdd==310||pdgdd==130){
658 nprongsD0++;
659 }
660 if(pdgd==211||pdgd==321){
661 nprongsD0charged++;
662 }
d33db6d9 663 }
d33db6d9
AR		 664 }
665 }
666 }
3333be85 667 }
668 if(!semilept){
669 fhNprongsD0->Fill(nprongsD0);
670 fhNprongsD0chargedOnly->Fill(nprongsD0charged);
671 fhNprongsD0chargedRef->Fill(nprongsD0charged);
d33db6d9 672 }
d33db6d9
AR		 673 }
674 }
d33db6d9
AR		 675 }
676
677 if(fDoHadronHadron){
678 for(Int_t jp=0;jp<fLastAss;jp++){
679
680 AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->At(fArrayAssoc->At(jp));
681 Int_t pdg=TMath::Abs(part->GetPdgCode());
682 if(pdg==11||pdg==13||pdg==211||pdg==2212||pdg==321){
683 if(TMath::Abs(part->Eta())>fEtarangeEleTrig)continue;
684 if(TMath::Abs(part->Pt())<fminDpt)continue;
685 FillSkipParticleArray(part,arrayMC,fArrayAssoc->At(jp));
686 FillCorrelationPlotsHadrons(part,arrayMC);
687 }
688 }
689
690 }
691
692 }
693
694
695 PostData(1,fNentries);
696 PostData(2,fhNprongsD0);
697 PostData(3,fhNprongsD0chargedOnly);
698 PostData(4,fhNprongsD0chargedRef);
699 PostData(5,fhMCcorrel);
700 PostData(6,fhMCtrigPart);
701 PostData(7,fhMChadroncorrel);
702 PostData(8,fhMChadrontrigPart);
3333be85 703 PostData(9,fhDzeroDecay);
704 PostData(10,fhDplusDecay);
705 PostData(11,fhLambdaCDecay);
706 PostData(12,fhAllBDecay);
d33db6d9
AR		 707}
708
709
710
711void AliAnalysisTaskSEmcCorr::FillSkipParticleArray(AliAODMCParticle *part,TClonesArray *arrayMC,Int_t partPos){
712 fArrayDaugh->Reset(0);
713 fArrayDaugh->AddAt(partPos,0);
714 flastdaugh=1;
715
716 // Loop on daughters and nephews
717 if(part->GetNDaughters()>0){
718 Int_t pos1=part->GetDaughter(0);
719 Int_t pos2=part->GetDaughter(1);
720 if(pos2<0)pos2=pos1;
721 if(pos1>=0){
722 for(Int_t jd=pos1;jd<=pos2;jd++){
723 AliAODMCParticle *pd=(AliAODMCParticle*)arrayMC->At(jd);
724 Int_t pdgd=TMath::Abs(pd->GetPdgCode());
725 if(pdgd==11||pdgd==13||pdgd==211||pdgd==321||pdgd==2212){
726 fArrayDaugh->AddAt(jd,flastdaugh);
727 flastdaugh++;
728 }
729 else if(pdgd!=12&&pdgd!=14&&pdgd!=16&&pdgd!=310&&pdgd!=111){// CHECK NEPHEWS (FOR RESONANT CHANNELS)
730 Int_t nephw=pd->GetNDaughters();
731 if(nephw>0){
732 Int_t posN1=pd->GetDaughter(0);
733 Int_t posN2=pd->GetDaughter(1);
734 if(posN2<0)posN2=posN1;
735 if(posN1>=0){
736 for(Int_t jN=posN1;jN<=posN2;jN++){
737 AliAODMCParticle *pN=(AliAODMCParticle*)arrayMC->At(jN);
738 Int_t pdgN=TMath::Abs(pN->GetPdgCode());
739 if(pdgN==11||pdgN==13||pdgN==211||pdgN==321||pdgN==2212){
740 fArrayDaugh->AddAt(jN,flastdaugh);
741 flastdaugh++;
742 }
743 else if(pdgN!=12&&pdgN!=14&&pdgN!=16&&pdgN!=310&&pdgN!=111){// CHECK one step more (FOR RESONANT CHANNELS)
744
745 Int_t gnephw=pN->GetNDaughters();
746 if(gnephw>0){
747 Int_t posGN1=pN->GetDaughter(0);
748 Int_t posGN2=pN->GetDaughter(1);
749 if(posGN2<0)posGN2=posGN1;
750 if(posGN1>=0){
751 for(Int_t jGN=posGN1;jGN<=posGN2;jGN++){
752 AliAODMCParticle *pGN=(AliAODMCParticle*)arrayMC->At(jGN);
753 Int_t pdgGN=TMath::Abs(pGN->GetPdgCode());
754 if(pdgGN==11||pdgGN==13||pdgGN==211||pdgGN==321||pdgGN==2212){
755 fArrayDaugh->AddAt(jGN,flastdaugh);
756 flastdaugh++;
757 }
758 }
759 }
760 }
761 }
762 }
763 }
764 }
765 }
766 }
767 }
768 }
769 // printf("particles to be skipped %d\n",flastdaugh++);
770}
771
772void AliAnalysisTaskSEmcCorr::SelectAssociatedParticles(TClonesArray *arrayMC){
773 fLastAss=0;
774 fArrayAssoc->Reset(0);
775 Double_t pt,eta;
776 for(Int_t jAss=0;jAss<arrayMC->GetEntriesFast();jAss++){
777 AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(jAss);
778 Int_t pdgAss=TMath::Abs(partAss->GetPdgCode());
779 if(pdgAss==11||pdgAss==13||pdgAss==211||pdgAss==321||pdgAss==2212){// check type: keep only e,mu,pi,K+,protons
780 pt=partAss->Pt();
781 eta=TMath::Abs(partAss->Eta());
782 if(pt<fminPtass)continue;
783 if(eta>fmaxEtaAss)continue;
784 if(!partAss->IsPhysicalPrimary()){
785 if(pdgAss!=11) continue;// keep ele
786 }
787 fArrayAssoc->AddAt(jAss,fLastAss);
788 fLastAss++;
3333be85 789 if(fLastAss==fArrayAssoc->GetSize()-1){
790 fArrayAssoc->Set(fArrayAssoc->GetSize()+200);
791 }
d33db6d9
AR		 792 }
793 }
794 // printf("n ass part: %d\n",fLastAss++);
795}
796
797void AliAnalysisTaskSEmcCorr::FillCorrelationPlots(AliAODMCParticle *part,TClonesArray *arrayMC){
798
799
800 Bool_t hasToSkip=kFALSE;
801 Int_t index=0,softpi=-1;
802 Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
803 if(pdgTrig==421){
804 pdgTrig=1;
805 }
806 else if(pdgTrig==411){
807 pdgTrig=2;
808 }
809 else if(pdgTrig==413){
810 pdgTrig=3;
811 }
812 else if(pdgTrig==11){
813
814 pdgTrig=4;
815 }
816 else pdgTrig=0;
817 // check if it is prompt or from B and, for electrons the origin
818 Int_t labmum=part->GetMother();
819 if(labmum>=0){
820 AliAODMCParticle *moth=(AliAODMCParticle*)arrayMC->At(labmum);
821 Int_t pdgmum=TMath::Abs(moth->GetPdgCode());
822 if(pdgTrig==4){// for electrons go to the parent hadron
823 // labmum=part->GetMother();
824 // moth=(AliAODMCParticle*)arrayMC->At(labmum);
825 // pdgmum=TMath::Abs(moth->GetPdgCode());
826 if(!(pdgmum==5||pdgmum==4||(400<pdgmum&&pdgmum<600)||(4000<pdgmum&&pdgmum<6000))){// NON HF electron
827 if(pdgmum==22){// GAMMA CONV
828 pdgTrig=5;
829 }
830 else if(pdgmum==111||pdgmum==113||pdgmum==221||pdgmum==223){ //DALITZ DECAY
831 pdgTrig=6;
832 }
833 else pdgTrig=7;
834 }
835 }
836 if(pdgmum==413){ // D from Dstar: check soft pion and go to the grandmother
837 if(pdgTrig==1){
838 for(Int_t isp=moth->GetDaughter(0);isp<=moth->GetDaughter(1);isp++){
839 AliAODMCParticle *sfp=(AliAODMCParticle*)arrayMC->At(isp);
840 Int_t pdgsp=TMath::Abs(sfp->GetPdgCode());
841 if(pdgsp==211)softpi=isp;
842 }
843 }
844 labmum=moth->GetMother();
3333be85 845 if(labmum>=0){
846 moth=(AliAODMCParticle*)arrayMC->At(labmum);
847 pdgmum=TMath::Abs(moth->GetPdgCode());
848 }
849 else pdgmum=-1;
d33db6d9
AR		 850 }
851 else if(pdgmum==423){// D*0 -> go to the mother
852 labmum=moth->GetMother();
3333be85 853 if(labmum>=0){
854 moth=(AliAODMCParticle*)arrayMC->At(labmum);
855 pdgmum=TMath::Abs(moth->GetPdgCode());
856 }
857 else pdgmum=-1;
d33db6d9
AR		 858 }
859 if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
860 pdgTrig*=-1;
861 }
862 else if(pdgmum==4){
863 labmum=moth->GetMother();
864 if(labmum>=0){
865 moth=(AliAODMCParticle*)arrayMC->At(labmum);
866 pdgmum=TMath::Abs(moth->GetPdgCode());
867 if(pdgmum==5||(500<pdgmum&&pdgmum<600)||(5000<pdgmum&&pdgmum<6000)){
868 // Printf("c quark from b quark");
869 pdgTrig*=-1;
870 }
871 }
872 }
873 else {
874 // printf("Charm not from charm or beauty \n");
875 }
876 }
877 Double_t phiTrig=part->Phi();
878 Double_t etaTrig=part->Eta();
879 Double_t deltaPhi;
2942f542 880 Double_t point[8]={ static_cast<Double_t>(pdgTrig),part->Pt(),etaTrig,0,0,0,0,0};//pdg,ptTrig,etaTrig,ptAss,etaAss,deltaPhi,deltaEta,pdgAss
d33db6d9
AR		 881 fhMCtrigPart->Fill(point);
882 for(Int_t jAss=0;jAss<fLastAss;jAss++){
883 index=fArrayAssoc->At(jAss);
884 hasToSkip=kFALSE;
885 for(Int_t jd=0;jd<flastdaugh;jd++){
886 if(index==fArrayDaugh->At(jd)){
887 hasToSkip=kTRUE;
888 break;
889 }
890 }
891 if(hasToSkip)continue;
892 AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
893 point[3]=partAss->Pt();
894 point[4]=partAss->Eta();
895 deltaPhi=partAss->Phi()-phiTrig;
896 if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
897 else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
898 point[5]=deltaPhi;
899 point[6]=partAss->Eta()-etaTrig;
900 UInt_t pdgassmum=TMath::Abs(partAss->GetPdgCode());
901 if(index==softpi){
902 point[7]=-1;
903 }
904 else if(pdgassmum==11){
905 if(partAss->IsPhysicalPrimary())point[7]=1;
906 Int_t elemum=partAss->GetMother();
907 if(elemum>0){
908 AliAODMCParticle *partEleMum=(AliAODMCParticle*)arrayMC->At(elemum);
909 UInt_t pdgelemum=TMath::Abs(partEleMum->GetPdgCode());
910 if(pdgelemum==111||pdgelemum==113||pdgelemum==221||pdgelemum==223){// DALITZ DECAY
911 point[7]=6;
912 }
913 else if(pdgelemum==22){// GAMMA CONV
914 point[7]=7;
915 }
916 else if(pdgelemum==4||pdgelemum==5||(pdgelemum>400&&pdgelemum<600)||(pdgelemum>4000&&pdgelemum<6000)){// HF e
917 point[7]=8;
918 }
919 else point[7]=9;
920 }
921 }
922 else if(pdgassmum==13){
923 point[7]=2;
924 }
925 else if(pdgassmum==211){
926 point[7]=3;
927 }
928 else if(pdgassmum==321){
929 point[7]=4;
930 }
931 else if(pdgassmum==2212){
932 point[7]=5;
933 }
934 else point[7]=0;
935
936 fhMCcorrel->Fill(point);
937 }
938
939}
940
941
942
943//______________________________________________________________
944void AliAnalysisTaskSEmcCorr::FillCorrelationPlotsHadrons(AliAODMCParticle *part,TClonesArray *arrayMC){
945
946
947 Bool_t hasToSkip=kFALSE;
948 Int_t index=0;
949 Int_t pdgTrig=TMath::Abs(part->GetPdgCode());
950 if(pdgTrig==11){
951 pdgTrig=0;
952 }
953 else if(pdgTrig==13){
954 pdgTrig=1;
955 }
956 else if(pdgTrig==211){
957 pdgTrig=2;
958 }
959 else if(pdgTrig==321){
960 pdgTrig=3;
961 }
962 else if(pdgTrig==2212){
963 pdgTrig=4;
964 }
965 else pdgTrig=-1;
966
967
968 Double_t phiTrig=part->Phi();
969 Double_t etaTrig=part->Eta();
970 Double_t deltaPhi;
2942f542 971 Double_t point[8]={ static_cast<Double_t>(pdgTrig),part->Pt(),etaTrig,0,0,0,0,0};//pdg,ptTrig,etaTrig,ptAss,etaAss,deltaPhi,deltaEta,pdgAss
d33db6d9
AR		 972 fhMChadrontrigPart->Fill(point);
973 for(Int_t jAss=0;jAss<fLastAss;jAss++){
974 index=fArrayAssoc->At(jAss);
975 hasToSkip=kFALSE;
976 for(Int_t jd=0;jd<flastdaugh;jd++){
977 if(index==fArrayDaugh->At(jd)){
978 hasToSkip=kTRUE;
979 break;
980 }
981 }
982 if(hasToSkip)continue;
983 AliAODMCParticle *partAss=(AliAODMCParticle*)arrayMC->At(index);
984 point[3]=partAss->Pt();
985 if(point[3]>point[1])continue;
986 point[4]=partAss->Eta();
987 deltaPhi=partAss->Phi()-phiTrig;
988 if(deltaPhi<-0.5*TMath::Pi())deltaPhi+=2.*TMath::Pi();
989 else if(deltaPhi>1.5*TMath::Pi())deltaPhi-=2.*TMath::Pi();
990 point[5]=deltaPhi;
991 point[6]=partAss->Eta()-etaTrig;
992 point[7]=TMath::Abs(partAss->GetPdgCode());
993 if(point[7]==11){
994 point[7]=0;
995 }
996 else if(point[7]==13){
997 point[7]=1;
998 }
999 else if(point[7]==211){
1000 point[7]=2;
1001 }
1002 else if(point[7]==321){
1003 point[7]=3;
1004 }
1005 else if(point[7]==2212){
1006 point[7]=4;
1007 }
1008 else point[7]=-1;
1009
1010 fhMChadroncorrel->Fill(point);
1011 }
1012
1013}
1014
1015//_______________________________________________________________
1016void AliAnalysisTaskSEmcCorr::Terminate(const Option_t*){
1017 //TERMINATE METHOD: NOTHING TO DO
1018
1019
1020
1021}
1022