]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
consolidate zero-length arrays (aka struct hack)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskCombinHF.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2018, 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 /* $Id: $ */
17
18 //*************************************************************************
19 // Class AliAnalysisTaskCombinHF
20 // AliAnalysisTaskSE to build D meson candidates by combining tracks
21 //  background is computed LS and track rotations is
22 // Authors: F. Prino, A. Rossi
23 /////////////////////////////////////////////////////////////
24
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TH3F.h>
29 #include <THnSparse.h>
30
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 #include "AliPIDResponse.h"
34 #include "AliAODHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODVertex.h"
39 #include "AliAODTrack.h"
40 #include "AliVertexingHFUtils.h"
41 #include "AliAnalysisTaskCombinHF.h"
42
43 ClassImp(AliAnalysisTaskCombinHF)
44
45
46 //________________________________________________________________________
47 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
48   AliAnalysisTaskSE(),
49   fOutput(0x0), 
50   fHistNEvents(0x0),
51   fHistTrackStatus(0x0),
52   fHistCheckOrigin(0x0),
53   fHistCheckOriginSel(0x0),
54   fHistCheckDecChan(0x0),
55   fHistCheckDecChanAcc(0x0),
56   fPtVsYGen(0x0),
57   fPtVsYGenLimAcc(0x0),
58   fPtVsYGenAcc(0x0),
59   fPtVsYReco(0x0),
60   fMassVsPtVsY(0x0),
61   fMassVsPtVsYRot(0x0),
62   fMassVsPtVsYLSpp(0x0),
63   fMassVsPtVsYLSmm(0x0),
64   fMassVsPtVsYSig(0x0),
65   fMassVsPtVsYRefl(0x0),
66   fMassVsPtVsYBkg(0x0),
67   fNSelected(0x0),
68   fNormRotated(0x0),
69   fDeltaMass(0x0),
70   fDeltaMassFullAnalysis(0x0),
71   fFilterMask(BIT(4)),
72   fTrackCutsAll(0x0),
73   fTrackCutsPion(0x0),
74   fTrackCutsKaon(0x0),
75   fPidHF(new AliAODPidHF()),
76   fAnalysisCuts(0x0),
77   fMinMass(1.750),
78   fMaxMass(2.150),
79   fEtaAccCut(0.9),
80   fPtAccCut(0.1),
81   fNRotations(9),
82   fMinAngleForRot(5*TMath::Pi()/6),
83   fMaxAngleForRot(7*TMath::Pi()/6),
84   fMinAngleForRot3(2*TMath::Pi()/6),
85   fMaxAngleForRot3(4*TMath::Pi()/6),
86   fCounter(0x0),
87   fMeson(kDzero),
88   fReadMC(kFALSE),
89   fPromptFeeddown(kPrompt),
90   fGoUpToQuark(kTRUE),
91   fFullAnalysis(0),
92   fmaxPforIDPion(0.6),
93   fmaxPforIDKaon(2.),  
94   fKeepNegID(kFALSE)  
95 {
96   // default constructor
97 }
98
99 //________________________________________________________________________
100 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
101   AliAnalysisTaskSE("DmesonCombin"),
102   fOutput(0x0), 
103   fHistNEvents(0x0),
104   fHistTrackStatus(0x0),
105   fHistCheckOrigin(0x0),
106   fHistCheckOriginSel(0x0),
107   fHistCheckDecChan(0x0),
108   fHistCheckDecChanAcc(0x0),
109   fPtVsYGen(0x0),
110   fPtVsYGenLimAcc(0x0),
111   fPtVsYGenAcc(0x0),
112   fPtVsYReco(0x0),
113   fMassVsPtVsY(0x0),
114   fMassVsPtVsYRot(0x0),
115   fMassVsPtVsYLSpp(0x0),
116   fMassVsPtVsYLSmm(0x0),
117   fMassVsPtVsYSig(0x0),
118   fMassVsPtVsYRefl(0x0),
119   fMassVsPtVsYBkg(0x0),
120   fNSelected(0x0),
121   fNormRotated(0x0),
122   fDeltaMass(0x0),
123   fDeltaMassFullAnalysis(0x0),
124   fFilterMask(BIT(4)),
125   fTrackCutsAll(0x0),
126   fTrackCutsPion(0x0),
127   fTrackCutsKaon(0x0),
128   fPidHF(new AliAODPidHF()),
129   fAnalysisCuts(analysiscuts),
130   fMinMass(1.750),
131   fMaxMass(2.150),
132   fEtaAccCut(0.9),
133   fPtAccCut(0.1),
134   fNRotations(9),
135   fMinAngleForRot(5*TMath::Pi()/6),
136   fMaxAngleForRot(7*TMath::Pi()/6),
137   fMinAngleForRot3(2*TMath::Pi()/6),
138   fMaxAngleForRot3(4*TMath::Pi()/6),
139   fCounter(0x0),
140   fMeson(meson),
141   fReadMC(kFALSE),
142   fPromptFeeddown(1),
143   fGoUpToQuark(kTRUE),
144   fFullAnalysis(0),
145   fmaxPforIDPion(0.6),
146   fmaxPforIDKaon(2.),  
147   fKeepNegID(kFALSE)  
148
149 {
150   // standard constructor
151
152   DefineOutput(1,TList::Class());  //My private output
153   DefineOutput(2,AliNormalizationCounter::Class());
154  }
155
156 //________________________________________________________________________
157 AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
158 {
159   //
160   // Destructor
161   //
162   delete fOutput;
163   delete fHistNEvents;
164   delete fHistTrackStatus;
165   delete fHistCheckOrigin;
166   delete fHistCheckOriginSel;
167   delete fHistCheckDecChan;
168   delete fHistCheckDecChanAcc;
169   delete fPtVsYGen;
170   delete fPtVsYGenLimAcc;
171   delete fPtVsYGenAcc;
172   delete fPtVsYReco;
173   delete fMassVsPtVsY; 
174   delete fMassVsPtVsYLSpp;
175   delete fMassVsPtVsYLSmm;
176   delete fMassVsPtVsYRot;
177   delete fMassVsPtVsYSig;
178   delete fMassVsPtVsYRefl;
179   delete fMassVsPtVsYBkg;
180   delete fNSelected;
181   delete fNormRotated;
182   delete fDeltaMass;
183   delete fCounter;
184   delete fTrackCutsAll;
185   delete fTrackCutsPion;
186   delete fTrackCutsKaon;
187   delete fPidHF;
188   delete fAnalysisCuts;
189 }
190
191 //________________________________________________________________________
192 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
193 {
194   // Create the output container
195   //
196   if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
197
198   fOutput = new TList();
199   fOutput->SetOwner();
200   fOutput->SetName("OutputHistos");
201
202   fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5);
203   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
204   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
205   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
206   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
207   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
208   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
209   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
210   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
211
212   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
213   fHistNEvents->Sumw2();
214   fHistNEvents->SetMinimum(0);
215   fOutput->Add(fHistNEvents);
216
217   fHistTrackStatus  = new TH1F("hTrackStatus", "",8,-0.5,7.5);
218   fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
219   fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
220   fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
221   fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
222   fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
223   fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
224   fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
225   fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
226   fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
227   fHistTrackStatus->Sumw2();
228   fHistTrackStatus->SetMinimum(0);
229   fOutput->Add(fHistTrackStatus);
230
231   if(fReadMC){
232
233     fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
234     fHistCheckOrigin->Sumw2();
235     fHistCheckOrigin->SetMinimum(0);
236     fOutput->Add(fHistCheckOrigin);
237
238     fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
239     fHistCheckOriginSel->Sumw2();
240     fHistCheckOriginSel->SetMinimum(0);
241     fOutput->Add(fHistCheckOriginSel);
242
243     fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
244     fHistCheckDecChan->Sumw2();
245     fHistCheckDecChan->SetMinimum(0);
246     fOutput->Add(fHistCheckDecChan);
247
248     fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
249     fHistCheckDecChanAcc->Sumw2();
250     fHistCheckDecChanAcc->SetMinimum(0);
251     fOutput->Add(fHistCheckDecChanAcc);
252
253     fPtVsYGen= new TH2F("hPtVsYGen","",20,0.,10.,20,-1.,1.);
254     fPtVsYGen->Sumw2();
255     fPtVsYGen->SetMinimum(0);
256     fOutput->Add(fPtVsYGen);
257
258     fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",20,0.,10.,20,-1.,1.);
259     fPtVsYGenLimAcc->Sumw2();
260     fPtVsYGenLimAcc->SetMinimum(0);
261     fOutput->Add(fPtVsYGenLimAcc);
262
263     fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",20,0.,10.,20,-1.,1.);
264     fPtVsYGenAcc->Sumw2();
265     fPtVsYGenAcc->SetMinimum(0);
266     fOutput->Add(fPtVsYGenAcc);
267
268     fPtVsYReco= new TH2F("hPtVsYReco","",20,0.,10.,20,-1.,1.);
269     fPtVsYReco->Sumw2();
270     fPtVsYReco->SetMinimum(0);
271     fOutput->Add(fPtVsYReco);
272   }
273
274
275  Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
276   Double_t maxm=fMinMass+nMassBins*0.001;
277   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
278   fMassVsPtVsY->Sumw2();
279   fMassVsPtVsY->SetMinimum(0);
280   fOutput->Add(fMassVsPtVsY);
281
282   fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
283   fMassVsPtVsYRot->Sumw2();
284   fMassVsPtVsYRot->SetMinimum(0);
285   fOutput->Add(fMassVsPtVsYRot);
286
287   if(fMeson==kDzero){
288     fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
289     fMassVsPtVsYLSpp->Sumw2();
290     fMassVsPtVsYLSpp->SetMinimum(0);
291     fOutput->Add(fMassVsPtVsYLSpp);
292     fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
293     fMassVsPtVsYLSmm->Sumw2();
294     fMassVsPtVsYLSmm->SetMinimum(0);
295     fOutput->Add(fMassVsPtVsYLSmm);
296   }
297
298   fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
299   fMassVsPtVsYSig->Sumw2();
300   fMassVsPtVsYSig->SetMinimum(0);
301   fOutput->Add(fMassVsPtVsYSig);
302
303   fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
304   fMassVsPtVsYRefl->Sumw2();
305   fMassVsPtVsYRefl->SetMinimum(0);
306   fOutput->Add(fMassVsPtVsYRefl);
307
308   fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,20,0.,10.,20,-1.,1.);
309   fMassVsPtVsYBkg->Sumw2();
310   fMassVsPtVsYBkg->SetMinimum(0);
311   fOutput->Add(fMassVsPtVsYBkg);
312
313   fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
314   fNSelected->Sumw2();
315   fNSelected->SetMinimum(0);
316   fOutput->Add(fNSelected);
317
318   fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
319   fNormRotated->Sumw2();
320   fNormRotated->SetMinimum(0);
321   fOutput->Add(fNormRotated);
322
323   fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
324   fDeltaMass->Sumw2();
325   fDeltaMass->SetMinimum(0);
326   fOutput->Add(fDeltaMass);
327
328   Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
329   Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
330   Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
331   fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot);
332   fOutput->Add(fDeltaMassFullAnalysis);
333
334   //Counter for Normalization
335   fCounter = new AliNormalizationCounter("NormalizationCounter");
336   fCounter->Init();
337
338   PostData(1,fOutput); 
339   PostData(2,fCounter);   
340 }
341
342 //________________________________________________________________________
343 void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
344   //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
345
346   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
347   if(!aod && AODEvent() && IsStandardAOD()) {
348     // In case there is an AOD handler writing a standard AOD, use the AOD 
349     // event in memory rather than the input (ESD) event.    
350     aod = dynamic_cast<AliAODEvent*> (AODEvent());
351   }
352   if(!aod){
353     printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
354     return;
355   }
356
357   // fix for temporary bug in ESDfilter 
358   // the AODs with null vertex pointer didn't pass the PhysSel
359   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
360   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
361   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
362   AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
363   fPidHF->SetPidResponse(pidResp);
364   
365  
366   fHistNEvents->Fill(0); // count event
367   // Post the data already here
368   PostData(1,fOutput);
369
370   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
371  
372   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
373   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
374   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
375   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
376   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
377   //  if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
378   if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7); 
379
380   if(!isEvSel)return;
381  
382   fHistNEvents->Fill(1);
383   
384   TClonesArray *arrayMC=0;
385   AliAODMCHeader *mcHeader=0;
386   if(fReadMC){
387     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
388     if(!arrayMC) {
389       printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
390       return;
391     }
392     
393     // load MC header
394     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
395     if(!mcHeader) {
396       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
397       return;
398     }
399     FillGenHistos(arrayMC);
400   }
401
402
403   Int_t ntracks=aod->GetNTracks();
404
405   // select and flag tracks
406   UChar_t* status = new UChar_t[ntracks];
407   for(Int_t iTr=0; iTr<ntracks; iTr++){
408     status[iTr]=0;
409     AliAODTrack* track=aod->GetTrack(iTr);
410     if(IsTrackSelected(track)) status[iTr]+=1;
411     if(IsKaon(track)) status[iTr]+=2;
412     if(IsPion(track)) status[iTr]+=4;
413     fHistTrackStatus->Fill(status[iTr]);
414   }
415
416   // build the combinatorics
417   Int_t nSelected=0;
418   Int_t nFiltered=0;
419   Double_t dummypos[3]={0.,0.,0.};
420   AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
421   AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
422   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
423   Double_t d02[2]={0.,0.};
424   Double_t d03[3]={0.,0.,0.};
425   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
426   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
427   UInt_t pdg0[2]={321,211};
428   UInt_t pdgp[3]={321,211,211};
429   //  UInt_t pdgs[3]={321,321,211};
430   Double_t tmpp[3];
431   Double_t px[3],py[3],pz[3];
432   Int_t dgLabels[3];
433
434   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
435     AliAODTrack* trK=aod->GetTrack(iTr1);
436     if((status[iTr1] & 1)==0) continue;
437     if((status[iTr1] & 2)==0) continue;
438     Int_t chargeK=trK->Charge();
439     trK->GetPxPyPz(tmpp);
440     px[0] = tmpp[0]; 
441     py[0] = tmpp[1]; 
442     pz[0] = tmpp[2]; 
443     dgLabels[0]=trK->GetLabel();
444     for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
445       if((status[iTr2] & 1)==0) continue;
446       if((status[iTr2] & 4)==0) continue;
447       if(iTr1==iTr2) continue;
448       AliAODTrack* trPi1=aod->GetTrack(iTr2);
449       Int_t chargePi1=trPi1->Charge();
450       trPi1->GetPxPyPz(tmpp);
451       px[1] = tmpp[0]; 
452       py[1] = tmpp[1]; 
453       pz[1] = tmpp[2]; 
454       dgLabels[1]=trPi1->GetLabel();
455       if(chargePi1==chargeK){
456         if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
457         continue;
458       }
459       if(fMeson==kDzero){
460         nFiltered++;
461         v2->AddDaughter(trK);
462         v2->AddDaughter(trPi1);
463         tmpRD2->SetSecondaryVtx(v2);
464         Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
465         v2->RemoveDaughters();
466         if(ok) nSelected++;
467       }else{
468         for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
469           if((status[iTr3] & 1)==0) continue;
470           if((status[iTr3] & 4)==0) continue;
471           if(iTr1==iTr3) continue;
472           AliAODTrack* trPi2=aod->GetTrack(iTr3);
473           Int_t chargePi2=trPi2->Charge();
474           if(chargePi2==chargeK) continue;
475           nFiltered++;
476           trPi2->GetPxPyPz(tmpp);
477           px[2] = tmpp[0]; 
478           py[2] = tmpp[1]; 
479           pz[2] = tmpp[2]; 
480           dgLabels[2]=trPi2->GetLabel();
481           v3->AddDaughter(trK);
482           v3->AddDaughter(trPi1);
483           v3->AddDaughter(trPi2);
484           tmpRD3->SetSecondaryVtx(v3);
485           Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
486           v3->RemoveDaughters();
487           if(ok) nSelected++;
488         }
489       }
490     }
491   }
492
493   delete [] status;
494   delete v2;
495   delete v3;
496   delete tmpRD2;
497   delete tmpRD3;
498
499   fNSelected->Fill(nSelected);
500
501   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
502   fCounter->StoreCandidates(aod,nSelected,kFALSE);
503
504   PostData(1,fOutput); 
505   PostData(2,fCounter);    
506
507   return;
508 }
509
510 //________________________________________________________________________
511 void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
512   // Fill histos for LS candidates
513     
514   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
515   Double_t pt = tmpRD->Pt();
516   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
517   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
518     Double_t rapid = tmpRD->Y(pdgD);
519     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
520       if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
521       else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
522     }
523   }
524   return;
525 }
526
527 //________________________________________________________________________
528 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
529   // Fill histos with generated quantities
530   Int_t totPart=arrayMC->GetEntriesFast();
531   Int_t thePDG=411;
532   Int_t nProng=3;
533   if(fMeson==kDzero){
534     thePDG=421;
535     nProng=2;
536   }
537   for(Int_t ip=0; ip<totPart; ip++){
538     AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
539     if(TMath::Abs(part->GetPdgCode())==thePDG){
540       Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
541       fHistCheckOrigin->Fill(orig);
542       if(fPromptFeeddown==kFeeddown && orig!=5) continue;
543       else if(fPromptFeeddown==kPrompt && orig!=4) continue;
544       else if(fPromptFeeddown==kBoth && orig<4) continue;
545       fHistCheckOriginSel->Fill(orig);
546       Int_t deca=0;
547       Bool_t isGoodDecay=kFALSE;
548       Int_t labDau[4]={-1,-1,-1,-1};
549       if(fMeson==kDzero){
550         deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
551         if(part->GetNDaughters()!=2) continue;
552         if(deca==1) isGoodDecay=kTRUE;
553       }else if(fMeson==kDplus){ 
554         deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
555         if(deca>0) isGoodDecay=kTRUE;
556       }
557       fHistCheckDecChan->Fill(deca);
558       if(labDau[0]==-1){
559         //      printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
560         continue; //protection against unfilled array of labels
561       }
562       Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
563       if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
564       if(isGoodDecay){
565         Double_t ptgen=part->Pt();
566         Double_t ygen=part->Y();
567         if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
568           fPtVsYGen->Fill(ptgen,ygen);
569           if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
570           if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
571         }
572       }
573     }
574   }
575 }
576
577 //________________________________________________________________________
578 Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels){
579   // Fill histos for candidates with proper charge sign
580
581   Bool_t accept=kFALSE;
582
583   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
584   Double_t pt = tmpRD->Pt();
585   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
586   Double_t mass=TMath::Sqrt(minv2);
587
588   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
589     Double_t rapid = tmpRD->Y(pdgD);
590     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
591       fMassVsPtVsY->Fill(mass,pt,rapid);
592       accept=kTRUE;
593       if(fReadMC){
594         Int_t signPdg[3]={0,0,0};
595         for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
596         Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
597         if(labD>=0){
598           AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
599           if(part){
600             Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
601             if(pdgCode==321){
602               AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
603               if(dmes){
604                 Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,dmes,fGoUpToQuark);
605                 if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
606                   fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
607                 }
608               }
609               fMassVsPtVsYSig->Fill(mass,pt,rapid);
610             }else{
611               fMassVsPtVsYRefl->Fill(mass,pt,rapid);
612             }
613           }
614         }else{
615           fMassVsPtVsYBkg->Fill(mass,pt,rapid);
616         }
617       }
618     }
619   }
620
621   Int_t nRotated=0;
622   Double_t massRot=0;// calculated later only if candidate is acceptable
623   Double_t angleProngXY;
624   if(TMath::Abs(pdgD)==421)angleProngXY=TMath::ACos((px[0]*px[1]+py[0]*py[1])/TMath::Sqrt((px[0]*px[0]+py[0]*py[0])*(px[1]*px[1]+py[1]*py[1])));
625   else {
626     angleProngXY=TMath::ACos(((px[0]+px[1])*px[2]+(py[0]+py[1])*py[2])/TMath::Sqrt(((px[0]+px[1])*(px[0]+px[1])+(py[0]+py[1])*(py[0]+py[1]))*(px[2]*px[2]+py[2]*py[2])));
627   }
628   Double_t ptOrig=pt;
629
630
631   Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
632   Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
633
634   for(Int_t irot=0; irot<fNRotations; irot++){
635     Double_t phirot=fMinAngleForRot+rotStep*irot;
636     Double_t tmpx=px[0];
637     Double_t tmpy=py[0];
638     Double_t tmpx2=px[2];
639     Double_t tmpy2=py[2];
640     px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
641     py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
642     if(pdgD==411 || pdgD==431){
643       Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
644       px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
645       py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
646     }
647     tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
648     pt = tmpRD->Pt();
649     minv2 = tmpRD->InvMass2(nProngs,pdgdau);
650     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
651       Double_t rapid = tmpRD->Y(pdgD);
652       if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
653         massRot=TMath::Sqrt(minv2);
654         fMassVsPtVsYRot->Fill(massRot,pt,rapid);
655         nRotated++;
656         fDeltaMass->Fill(massRot-mass);
657         if(fFullAnalysis){
658           Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
659           fDeltaMassFullAnalysis->Fill(pointRot);
660         }
661       }
662     }
663     px[0]=tmpx;
664     py[0]=tmpy;
665     if(pdgD==411 || pdgD==431){
666       px[2]=tmpx2;
667       py[2]=tmpy2;
668     }
669   }
670   fNormRotated->Fill(nRotated);
671
672   return accept;
673
674 }
675 //________________________________________________________________________
676 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
677   // track selection cuts
678
679   if(track->Charge()==0) return kFALSE;
680   if(track->GetID()<0&&!fKeepNegID)return kFALSE;
681   if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
682   if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;    
683   return kTRUE;
684 }
685 //________________________________________________________________________
686 Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
687   // kaon selection cuts
688
689   if(!fPidHF) return kTRUE;
690   Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);  
691
692   if(SelectAODTrack(track,fTrackCutsKaon)) {
693     if(isKaon>=1)    return kTRUE;
694     else if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;
695   }
696   return kFALSE;
697 }
698 //________________________________________________________________________
699 Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
700   // pion selection cuts
701
702   if(!fPidHF) return kTRUE;
703   Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
704   if(SelectAODTrack(track,fTrackCutsPion)) {
705     if(isPion>=1)    return kTRUE;
706     else if(isPion>=0 && track->P()>fmaxPforIDPion)return kTRUE;
707   }
708   return kFALSE;
709 }
710 //________________________________________________________________________
711 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
712   // AOD track selection
713
714   if(!cuts) return kTRUE;
715
716   AliESDtrack esdTrack(track);
717   // set the TPC cluster info
718   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
719   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
720   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
721   if(!cuts->IsSelected(&esdTrack)) return kFALSE; 
722
723   return kTRUE;  
724 }
725
726 //_________________________________________________________________
727 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
728   // check if the decay products are in the good eta and pt range
729   for (Int_t iProng = 0; iProng<nProng; iProng++){
730     AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
731     if(!mcPartDaughter) return kFALSE;
732     Double_t eta = mcPartDaughter->Eta();
733     Double_t pt = mcPartDaughter->Pt();
734     if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
735   }
736   return kTRUE;
737 }
738
739 //_________________________________________________________________
740 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
741 {
742   // Terminate analysis
743   //
744   if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
745   fOutput = dynamic_cast<TList*> (GetOutputData(1));
746   if (!fOutput) {     
747     printf("ERROR: fOutput not available\n");
748     return;
749   }
750   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
751   if(fHistNEvents){
752     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
753   }else{
754     printf("ERROR: fHistNEvents not available\n");
755     return;
756   }
757   return;
758 }
759