4f9b192a836153e8f3429dbbf58dd772bd24e75e
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.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 // AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
19 // decay candidates with the MC truth.
20 // Authors: Renu Bala, bala@to.infn.it
21 // F. Prino, prino@to.infn.it
22 // G. Ortona, ortona@to.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TClonesArray.h>
26 #include <TNtuple.h>
27 #include <TCanvas.h>
28 #include <TList.h>
29 #include <TString.h>
30 #include <TH1F.h>
31 #include <TDatabasePDG.h>
32
33 #include "AliAnalysisManager.h"
34 #include "AliAODHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODVertex.h"
37 #include "AliAODTrack.h"
38 #include "AliAODMCHeader.h"
39 #include "AliAODMCParticle.h"
40 #include "AliAODRecoDecayHF3Prong.h"
41 #include "AliAnalysisVertexingHF.h"
42 #include "AliAnalysisTaskSE.h"
43 #include "AliAnalysisTaskSEDplus.h"
44
45 ClassImp(AliAnalysisTaskSEDplus)
46
47
48 //________________________________________________________________________
49 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
50 AliAnalysisTaskSE(),
51 fOutput(0), 
52 fHistNEvents(0),
53 fNtupleDplus(0),
54 fUpmasslimit(1.965),
55 fLowmasslimit(1.765),
56 fNPtBins(0),
57 fFillNtuple(kFALSE),
58 fReadMC(kFALSE),
59 fDoLS(kFALSE),
60 fVHF(0)
61 {
62    // Default constructor
63 }
64
65 //________________________________________________________________________
66 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
67 AliAnalysisTaskSE(name),
68 fOutput(0),
69 fHistNEvents(0),
70 fNtupleDplus(0),
71 fUpmasslimit(1.965),
72 fLowmasslimit(1.765),
73 fNPtBins(0),
74 fFillNtuple(fillNtuple),
75 fReadMC(kFALSE),
76 fDoLS(kFALSE),
77 fVHF(0)
78 {
79   Double_t ptlim[5]={0.,2.,3.,5,9999999.};
80   SetPtBinLimit(5, ptlim);
81   // Default constructor
82    // Output slot #1 writes into a TList container
83   DefineOutput(1,TList::Class());  //My private output
84
85   if(fFillNtuple){
86     // Output slot #2 writes into a TNtuple container
87     DefineOutput(2,TNtuple::Class());  //My private output
88   }
89 }
90
91 //________________________________________________________________________
92 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
93 {
94   // Destructor
95   if (fOutput) {
96     delete fOutput;
97     fOutput = 0;
98   }
99   if (fVHF) {
100     delete fVHF;
101     fVHF = 0;
102   }
103   
104   // if(fArrayBinLimits) {
105   //delete fArrayBinLimits;
106   //fArrayBinLimits= 0;
107   //} 
108   
109 }  
110 //_________________________________________________________________
111 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
112   fUpmasslimit = 1.865+range;
113   fLowmasslimit = 1.865-range;
114 }
115 //_________________________________________________________________
116 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
117   if(uplimit>lowlimit)
118     {
119       fUpmasslimit = lowlimit;
120       fLowmasslimit = uplimit;
121     }
122 }
123
124
125 //________________________________________________________________________
126 void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Double_t* lim){
127   // define pt bins for analysis
128   if(n>kMaxPtBins){
129     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
130     fNPtBins=kMaxPtBins;
131     fArrayBinLimits[0]=0.;
132     fArrayBinLimits[1]=2.;
133     fArrayBinLimits[2]=3.;
134     fArrayBinLimits[3]=5.;
135     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
136   }else{
137     fNPtBins=n-1;
138     fArrayBinLimits[0]=lim[0];
139     for(Int_t i=1; i<fNPtBins+1; i++) 
140       if(lim[i]>fArrayBinLimits[i-1]){
141         fArrayBinLimits[i]=lim[i];
142       }
143       else {
144         fArrayBinLimits[i]=fArrayBinLimits[i-1];
145       }
146     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
147   }
148   if(fDebug > 1){
149     printf("Number of Pt bins = %d\n",fNPtBins);
150     for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);    
151   }
152 }
153 //_________________________________________________________________
154 Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
155   if(ibin>fNPtBins)return -1;
156   return fArrayBinLimits[ibin];
157
158
159 //_________________________________________________________________
160 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
161
162 /*
163  * Fill the Like Sign histograms
164  */
165
166   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
167
168   //count pos/neg tracks
169   Int_t nPosTrks=0,nNegTrks=0;
170  //counter for particles passing single particle cuts
171   Int_t nspcplus=0;
172   Int_t nspcminus=0;
173
174   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
175     AliAODTrack *track = aod->GetTrack(it);
176     if(track->Charge()>0){
177       nPosTrks++;
178       if(track->Pt()>=0.4){
179         nspcplus++;
180       }
181     }
182     if(track->Charge()<0)
183       {
184         nNegTrks++;
185         if(track->Pt()>=0.4){
186           nspcminus++;
187         }
188       }
189   }
190
191   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
192
193   Int_t nDplusLS=0;
194   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
195   Int_t index; 
196
197   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
198     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
199     Bool_t unsetvtx=kFALSE;
200     if(!d->GetOwnPrimaryVtx()) {
201       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
202       unsetvtx=kTRUE;
203     }
204     if(d->SelectDplus(fVHF->GetDplusCuts()))nDplusLS++;
205     if(unsetvtx) d->UnsetOwnPrimaryVtx();
206   }
207
208  Float_t wei2=0;
209  if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
210  Float_t wei3=0;
211  if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
212
213  // loop over like sign candidates
214   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
215     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
216     Bool_t unsetvtx=kFALSE;
217     if(!d->GetOwnPrimaryVtx()) {
218       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
219       unsetvtx=kTRUE;
220     }
221  
222     if(d->SelectDplus(fVHF->GetDplusCuts())){
223
224       //set tight cuts values
225       Int_t iPtBin=-1;
226       Double_t ptCand = d->Pt();
227       if(ptCand<2.){
228         //iPtBin=0;
229         cutsDplus[7]=0.08;
230         cutsDplus[8]=0.5;
231         cutsDplus[9]=0.979;
232         cutsDplus[10]=0.0055;
233       }
234       else if(ptCand>2. && ptCand<3){ 
235         //iPtBin=1;
236         cutsDplus[7]=0.08;
237         cutsDplus[8]=0.5;
238         cutsDplus[9]=0.991;
239         cutsDplus[10]=0.005;
240       }else if(ptCand>3. && ptCand<5){ 
241         //iPtBin=2;
242         cutsDplus[7]=0.1;
243         cutsDplus[8]=0.5;
244         cutsDplus[9]=0.995;
245         cutsDplus[10]=0.0035;
246       }else{
247         //iPtBin=3;
248         cutsDplus[7]=0.1;
249         cutsDplus[8]=0.5;
250         cutsDplus[9]=0.997;
251         cutsDplus[10]=0.001;
252       }
253       
254       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
255         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
256       }
257       
258       if(iPtBin<0){
259         return;
260       }
261
262       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
263
264       Int_t sign= d->GetCharge();
265       Float_t wei=1;
266       Float_t wei4=1;
267       if(sign>0&&nPosTrks>2&&nspcplus>2) {  //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
268         
269         wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
270         wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
271       }
272       
273       if(sign<0&&nNegTrks>2&&nspcminus>2){     
274         wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
275         wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
276
277       }
278
279       Float_t invMass = d->InvMassDplus();
280       
281       
282       index=GetLSHistoIndex(iPtBin);
283       fMassHistLS[index]->Fill(invMass,wei);
284       fMassHistLS[index+1]->Fill(invMass);
285       fMassHistLS[index+2]->Fill(invMass,wei2);
286       fMassHistLS[index+3]->Fill(invMass,wei3);
287       fMassHistLS[index+4]->Fill(invMass,wei4);
288
289       
290       if(passTightCuts){
291         fMassHistLSTC[index]->Fill(invMass,wei);
292         fMassHistLSTC[index+1]->Fill(invMass);
293         fMassHistLSTC[index+2]->Fill(invMass,wei2);
294         fMassHistLSTC[index+3]->Fill(invMass,wei3);
295         fMassHistLSTC[index+4]->Fill(invMass,wei4);
296       }
297     }
298     if(unsetvtx) d->UnsetOwnPrimaryVtx();
299   }
300   
301   //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
302   //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
303
304   //  printf("LS analysis...done\n");
305
306 }
307
308 //________________________________________________________________________
309 void AliAnalysisTaskSEDplus::Init()
310 {
311   // Initialization
312
313   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
314
315   gROOT->LoadMacro("ConfigVertexingHF.C");
316
317   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
318   fVHF->PrintStatus();
319
320   return;
321 }
322
323 //________________________________________________________________________
324 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
325 {
326   // Create the output container
327   //
328   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
329
330   // Several histograms are more conveniently managed in a TList
331   fOutput = new TList();
332   fOutput->SetOwner();
333   fOutput->SetName("OutputHistos");
334
335   TString hisname;
336   Int_t index=0;
337   Int_t indexLS=0;
338   for(Int_t i=0;i<fNPtBins;i++){
339
340     index=GetHistoIndex(i);
341     indexLS=GetLSHistoIndex(i);
342     hisname.Form("hMassPt%d",i);
343     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
344     fMassHist[index]->Sumw2();
345     hisname.Form("hMassPt%dTC",i);
346     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
347     fMassHistTC[index]->Sumw2();
348     hisname.Form("hLSPt%dTC",i);
349     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
350     fMassHistLSTC[indexLS]->Sumw2();
351     hisname.Form("hLSPt%dLC",i);
352     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
353     fMassHistLS[indexLS]->Sumw2();
354     
355     index=GetSignalHistoIndex(i);    
356     indexLS++;
357     hisname.Form("hSigPt%d",i);
358     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
359     fMassHist[index]->Sumw2();
360     hisname.Form("hSigPt%dTC",i);
361     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
362     fMassHistTC[index]->Sumw2();
363     hisname.Form("hLSPt%dLCnw",i);
364     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
365     fMassHistLS[indexLS]->Sumw2();
366     hisname.Form("hLSPt%dTCnw",i);
367     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
368     fMassHistLSTC[indexLS]->Sumw2();
369
370     index=GetBackgroundHistoIndex(i); 
371     indexLS++;
372     hisname.Form("hBkgPt%d",i);
373     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
374     fMassHist[index]->Sumw2();
375     hisname.Form("hBkgPt%dTC",i);
376     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
377     fMassHistTC[index]->Sumw2();
378     hisname.Form("hLSPt%dLCntrip",i);
379     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
380     fMassHistLS[indexLS]->Sumw2();
381     hisname.Form("hLSPt%dTCntrip",i);
382     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
383     fMassHistLSTC[indexLS]->Sumw2();
384
385     indexLS++;
386     hisname.Form("hLSPt%dLCntripsinglecut",i);
387     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
388     fMassHistLS[indexLS]->Sumw2();
389     hisname.Form("hLSPt%dTCntripsinglecut",i);
390     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
391     fMassHistLSTC[indexLS]->Sumw2();
392
393     indexLS++;
394     hisname.Form("hLSPt%dLCspc",i);
395     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
396     fMassHistLS[indexLS]->Sumw2();
397     hisname.Form("hLSPt%dTCspc",i);
398     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),100,fLowmasslimit,fUpmasslimit);
399     fMassHistLSTC[indexLS]->Sumw2();
400   }
401   
402   for(Int_t i=0; i<3*fNPtBins; i++){
403     fOutput->Add(fMassHist[i]);
404     fOutput->Add(fMassHistTC[i]);
405   }
406   for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
407     fOutput->Add(fMassHistLS[i]);
408     fOutput->Add(fMassHistLSTC[i]);
409   }
410
411
412   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
413  fHistNEvents->Sumw2();
414  fHistNEvents->SetMinimum(0);
415  fOutput->Add(fHistNEvents);
416
417   if(fFillNtuple){
418     OpenFile(2); // 2 is the slot number of the ntuple
419     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2");
420   }
421
422   return;
423 }
424
425 //________________________________________________________________________
426 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
427 {
428   // Execute analysis for current event:
429   // heavy flavor candidates association to MC truth
430
431    AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
432  fHistNEvents->Fill(0); // count event
433   // Post the data already here
434   PostData(1,fOutput);
435   
436   TClonesArray *array3Prong = 0;
437   TClonesArray *arrayLikeSign =0;
438   if(!aod && AODEvent() && IsStandardAOD()) {
439     // In case there is an AOD handler writing a standard AOD, use the AOD 
440     // event in memory rather than the input (ESD) event.    
441     aod = dynamic_cast<AliAODEvent*> (AODEvent());
442     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
443     // have to taken from the AOD event hold by the AliAODExtension
444     AliAODHandler* aodHandler = (AliAODHandler*) 
445       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
446     if(aodHandler->GetExtensions()) {
447       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
448       AliAODEvent *aodFromExt = ext->GetAOD();
449       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
450       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
451     }
452   } else {
453     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
454     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
455   }
456
457   if(!array3Prong) {
458     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
459     return;
460   }
461   if(!arrayLikeSign) {
462     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
463     return;
464   }
465
466  
467   TClonesArray *arrayMC=0;
468   AliAODMCHeader *mcHeader=0;
469
470   // AOD primary vertex
471   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
472   //    vtx1->Print();
473   
474   // load MC particles
475   if(fReadMC){
476     
477     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
478     if(!arrayMC) {
479       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
480       //    return;
481     }
482     
483   // load MC header
484     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
485     if(!mcHeader) {
486     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
487     return;
488     }
489   }
490   
491   Int_t n3Prong = array3Prong->GetEntriesFast();
492   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
493   
494   
495   Int_t nOS=0;
496   Int_t index;
497   Int_t pdgDgDplustoKpipi[3]={321,211,211};
498   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
499   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
500     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
501     
502     
503     Bool_t unsetvtx=kFALSE;
504     if(!d->GetOwnPrimaryVtx()){
505       d->SetOwnPrimaryVtx(vtx1);
506       unsetvtx=kTRUE;
507     }
508
509     if(d->SelectDplus(fVHF->GetDplusCuts())) {
510       Int_t iPtBin = -1;
511       Double_t ptCand = d->Pt();
512       
513       if(ptCand<2.){
514         //iPtBin=0;
515         cutsDplus[7]=0.08;
516         cutsDplus[8]=0.5;
517         cutsDplus[9]=0.979;
518         cutsDplus[10]=0.0055;
519       }
520       else if(ptCand>2. && ptCand<3){ 
521         //iPtBin=1;
522         cutsDplus[7]=0.08;
523         cutsDplus[8]=0.5;
524         cutsDplus[9]=0.991;
525         cutsDplus[10]=0.005;
526       }else if(ptCand>3. && ptCand<5){ 
527         //iPtBin=2;
528         cutsDplus[7]=0.1;
529         cutsDplus[8]=0.5;
530         cutsDplus[9]=0.995;
531         cutsDplus[10]=0.0035;
532       }else{
533         //iPtBin=3;
534         cutsDplus[7]=0.1;
535         cutsDplus[8]=0.5;
536         cutsDplus[9]=0.997;
537         cutsDplus[10]=0.001;
538       }
539       
540       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
541         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
542       }
543       
544       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
545       Int_t labDp=-1;
546       Float_t deltaPx=0.;
547       Float_t deltaPy=0.;
548       Float_t deltaPz=0.;
549       Float_t truePt=0.;
550       Float_t xDecay=0.;
551       Float_t yDecay=0.;
552       Float_t zDecay=0.;
553       Float_t pdgCode=-2;
554       if(fReadMC){
555         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
556         if(labDp>=0){
557           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
558           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
559           deltaPx=partDp->Px()-d->Px();
560           deltaPy=partDp->Py()-d->Py();
561           deltaPz=partDp->Pz()-d->Pz();
562           truePt=partDp->Pt();
563           xDecay=dg0->Xv();       
564           yDecay=dg0->Yv();       
565           zDecay=dg0->Zv();
566           pdgCode=TMath::Abs(partDp->GetPdgCode());
567         }else{
568           pdgCode=-1;
569         }
570       }
571       Double_t invMass=d->InvMassDplus();
572
573       Float_t tmp[22];
574       if(fFillNtuple){            
575         tmp[0]=pdgCode;
576         tmp[1]=deltaPx;
577         tmp[2]=deltaPy;
578         tmp[3]=deltaPz;
579         tmp[4]=truePt;
580         tmp[5]=xDecay;    
581         tmp[6]=yDecay;    
582         tmp[7]=zDecay;    
583         tmp[8]=d->PtProng(0);
584         tmp[9]=d->PtProng(1);
585         tmp[10]=d->PtProng(2);
586         tmp[11]=d->Pt();
587         tmp[12]=d->CosPointingAngle();
588         tmp[13]=d->DecayLength();
589         tmp[14]=d->Xv();
590         tmp[15]=d->Yv();
591         tmp[16]=d->Zv();
592         tmp[17]=d->InvMassDplus();
593         tmp[18]=d->GetSigmaVert();
594         tmp[19]=d->Getd0Prong(0);
595         tmp[20]=d->Getd0Prong(1);
596         tmp[21]=d->Getd0Prong(2);         
597         fNtupleDplus->Fill(tmp);
598         PostData(2,fNtupleDplus);
599       }
600
601       if(iPtBin>=0){
602       
603         index=GetHistoIndex(iPtBin);
604         fMassHist[index]->Fill(invMass);
605         if(passTightCuts){
606           fMassHistTC[index]->Fill(invMass);
607
608         }
609         
610         if(fReadMC){
611           if(labDp>=0) {
612             index=GetSignalHistoIndex(iPtBin);
613             fMassHist[index]->Fill(invMass);
614
615             if(passTightCuts){
616               fMassHistTC[index]->Fill(invMass);
617
618             }
619             
620           }else{
621             index=GetBackgroundHistoIndex(iPtBin);
622             fMassHist[index]->Fill(invMass);
623
624             if(passTightCuts){
625               fMassHistTC[index]->Fill(invMass);
626
627             }   
628           }
629         }
630       }
631       /*
632       //start OS analysis
633       if(labDp<0)fHistOSbkg->Fill(d->InvMassDplus());
634       fHistOS->Fill(d->InvMassDplus());
635       */
636       nOS++;
637     }
638     if(unsetvtx) d->UnsetOwnPrimaryVtx();
639   }
640  
641   //start LS analysis
642   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
643   
644   PostData(1,fOutput);    
645   return;
646 }
647
648
649
650 //________________________________________________________________________
651 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
652 {
653   // Terminate analysis
654   //
655   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
656
657   fOutput = dynamic_cast<TList*> (GetOutputData(1));
658   if (!fOutput) {     
659     printf("ERROR: fOutput not available\n");
660     return;
661   }
662  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
663
664  TString hisname;
665  Int_t index=0;
666  Int_t indexLS=0;
667  for(Int_t i=0;i<fNPtBins;i++){
668     index=GetHistoIndex(i);
669     if(fDoLS)indexLS=GetLSHistoIndex(i);
670     hisname.Form("hMassPt%d",i);
671     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
672  
673     hisname.Form("hMassPt%dTC",i);
674     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
675     if(fDoLS){
676       hisname.Form("hLSPt%dTC",i);
677       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
678  
679       hisname.Form("hLSPt%dLC",i);
680       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
681     } 
682     
683     index=GetSignalHistoIndex(i);    
684     if(fDoLS)indexLS++;
685     hisname.Form("hSigPt%d",i);
686     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
687  
688     hisname.Form("hSigPt%dTC",i);
689     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
690     if(fDoLS){
691       hisname.Form("hLSPt%dLCnw",i);
692       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
693  
694       hisname.Form("hLSPt%dTCnw",i);
695       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
696     }
697
698     index=GetBackgroundHistoIndex(i); 
699     if(fDoLS)indexLS++;
700     hisname.Form("hBkgPt%d",i);
701     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
702  
703     hisname.Form("hBkgPt%dTC",i);
704     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
705     if(fDoLS){
706       hisname.Form("hLSPt%dLCntrip",i);
707       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
708  
709       hisname.Form("hLSPt%dTCntrip",i);
710       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
711     
712
713       indexLS++;
714       hisname.Form("hLSPt%dLCntripsinglecut",i);
715       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
716       
717       hisname.Form("hLSPt%dTCntripsinglecut",i);
718       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
719       
720       
721       indexLS++;
722       hisname.Form("hLSPt%dLCspc",i);
723       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
724       
725       hisname.Form("hLSPt%dTCspc",i);
726       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
727     }
728  
729   }
730
731   if(fFillNtuple){
732     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
733   }
734
735   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
736   c1->cd();
737   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
738   hMassPt->SetLineColor(kBlue);
739   hMassPt->SetXTitle("M[GeV/c^{2}]"); 
740   hMassPt->Draw();
741  
742  return;
743 }