cf0a0e61449098395fdc8159e1bcbe56b51fb14a
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx
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 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A.Grelli (Utrecht University) a.grelli@uu.nl
23 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
24 //-----------------------------------------------------------------------
25
26 #include <TDatabasePDG.h>
27 #include <TParticle.h>
28 #include <TVector3.h>
29 #include "TROOT.h"
30 #include <TH3F.h>
31
32 #include "AliRDHFCutsDStartoKpipi.h"
33 #include "AliRDHFCutsD0toKpi.h"
34 #include "AliAODMCHeader.h"
35 #include "AliAODHandler.h"
36 #include "AliAnalysisManager.h"
37 #include "AliLog.h"
38 #include "AliAODVertex.h"
39 #include "AliAODRecoDecay.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAODRecoDecayHF2Prong.h"
42 #include "AliESDtrack.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
45
46 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
47
48 //_______________________________________________________________________________
49
50 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
51 AliAnalysisTaskSE(),
52 fUseMCInfo(kFALSE),
53 fUseReco(kTRUE),
54 fCandidateType(0),
55 fCandidateName(""),
56 fPDGmother(0),
57 fNProngs(0),
58 fBranchName(""),
59 fOutput(0),
60 fCuts(0),
61 fMinMass(0.),
62 fMaxMass(0.),
63 fCandidateArray(0)
64
65 {
66    //
67    // Default constructor
68    //
69    
70    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
71    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
72 }
73
74 //_______________________________________________________________________________
75
76 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
77 AliAnalysisTaskSE(name),
78 fUseMCInfo(kFALSE),
79 fUseReco(kTRUE),
80 fCandidateType(candtype),
81 fCandidateName(""),
82 fPDGmother(0),
83 fNProngs(0),
84 fBranchName(""),
85 fOutput(0),
86 fCuts(cuts),
87 fMinMass(0.),
88 fMaxMass(0.),
89 fCandidateArray(0)
90 {
91    //
92    // Constructor. Initialization of Inputs and Outputs
93    //
94    
95    Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
96    
97    for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
98    for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
99    
100    const Int_t nptbins = fCuts->GetNPtBins();
101    Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
102    
103    switch (fCandidateType) {
104    case 0 :
105       fCandidateName = "D0";
106       fPDGmother = 421;
107       fNProngs = 2;
108       fPDGdaughters[0] = 211;  // pi 
109       fPDGdaughters[1] = 321;  // K
110       fPDGdaughters[2] = 0;    // empty
111       fPDGdaughters[3] = 0;    // empty
112       fBranchName = "D0toKpi";
113       break;
114    case 1 :
115       fCandidateName = "Dstar";
116       fPDGmother = 413;
117       fNProngs = 3;
118       fPDGdaughters[1] = 211; // pi soft
119       fPDGdaughters[0] = 421; // D0
120       fPDGdaughters[2] = 211; // pi fromD0
121       fPDGdaughters[3] = 321; // K from D0
122       fBranchName = "Dstar";
123       
124       if (nptbins<=13) {
125          for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
126       } else {
127          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
128       }
129       break;
130    default :
131       printf("%d not accepted!!\n",fCandidateType);
132       break;
133    }
134    
135    if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
136    if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
137    
138    DefineOutput(1, TList::Class());       // histos
139    DefineOutput(2, AliRDHFCuts::Class()); // my cuts
140 }
141
142 //_______________________________________________________________________________
143
144 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
145 {
146    //
147    // destructor
148    //
149    
150    Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  
151    
152    if (fOutput) { delete fOutput; fOutput = 0; }
153    if (fCuts)   { delete fCuts;   fCuts   = 0; }
154    if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
155    
156 }
157
158 //_______________________________________________________________________________
159
160 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
161 {
162    //
163    // Initialization
164    //
165    
166    if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
167    
168    switch (fCandidateType) {
169    case 0: {
170          AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
171          copyfCutsDzero->SetName("AnalysisCutsDzero");
172          PostData(2, copyfCutsDzero);  // Post the data
173    } break;
174 case 1: {
175       AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
176       copyfCutsDstar->SetName("AnalysisCutsDStar");
177       PostData(2, copyfCutsDstar); // Post the cuts
178 } break;
179 default: return;
180    }
181    
182    return;
183 }
184
185 //_______________________________________________________________________________
186
187 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
188
189    //
190    // output 
191    //
192    
193    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
194    
195    fOutput = new TList(); fOutput->SetOwner();
196    DefineHistoForAnalysis(); // define histograms
197    
198    fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
199    fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
200    
201    if (fCandidateType==kDstartoKpipi){
202       fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!
203       fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
204    }
205    
206    PostData(1, fOutput);
207    return;
208 }
209
210 //_______________________________________________________________________________
211
212 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
213    //
214    // user exec
215    //
216    
217    // add cadidate branch
218    fCandidateArray->Delete();
219    if (!(InputEvent()->FindListObject(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fCandidateArray);
220    if (fCandidateType==kDstartoKpipi){
221       fSideBandArray->Delete();
222       if (!(InputEvent()->FindListObject(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen")))) InputEvent()->AddObject(fSideBandArray);
223    }
224    //Printf("Arr names %s, %s",fCandidateArray->GetName(),fSideBandArray->GetName());
225    // Load the event
226    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
227    
228    TClonesArray *arrayDStartoD0pi = 0;
229    if (!aodEvent && AODEvent() && IsStandardAOD()) {
230       
231       // In case there is an AOD handler writing a standard AOD, use the AOD 
232       // event in memory rather than the input (ESD) event.    
233       aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
234       
235       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
236       // have to taken from the AOD event hold by the AliAODExtension
237       AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
238       if(aodHandler->GetExtensions()) {
239          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
240          AliAODEvent *aodFromExt = ext->GetAOD();
241          arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
242       }
243    } else {
244       arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
245    }
246    
247    if (!arrayDStartoD0pi) {
248       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
249       return;
250    } else {
251       AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
252    }
253    
254    TClonesArray* mcArray = 0x0;
255    if (fUseMCInfo) {
256       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
257       if (!mcArray) {
258          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
259          return;
260       }
261    }
262    
263    //Histograms
264    TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
265    TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
266    TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
267    TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
268    
269    TH1F* hPtPion=0x0;
270    if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
271    hstat->Fill(0);
272    
273    // fix for temporary bug in ESDfilter 
274    // the AODs with null vertex pointer didn't pass the PhysSel
275    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
276    
277    //Event selection
278    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
279    //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
280    if (!iseventselected) return;
281    hstat->Fill(1);
282    
283    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
284    if(fUseReco) hstat->Fill(2, nD);
285    
286    //D* and D0 prongs needed to MatchToMC method
287    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
288    Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
289    
290    //D0 from D0 bar
291    Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K 
292    Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi 
293    
294    Int_t iCand =0;
295    Int_t iSBCand=0;
296    Int_t isSelected = 0;
297    AliAODRecoDecayHF *charmCand = 0;
298    AliAODMCParticle *charmPart = 0;
299    Bool_t isMCBkg=kFALSE;
300    
301    Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
302    
303    Int_t mcLabel = -9999;
304    Int_t pdgMeson = 413;
305    if (fCandidateType==kD0toKpi) pdgMeson = 421;
306    
307    for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
308       charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
309       if (!charmCand) continue;
310       
311       
312       if (fUseMCInfo) { // Look in MC, try to simulate the z
313          if (fCandidateType==kDstartoKpipi) {
314             AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
315             mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
316          }
317          
318          if (fCandidateType==kD0toKpi) 
319             mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
320          
321          if (mcLabel<=0) isMCBkg=kTRUE;
322          else hstat->Fill(2);
323          if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
324       }
325       
326       Double_t ptD = charmCand->Pt();
327       
328       // region of interest + cuts
329       if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;    
330       
331       if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
332          //select by track cuts the side band candidates (don't want mass cut)
333          isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent); 
334          if (!isSelected) continue;
335          //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
336          Int_t bin = fCuts->PtBin(ptD);
337          if(bin<0 || bin>=fCuts->GetNPtBins()) {
338             AliError(Form("Pt %.3f out of bounds",ptD));
339             continue;
340          }
341          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
342          //if data and Dstar from D0 side band
343          if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){   
344             
345             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);
346             iSBCand++;
347          }
348       }
349       //candidate selected by cuts and PID
350       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
351       if (!isSelected) continue;
352       
353       //for data and MC signal fill fCandidateArray
354       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
355          // for data or MC with the requirement fUseReco fill with candidates
356          if(fUseReco) {
357             new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
358             //Printf("Filling reco");
359             hstat->Fill(3);
360          }
361          // for MC with requirement particle level fill with AliAODMCParticle
362          else if (fUseMCInfo) {
363             new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
364             //Printf("Filling gen");
365             hstat->Fill(3);
366          }
367          
368          iCand++;       
369       }
370       //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
371       else if(fUseReco){
372          new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
373          iSBCand++;
374       }
375       
376       
377       Double_t masses[2];
378       if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
379          
380          //softpion from D* decay
381          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
382          AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();  
383          
384          // select D* in the D0 window.
385          // In the cut object window is loose to allow for side bands
386          
387          
388          // retrieve the corresponding pt bin for the candidate
389          // and set the expected D0 width (x3)
390          //    static const Int_t n = fCuts->GetNPtBins();
391          Int_t bin = fCuts->PtBin(ptD);
392          if(bin<0 || bin>=fCuts->GetNPtBins()) {
393             AliError(Form("Pt %.3f out of bounds",ptD));
394             continue;
395          }
396          
397          AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
398          //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
399          if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {        
400             masses[0] = temp->DeltaInvMass(); //D*
401             masses[1] = 0.; //dummy for D*
402             
403             //D*  delta mass
404             hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
405             
406             // D* pt and soft pion pt for good candidates               
407             Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
408             Double_t invmassDelta = temp->DeltaInvMass();
409             if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
410          }
411       } //Dstar specific
412       
413       if (fCandidateType==kD0toKpi) { //D0->Kpi
414          
415          //needed quantities
416          masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
417          masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
418          hstat->Fill(3);
419          
420          // mass vs pt
421          if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
422          if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
423       } //D0 specific
424       
425       charmCand = 0;
426       isMCBkg=kFALSE;
427    } // end of D cand loop
428    
429    hnCandEv->Fill(fCandidateArray->GetEntriesFast());
430    if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
431       Int_t nsbcand=fSideBandArray->GetEntriesFast();
432       hstat->Fill(4,nsbcand);
433       hnSBCandEv->Fill(nsbcand);
434    }
435    
436    return;
437 }
438
439 //_______________________________________________________________________________
440
441 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
442 {
443    // The Terminate() function is the last function to be called during
444    // a query. It always runs on the client, it can be used to present
445    // the results graphically or save the results to file.
446    
447    Info("Terminate"," terminate");
448    AliAnalysisTaskSE::Terminate();
449    
450    fOutput = dynamic_cast<TList*>(GetOutputData(1));
451    if (!fOutput) {
452       printf("ERROR: fOutput not available\n");
453       return;
454    }
455    
456    return;
457 }
458
459 //_______________________________________________________________________________
460
461 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
462 {
463    //
464    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
465    //
466    
467    Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
468    
469    // compute the Delta mass for the D*
470    if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
471    
472    
473    fMinMass = mass - range;
474    fMaxMass = mass + range;
475    
476    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
477    if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
478    
479    return;
480 }
481
482 //_______________________________________________________________________________
483
484 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
485 {
486    //
487    // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
488    //
489    
490    if (uplimit>lowlimit) {
491       fMinMass = lowlimit;
492       fMaxMass = uplimit;
493    } else {
494       printf("Error! Lower limit larger than upper limit!\n");
495       fMinMass = uplimit - uplimit*0.2;
496       fMaxMass = uplimit;
497    }
498    
499    return;
500 }
501
502 //_______________________________________________________________________________
503
504 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
505 {
506    //
507    // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
508    //
509    
510    if (nptbins>30) {
511       AliInfo("Maximum number of bins allowed is 30!");
512       return kFALSE;
513    }
514    
515    if (!width) return kFALSE;
516    for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
517    
518    return kTRUE;
519 }
520
521 //_______________________________________________________________________________
522
523 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
524 {
525    //
526    // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
527    //
528    
529    // Statistics 
530    TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
531    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
532    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
533    if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
534    else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
535    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
536    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");  
537    hstat->SetNdivisions(1);
538    fOutput->Add(hstat);
539    
540    TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
541    fOutput->Add(hnCandEv);
542    
543    // Invariant mass related histograms
544    const Int_t nbinsmass = 200;
545    TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
546    hInvMass->SetStats(kTRUE);
547    hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
548    hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
549    fOutput->Add(hInvMass);
550    
551    if (fCandidateType==kDstartoKpipi) {
552       TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
553       fOutput->Add(hnSBCandEv);
554       
555       TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
556       hPtPion->SetStats(kTRUE);
557       hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
558       hPtPion->GetYaxis()->SetTitle("entries");
559       fOutput->Add(hPtPion);
560    }
561    
562    return kTRUE; 
563 }