Compatibility with the Root trunk
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSESignificance.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////
19 //
20 // AliAnalysisTaskSESignificane to calculate effects on 
21 // significance of D mesons  cut 
22 // Authors: G. Ortona, ortona@to.infn.it
23 // F. Prino, prino@to.infn.it
24 // Renu Bala, bala@to.infn.it
25 // Chiara Bianchin, cbianchi@pd.infn.it
26 /////////////////////////////////////////////////////////////
27
28 #include <Riostream.h>
29 #include <TClonesArray.h>
30 #include <TCanvas.h>
31 #include <TList.h>
32 #include <TFile.h>
33 #include <TString.h>
34 #include <TH1F.h>
35 #include <TDatabasePDG.h>
36
37 #include <AliLog.h>
38 #include "AliAnalysisManager.h"
39 #include "AliAODHandler.h"
40 #include "AliAODEvent.h"
41 #include "AliAODVertex.h"
42 #include "AliAODTrack.h"
43 #include "AliAODMCHeader.h"
44 #include "AliAODMCParticle.h"
45 #include "AliAODRecoDecayHF3Prong.h"
46 #include "AliAODRecoDecayHF.h"
47 #include "AliAODRecoDecayHF2Prong.h"
48 #include "AliAODRecoDecayHF4Prong.h"
49 #include "AliAODRecoCascadeHF.h"
50
51 #include "AliAnalysisTaskSE.h"
52 #include "AliRDHFCutsDplustoKpipi.h"
53 #include "AliRDHFCutsD0toKpipipi.h"
54 #include "AliRDHFCutsDstoKKpi.h"
55 #include "AliRDHFCutsDStartoKpipi.h"
56 #include "AliRDHFCutsD0toKpi.h"
57 #include "AliRDHFCutsLctopKpi.h"
58 #include "AliMultiDimVector.h"
59
60 #include "AliAnalysisTaskSESignificance.h"
61
62 using std::cout;
63 using std::endl;
64
65 ClassImp(AliAnalysisTaskSESignificance)
66
67
68 //________________________________________________________________________
69 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
70   AliAnalysisTaskSE(),
71   fOutput(0),
72   fCutList(0),
73   fHistNEvents(0),
74   fUpmasslimit(1.965),
75   fLowmasslimit(1.765),
76   fRDCuts(0),
77   fNPtBins(0),
78   fReadMC(kFALSE),
79   fUseSelBit(kFALSE),
80   fBFeedDown(kBoth),
81   fDecChannel(0),
82   fPDGmother(0),
83   fNProngs(0),
84   fBranchName(""),
85   fSelectionlevel(0),
86   fNVars(0),
87   fNBins(100),
88   fPartOrAndAntiPart(0),
89   fDsChannel(0)
90 {
91   // Default constructor
92   SetPDGCodes();
93   SetDsChannel(kPhi);
94
95   for(Int_t i=0;i<4;i++) fPDGdaughters[i]=0;
96   for(Int_t i=0;i<2;i++) {fPDGD0ToKpi[i]=0;fPDGDStarToD0pi[i]=0;}
97   for(Int_t i=0;i<kMaxCutVar;i++) fVars[i]=0.;
98   for(Int_t i=0;i<kMaxNHist;i++) {
99     fMassHist[i]=0;
100     fSigHist[i]=0;
101     fBkgHist[i]=0;
102     fRflHist[i]=0;
103   }
104
105 }
106
107 //________________________________________________________________________
108 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
109   AliAnalysisTaskSE(name),
110   fOutput(0),
111   fCutList(listMDV),
112   fHistNEvents(0),
113   fUpmasslimit(0),
114   fLowmasslimit(0),
115   fRDCuts(rdCuts),
116   fNPtBins(0),
117   fReadMC(kFALSE),
118   fUseSelBit(kFALSE),
119   fBFeedDown(kBoth),
120   fDecChannel(decaychannel),
121   fPDGmother(0),
122   fNProngs(0),
123   fBranchName(""),
124   fSelectionlevel(selectionlevel),
125   fNVars(0),
126   fNBins(100),
127   fPartOrAndAntiPart(0),
128   fDsChannel(0)
129 {
130
131   SetPDGCodes();
132   SetDsChannel(kPhi);
133   if (fDecChannel!=2) SetMassLimits(0.15,fPDGmother); //check range
134   else {
135     Float_t min = 0.13;
136     Float_t max = 0.19;
137     SetMassLimits(min, max);
138   }
139   fNPtBins=fRDCuts->GetNPtBins();
140
141   fNVars=fRDCuts->GetNVarsForOpt();
142   if(fNVars>kMaxCutVar) AliFatal(Form("Too large number of cut variables, maximum is %d",kMaxCutVar));
143   
144   if(fDebug>1)fRDCuts->PrintAll();
145    // Output slot #1 writes into a TList container
146   DefineOutput(1,TList::Class());  //My private output
147   DefineOutput(2,TList::Class());
148   DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
149   CheckConsistency();
150 }
151
152  //________________________________________________________________________
153 AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
154 {
155   // Destructor
156   if (fOutput) {
157     delete fOutput;
158     fOutput = 0;
159   }
160   if (fCutList) {
161     delete fCutList;
162     fCutList = 0;
163   }
164   if(fHistNEvents){
165     delete fHistNEvents;
166     fHistNEvents=0;
167   }
168 /*
169   if(fRDCuts) {
170     delete fRDCuts;
171     fRDCuts= 0;
172   } 
173 */
174   
175 }
176 //_________________________________________________________________
177 void AliAnalysisTaskSESignificance::SetPDGCodes(){
178   // sets channel dependent quantities
179  
180   switch(fDecChannel){
181   case 0:
182     //D+
183     fPDGmother=411;
184     fNProngs=3;
185     fPDGdaughters[0]=211;//pi
186     fPDGdaughters[1]=321;//K
187     fPDGdaughters[2]=211;//pi
188     fPDGdaughters[3]=0; //empty
189     fBranchName="Charm3Prong";
190     break;
191   case 1:
192     //D0
193     fPDGmother=421;
194     fNProngs=2;
195     fPDGdaughters[0]=211;//pi 
196     fPDGdaughters[1]=321;//K
197     fPDGdaughters[2]=0; //empty
198     fPDGdaughters[3]=0; //empty
199     fBranchName="D0toKpi";
200     break;
201   case 2:
202     //D*
203     fPDGmother=413;
204     fNProngs=3;
205     fPDGdaughters[1]=211;//pi
206     fPDGdaughters[0]=321;//K
207     fPDGdaughters[2]=211;//pi (soft?)
208     fPDGdaughters[3]=0; //empty
209     fBranchName="Dstar";
210     break;
211   case 3:
212     //Ds
213     fPDGmother=431;
214     fNProngs=3;
215     fPDGdaughters[0]=321;//K
216     fPDGdaughters[1]=321;//K
217     fPDGdaughters[2]=211;//pi
218     fPDGdaughters[3]=0; //empty
219     fBranchName="Charm3Prong";
220     break;
221   case 4:
222     //D0 in 4 prongs
223     fPDGmother=421;
224     fNProngs=4;
225     fPDGdaughters[0]=321;
226     fPDGdaughters[1]=211;
227     fPDGdaughters[2]=211;
228     fPDGdaughters[3]=211;
229     fBranchName="Charm4Prong";
230     break;
231   case 5:
232     //Lambda_c
233     fPDGmother=4122;
234     fNProngs=3;
235     fPDGdaughters[0]=2212;//p
236     fPDGdaughters[1]=321;//K
237     fPDGdaughters[2]=211;//pi
238     fPDGdaughters[3]=0; //empty
239     fBranchName="Charm3Prong";
240     break;
241   }
242 }
243
244 //_________________________________________________________________
245 Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
246
247   Bool_t result = kTRUE;
248
249   const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
250   //Float_t *vars = new Float_t[nvars];
251   Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
252   Bool_t *uppervars = fRDCuts->GetIsUpperCut();
253   TString *names = fRDCuts->GetVarNames();
254
255   for(Int_t i=0;i<fNPtBins;i++){
256     TString mdvname=Form("multiDimVectorPtBin%d",i);
257     Int_t ic=0;
258  
259     for(Int_t ivar=0;ivar<nvars;ivar++){
260       if(varsforopt[ivar]){
261         Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
262         Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
263         if(min==max){
264           AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
265           return kFALSE;
266         }
267         Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
268         if(uppervars[ivar]&&lowermdv){
269           AliFatal(Form("%s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---please check your cuts \n ",names[ivar].Data(),i));
270           return kFALSE;
271         }
272         if(!uppervars[ivar]&&!lowermdv){
273           AliFatal(Form("%s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---please check your cuts \n",names[ivar].Data(),i));
274           return kFALSE;
275         }
276         ic++;
277       }
278     }
279   }
280   return result;
281 }
282 //_________________________________________________________________
283 void AliAnalysisTaskSESignificance::SetBFeedDown(FeedDownEnum flagB){
284   if(fReadMC)fBFeedDown=flagB;
285   else {
286     AliInfo("B feed down not allowed without MC info\n");
287     fBFeedDown=kBoth;
288   }
289 }
290 //_________________________________________________________________
291 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
292   Float_t mass=0;
293   Int_t abspdg=TMath::Abs(pdg);
294   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
295   fUpmasslimit = mass+range;
296   fLowmasslimit = mass-range;
297 }
298 //_________________________________________________________________
299 void  AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
300   if(uplimit>lowlimit)
301     {
302       fUpmasslimit = uplimit;
303       fLowmasslimit = lowlimit;
304     }
305 }
306
307
308
309 //________________________________________________________________________
310 void AliAnalysisTaskSESignificance::LocalInit()
311 {
312   // Initialization
313
314   if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
315
316   switch(fDecChannel){
317   case 0:
318     {
319       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
320       // Post the data
321       PostData(3,copycut);
322     }
323     break;
324   case 1:
325     {
326       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
327       // Post the data
328       PostData(3,copycut);
329     }
330     break;
331   case 2:
332     {
333       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
334       // Post the data
335       PostData(3,copycut);
336     }
337     break;
338   case 3:
339     {
340       AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
341       // Post the data
342       PostData(3,copycut);
343     }
344     break;
345   case 4:
346     {
347       AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fRDCuts)));
348       // Post the data
349       PostData(3,copycut);
350     }
351     break;
352   case 5:
353     {
354       AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fRDCuts)));
355       // Post the data
356       PostData(3,copycut);
357     }
358     break;
359
360   default:
361     return;
362   }
363
364   TList *mdvList =  new TList();
365   mdvList->SetOwner();
366   mdvList = fCutList;
367   
368   PostData(2,mdvList);
369
370
371 }
372 //________________________________________________________________________
373 void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
374 {
375   // Create the output container
376  
377   if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
378
379   // Several histograms are more conveniently managed in a TList
380   fOutput = new TList();
381   fOutput->SetOwner();
382   fOutput->SetName("OutputHistos");
383
384   //same number of steps in each multiDimVectorPtBin%d !
385   Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
386   cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
387   nHist=nHist*fNPtBins;
388   cout<<"Total = "<<nHist<<endl;
389   for(Int_t i=0;i<nHist;i++){
390
391     TString hisname;
392     TString signame;
393     TString bkgname;
394     TString rflname;
395     TString title;
396     
397     hisname.Form("hMass_%d",i);
398     signame.Form("hSig_%d",i);
399     bkgname.Form("hBkg_%d",i);
400     rflname.Form("hRfl_%d",i);
401     
402     title.Form("Invariant mass;M[GeV/c^{2}];Entries");
403
404     fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
405     fMassHist[i]->Sumw2();
406     fOutput->Add(fMassHist[i]);
407
408     if(fReadMC){
409       fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
410       fSigHist[i]->Sumw2();
411       fOutput->Add(fSigHist[i]);
412       
413       fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
414       fBkgHist[i]->Sumw2();
415       fOutput->Add(fBkgHist[i]);
416
417       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi){
418         fRflHist[i]=new TH1F(rflname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
419         fRflHist[i]->Sumw2();
420         fOutput->Add(fRflHist[i]);
421       }
422     }
423   }
424
425   fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
426   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
427   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
428   fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
429   fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
430   fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
431   fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
432   if(fReadMC){
433     fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
434     fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
435   } else{
436     fHistNEvents->GetXaxis()->SetBinLabel(7,"N candidates");
437   }
438   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
439   fOutput->Add(fHistNEvents);
440
441  
442   PostData(1,fOutput);    
443   return;
444 }
445
446 //________________________________________________________________________
447 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
448 {
449   // Execute analysis for current event:
450   // heavy flavor candidates association to MC truth
451    
452   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
453   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
454   // Post the data already here
455   PostData(1,fOutput);
456   TClonesArray *arrayProng =0;
457
458   if(!aod && AODEvent() && IsStandardAOD()) { 
459     // In case there is an AOD handler writing a standard AOD, use the AOD 
460     // event in memory rather than the input (ESD) event.    
461     aod = dynamic_cast<AliAODEvent*> (AODEvent());
462     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
463     // have to taken from the AOD event hold by the AliAODExtension
464     AliAODHandler* aodHandler = (AliAODHandler*) 
465       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
466     if(aodHandler->GetExtensions()) {
467       
468       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
469       AliAODEvent *aodFromExt = ext->GetAOD();
470       arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
471       
472    }
473   } else if(aod) {
474     arrayProng=(TClonesArray*)aod->GetList()->FindObject(fBranchName.Data());
475   }
476   if(!aod || !arrayProng) {
477     AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
478     return;
479   }
480   
481   // fix for temporary bug in ESDfilter 
482   // the AODs with null vertex pointer didn't pass the PhysSel
483   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
484   TClonesArray *arrayMC=0;
485   AliAODMCHeader *mcHeader=0;
486
487   // load MC particles
488   if(fReadMC){
489     
490     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
491     if(!arrayMC) {
492       AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
493       return;
494     }
495     
496     // load MC header
497     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
498     if(!mcHeader) {
499       AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
500       return;
501     }
502   }
503
504   fHistNEvents->Fill(0); // count event
505
506   AliAODRecoDecayHF *d=0;
507   
508
509
510   Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
511   Int_t nProng = arrayProng->GetEntriesFast();
512   if(fDebug>1) printf("Number of D2H: %d\n",nProng);
513
514   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
515   TString trigclass=aod->GetFiredTriggerClasses();
516   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
517
518   if(fRDCuts->IsEventSelected(aod)) {
519     fHistNEvents->Fill(1);
520   }else{
521     if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
522       fHistNEvents->Fill(4);
523     return;
524   }
525   
526   for (Int_t iProng = 0; iProng < nProng; iProng++) {
527     fHistNEvents->Fill(6);
528     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
529
530     AliAODRecoCascadeHF* DStarToD0pi = NULL;
531     AliAODRecoDecayHF2Prong* D0Particle = NULL;
532     fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
533     fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
534
535     Bool_t isSelBit=kTRUE;
536     if(fUseSelBit){
537       if(fDecChannel==0) {
538         isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
539       }else{ 
540         if(fDecChannel==1) {
541           isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
542         }else{
543           if(fDecChannel==2) {
544             isSelBit=d->HasSelectionBit(AliRDHFCuts::kDstarCuts);
545           }else{
546             if(fDecChannel==3) {
547               isSelBit=d->HasSelectionBit(AliRDHFCuts::kDsCuts);
548             }else{
549               if(fDecChannel==5) isSelBit=d->HasSelectionBit(AliRDHFCuts::kLcCuts);
550             }
551           }
552         }
553       }
554     }
555     if(!isSelBit) continue; 
556
557     if (fDecChannel==2) {
558       DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
559       if (!DStarToD0pi->GetSecondaryVtx()) continue;
560       D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
561       if (!D0Particle) continue;
562     }
563     
564     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
565     Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
566
567     if(fReadMC && fBFeedDown!=kBoth && isSelected){
568       Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
569       if(labD>=0){
570         AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
571         Int_t pdgGranma = CheckOrigin(partD, arrayMC);
572         Int_t abspdgGranma = TMath::Abs(pdgGranma);
573         if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
574           //feed down particle
575           AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
576           fHistNEvents->Fill(7);
577           if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
578         }
579         else { 
580           //prompt particle
581           fHistNEvents->Fill(6);
582           if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
583         } 
584       }
585     }
586     
587     if(isSelected&&isFidAcc) {
588       fHistNEvents->Fill(2); // count selected with loosest cuts
589       if(fDebug>1) printf("+++++++Is Selected\n");
590     
591       Int_t nVals=0;
592       if(fDecChannel==3) SetPDGdaughterDstoKKpi();
593       fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
594       Int_t ptbin=fRDCuts->PtBin(d->Pt());
595       if(ptbin==-1) continue;
596       TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
597       AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
598
599       ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
600       if(fDebug>1)printf("nvals = %d\n",nVals);
601       for(Int_t ivals=0;ivals<nVals;ivals++){
602         if(addresses[ivals]>=muvec->GetNTotCells()){
603           if (fDebug>1) printf("Overflow!!\n");
604           delete [] addresses;
605           return;
606         }
607         
608         fHistNEvents->Fill(3);
609         
610         //fill the histograms with the appropriate method
611         switch (fDecChannel){
612         case 0:
613           FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
614           break;
615         case 1:
616           FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
617           break;
618         case 2:
619           FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
620           break;
621         case 3:
622           if(isSelected&1){
623             FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
624           }
625           break;
626         case 4:
627           FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
628           break;
629         case 5:
630           FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
631           break;
632         default:
633           break;
634         }
635         
636       }
637       
638       if (fDecChannel==3 && isSelected&2){
639         SetPDGdaughterDstopiKK();
640         nVals=0;
641         fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
642         delete [] addresses;
643         addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
644         if(fDebug>1)printf("nvals = %d\n",nVals);
645         for(Int_t ivals=0;ivals<nVals;ivals++){
646           if(addresses[ivals]>=muvec->GetNTotCells()){
647             if (fDebug>1) printf("Overflow!!\n");
648             delete [] addresses;            
649             return;
650           }
651           FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,0);      
652           
653         }
654
655       }
656       
657       delete [] addresses;
658     }// end if selected
659     
660   }
661   
662
663   PostData(1,fOutput);    
664   return;
665 }
666
667 //***************************************************************************
668
669 // Methods used in the UserExec
670
671
672 //********************************************************************************************
673
674 //Methods to fill the istograms with MC information, one for each candidate
675 //NB: the implementation for each candidate is responsibility of the corresponding developer
676
677 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
678     //D+ channel
679   if(!isSel){
680     AliError("Candidate not selected\n");
681     return;
682   }
683
684   if(fPartOrAndAntiPart*d->GetCharge()<0)return;
685   Int_t pdgdaughters[3] = {211,321,211};
686   Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
687
688   fMassHist[index]->Fill(mass);
689
690
691   if(fReadMC){
692     Int_t lab=-1;
693     lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
694     if(lab>=0){ //signal
695       fSigHist[index]->Fill(mass);
696     } else{ //background
697       fBkgHist[index]->Fill(mass);
698     } 
699   }   
700 }
701
702
703 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
704
705   //D0->Kpi channel
706
707   //mass histograms
708   Int_t pdgdaughtersD0[2]={211,321};//pi,K 
709   Int_t pdgdaughtersD0bar[2]={321,211};//K,pi 
710   Double_t masses[2];
711   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
712   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
713
714   if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
715   if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
716
717
718
719   //MC histograms
720   if(fReadMC){
721
722     Int_t matchtoMC=-1;
723
724     //D0
725     Int_t pdgdaughters[2];
726     pdgdaughters[0]=211;//pi 
727     pdgdaughters[1]=321;//K
728
729     matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
730
731     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
732     if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
733       if(matchtoMC>=0){
734         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
735         Int_t pdgMC = dMC->GetPdgCode();
736         
737         if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
738         else fRflHist[index]->Fill(masses[0]);
739         
740       } else fBkgHist[index]->Fill(masses[0]);
741       
742     }
743     if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
744       if(matchtoMC>=0){
745         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
746         Int_t pdgMC = dMC->GetPdgCode();
747         
748         if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
749         else fRflHist[index]->Fill(masses[1]);
750       } else fBkgHist[index]->Fill(masses[1]);
751     }
752   }
753 }
754
755 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
756     //D* channel
757
758     AliInfo("Dstar selected\n");
759     
760     Double_t mass = dstarD0pi->DeltaInvMass();
761
762     if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(mass);
763     if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(mass);
764         
765     if(fReadMC) {
766        Int_t matchtoMC = -1; 
767        matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
768
769        Int_t prongPdgDStarPlus=413,prongPdgDStarMinus=(-1)*413;
770         
771        if ((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) { 
772           //D*+
773           if(matchtoMC>=0) {
774              AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
775              Int_t pdgMC = dMC->GetPdgCode();
776         
777              if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
778              else {
779                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
780                 mass =  dstarD0pi->DeltaInvMass();
781                 fRflHist[index]->Fill(mass);
782                 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
783               } 
784             } 
785           else fBkgHist[index]->Fill(mass);
786         }
787     
788         if (isSel>=2 && fPartOrAndAntiPart<=0) { 
789            //D*-
790            if (matchtoMC>=0) {
791               AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
792               Int_t pdgMC = dMC->GetPdgCode();
793         
794               if (pdgMC==prongPdgDStarMinus) fSigHist[index]->Fill(mass);
795               else {
796                  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
797                  mass = dstarD0pi->DeltaInvMass();
798                  fRflHist[index]->Fill(mass);
799                  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
800               }
801             } 
802             else fBkgHist[index]->Fill(mass);
803           }
804        }
805
806 }
807
808
809 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
810
811   // Fill Ds histos
812
813
814   Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi 
815   Int_t pdgDspiKK[3]={211,321,321};//pi,K,K 
816   Double_t masses[2];  
817   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
818   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar 
819
820   Int_t labDs=-1;
821   if(fReadMC){
822     labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
823   }
824
825   Int_t isKKpi=isSel&1;
826   Int_t ispiKK=isSel&2;
827   Int_t isPhiKKpi=isSel&4;
828   Int_t isPhipiKK=isSel&8;
829   Int_t isK0starKKpi=isSel&16;
830   Int_t isK0starpiKK=isSel&32;
831
832
833   if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
834   if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
835    
836   if (optDecay==1){ 
837     if(isKKpi  && fPartOrAndAntiPart*d->GetCharge()>=0) {
838       if(fDsChannel==kPhi && isPhiKKpi==0) return;
839       if(fDsChannel==kK0star && isK0starKKpi==0) return;
840       
841       fMassHist[index]->Fill(masses[0]); 
842       
843       if(fReadMC){
844         if(labDs>=0){
845           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
846           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
847           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
848           if(pdgCode0==321) {
849             fSigHist[index]->Fill(masses[0]); //signal
850           }else{
851             fRflHist[index]->Fill(masses[0]); //Reflected signal
852           }
853         }else{
854           fBkgHist[index]->Fill(masses[0]); // Background
855         }
856       }
857     }
858   }
859   
860   if (optDecay==0){ 
861     if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
862       if(fDsChannel==kPhi && isPhipiKK==0) return;
863       if(fDsChannel==kK0star && isK0starpiKK==0) return;
864       
865       fMassHist[index]->Fill(masses[1]);
866       
867       if(fReadMC){
868         if(labDs>=0){
869           Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
870           AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
871           Int_t pdgCode0=TMath::Abs(p->GetPdgCode());   
872           if(pdgCode0==211) {             
873             fSigHist[index]->Fill(masses[1]);
874           }else{
875             fRflHist[index]->Fill(masses[1]);
876           }
877         }else{
878           fBkgHist[index]->Fill(masses[1]);
879         }
880       }
881     }
882   }
883 }
884
885 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
886   //D0->Kpipipi channel
887   AliInfo("D0 in 4 prongs channel not implemented\n");
888
889 }
890
891 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
892
893   // Mass hypothesis
894   Double_t masses[2];
895   Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
896   masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
897   Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
898   masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
899
900
901   if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
902     
903     // isSel=1 : p K pi ; isSel=2 : pi K p ;
904     if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
905     if(isSel>=2) fMassHist[index]->Fill(masses[1]);
906     
907     // Check the MC truth
908     if(fReadMC){
909
910       Int_t ispKpi = 0;
911       Int_t ispiKp = 0;
912       ispKpi = isSel&1;
913       ispiKp = isSel&2;
914       Int_t matchtoMC = -1;
915       //
916       Int_t pPDG = 2212; // p 
917       Int_t kPDG = 321;  // K
918       Int_t piPDG = 211; // pi
919       Int_t absPdgMom = 4122;
920       matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
921
922       if(matchtoMC>=0){
923
924         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
925         Int_t pdgMC = dMC->GetPdgCode();
926         if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
927         Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
928         AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
929         Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
930         Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
931         AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
932         Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
933         Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
934         AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
935         Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
936
937         // Fill in the histograms in case of p K pi decays
938         if(ispKpi==1){
939           if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
940             fSigHist[index]->Fill(masses[0]);
941           } else {
942             fRflHist[index]->Fill(masses[0]);
943           }
944         }
945         // Fill in the histograms in case of pi K p decays
946         if(ispiKp==2){
947           if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
948             fSigHist[index]->Fill(masses[1]);
949           } else {
950             fRflHist[index]->Fill(masses[1]);
951           }
952         }
953       } else {
954         if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
955         if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
956       }
957     }
958   }
959
960
961 }
962
963
964 //________________________________________________________________________
965 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
966 {
967   // Terminate analysis
968   //
969   if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
970
971
972   fOutput = dynamic_cast<TList*> (GetOutputData(1));
973   if (!fOutput) {     
974     printf("ERROR: fOutput not available\n");
975     return;
976   }
977  
978   fCutList = dynamic_cast<TList*> (GetOutputData(2));
979   if (!fCutList) {     
980     printf("ERROR: fCutList not available\n");
981     return;
982   }
983
984   AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
985   if (!mdvtmp){
986     cout<<"multidimvec not found in TList"<<endl;
987     fCutList->ls();
988     return;
989   }
990   Int_t nHist=mdvtmp->GetNTotCells();
991   TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
992   Bool_t drawn=kFALSE;
993   for(Int_t i=0;i<nHist;i++){
994
995     TString hisname;
996     TString signame;
997     TString bkgname;
998     TString rflname;
999     
1000     hisname.Form("hMass_%d",i);
1001     signame.Form("hSig_%d",i);
1002     bkgname.Form("hBkg_%d",i);
1003     rflname.Form("hRfl_%d",i);
1004
1005     fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1006
1007     if (!drawn && fMassHist[i]->GetEntries() > 0){
1008       c1->cd();
1009       fMassHist[i]->Draw();
1010       drawn=kTRUE;
1011     }
1012     
1013     if(fReadMC){
1014       fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
1015       fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
1016       if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1017     }
1018     
1019   }
1020   
1021
1022  
1023   
1024   return;
1025 }
1026 //_________________________________________________________________________________________________
1027 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1028
1029         //
1030         // checking whether the very mother of the D0 is a charm or a bottom quark
1031         //
1032
1033         Int_t pdgGranma = 0;
1034         Int_t mother = 0;
1035         mother = mcPart->GetMother();
1036         Int_t istep = 0;
1037         while (mother >0 ){
1038                 istep++;
1039                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1040                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1041                 if(!mcGranma) break;
1042                 pdgGranma = mcGranma->GetPdgCode();
1043                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1044                 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1045                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1046                         break;
1047                 }
1048                 mother = mcGranma->GetMother();
1049         }
1050         return pdgGranma;
1051 }
1052 //_________________________________________________________________________________________________