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