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