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