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