6969eb8beeb24419705d12b096869d30e3d9bb23
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStarSpectra.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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  * appeuear 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 //                  Base class for DStar Analysis
21 //
22 //
23 //  The D* spectra study is done in pt bins:
24 //  [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
25 //  [8,10],[10,12], [12,16], [16,20] and [20,24]
26 //
27 //  Cuts arew centralized in AliRDHFCutsDStartoKpipi
28 //  Side Band and like sign background are implemented in the macro
29 //
30 //-----------------------------------------------------------------------
31 //
32 //                         Author A.Grelli 
33 //              ERC-QGP Utrecht University - a.grelli@uu.nl,
34 //                         Author Y.Wang
35 //        University of Heidelberg - yifei@physi.uni-heidelberg.de
36 //                         Author C.Ivan 
37 //             ERC-QGP Utrecht University - c.ivan@uu.nl,
38 //
39 //-----------------------------------------------------------------------
40
41 #include <TSystem.h>
42 #include <TParticle.h>
43 #include <TH1I.h>
44 #include "TROOT.h"
45 #include <TDatabasePDG.h>
46 #include <AliAnalysisDataSlot.h>
47 #include <AliAnalysisDataContainer.h>
48 #include "AliRDHFCutsDStartoKpipi.h"
49 #include "AliStack.h"
50 #include "AliMCEvent.h"
51 #include "AliAnalysisManager.h"
52 #include "AliAODMCHeader.h"
53 #include "AliAODHandler.h"
54 #include "AliLog.h"
55 #include "AliAODVertex.h"
56 #include "AliAODRecoDecay.h"
57 #include "AliAODRecoDecayHF.h"
58 #include "AliAODRecoCascadeHF.h"
59 #include "AliAODRecoDecayHF2Prong.h"
60 #include "AliAnalysisVertexingHF.h"
61 #include "AliESDtrack.h"
62 #include "AliAODMCParticle.h"
63 #include "AliAnalysisTaskSE.h"
64 #include "AliAnalysisTaskSEDStarSpectra.h"
65 #include "AliNormalizationCounter.h"
66
67 ClassImp(AliAnalysisTaskSEDStarSpectra)
68
69 //__________________________________________________________________________
70 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():  
71   AliAnalysisTaskSE(),
72   fEvents(0),
73   fAnalysis(0),
74   fD0Window(0),
75   fPeakWindow(0),
76   fUseMCInfo(kFALSE),
77   fDoSearch(kFALSE),
78   fOutput(0),
79   fOutputAll(0),
80   fOutputPID(0),
81   fNSigma(3),
82   fCuts(0),
83   fCEvents(0),     
84   fTrueDiff2(0),
85   fDeltaMassD1(0),
86   fCounter(0)
87 {
88   //
89   // Default ctor
90   //
91 }
92 //___________________________________________________________________________
93 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
94   AliAnalysisTaskSE(name),
95   fEvents(0),
96   fAnalysis(0),
97   fD0Window(0),
98   fPeakWindow(0),
99   fUseMCInfo(kFALSE),
100   fDoSearch(kFALSE),
101   fOutput(0),
102   fOutputAll(0),
103   fOutputPID(0),
104   fNSigma(3),
105   fCuts(0),
106   fCEvents(0),     
107   fTrueDiff2(0),
108   fDeltaMassD1(0),
109   fCounter(0)
110 {
111   //
112   // Constructor. Initialization of Inputs and Outputs
113   //
114   Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
115
116   fCuts=cuts;
117
118   DefineOutput(1,TList::Class());  //conters
119   DefineOutput(2,TList::Class());  //All Entries output
120   DefineOutput(3,TList::Class());  //3sigma PID output
121   DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());   //My private output
122   DefineOutput(5,AliNormalizationCounter::Class());   // normalization
123 }
124
125 //___________________________________________________________________________
126 AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
127   //
128   // destructor
129   //
130   Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
131   
132   if (fOutput) {
133     delete fOutput;
134     fOutput = 0;
135   }
136   if (fOutputAll) {
137     delete fOutputAll;
138     fOutputAll = 0;
139   }
140   if (fOutputPID) {
141     delete fOutputPID;
142     fOutputPID = 0;
143   }
144   if (fCuts) {
145     delete fCuts;
146     fCuts = 0;
147   }
148   if(fCEvents){
149     delete fCEvents;
150     fCEvents =0;
151   }
152   if(fDeltaMassD1){
153     delete fDeltaMassD1;
154     fDeltaMassD1 =0;
155   }
156 }
157 //_________________________________________________
158 void AliAnalysisTaskSEDStarSpectra::Init(){
159   //
160   // Initialization
161   //
162
163   if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
164    AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
165   // Post the data
166   PostData(4,copyfCuts);
167
168   return;
169 }
170
171 //_________________________________________________
172 void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
173 {
174   // user exec
175   if (!fInputEvent) {
176     Error("UserExec","NO EVENT FOUND!");
177     return;
178   }
179
180   fEvents++;
181
182   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
183   TClonesArray *arrayDStartoD0pi=0;
184
185   fCEvents->Fill(1);
186
187   if(!aodEvent && AODEvent() && IsStandardAOD()) {
188     // In case there is an AOD handler writing a standard AOD, use the AOD 
189     // event in memory rather than the input (ESD) event.    
190     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
191     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
192     // have to taken from the AOD event hold by the AliAODExtension
193     AliAODHandler* aodHandler = (AliAODHandler*) 
194       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
195     if(aodHandler->GetExtensions()) {
196       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
197       AliAODEvent *aodFromExt = ext->GetAOD();
198       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
199     }
200   } else {
201     arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
202   }
203
204   // fix for temporary bug in ESDfilter 
205   // the AODs with null vertex pointer didn't pass the PhysSel
206   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
207  
208   fCEvents->Fill(2);
209
210   fCounter->StoreEvent(aodEvent,fUseMCInfo);
211
212   if(!fCuts->IsEventSelected(aodEvent)) {
213     if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
214       fCEvents->Fill(6);
215     return;
216   }
217
218   // Load the event
219   AliInfo(Form("Event %d",fEvents));
220   if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
221
222   // counters for efficiencies
223   Int_t icountReco = 0;
224   
225   //D* and D0 prongs needed to MatchToMC method
226   Int_t pdgDgDStartoD0pi[2]={421,211};
227   Int_t pdgDgD0toKpi[2]={321,211};
228
229   // AOD primary vertex
230   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
231   if(!vtx1) return;
232   if(vtx1->GetNContributors()<1) return;
233
234   if (!arrayDStartoD0pi){
235     AliInfo("Could not find array of HF vertices, skipping the event");
236     return;
237   }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast())); 
238
239   Int_t nSelectedAna =0;
240   Int_t nSelectedProd =0;
241
242   // loop over the tracks to search for candidates soft pion 
243   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
244
245     // D* candidates and D0 from D*
246     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
247     if(!dstarD0pi->GetSecondaryVtx()) continue;
248     AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
249     if (!theD0particle) continue;
250     
251     Int_t isDStar = 0;   
252     TClonesArray *mcArray = 0; // fix coverity
253
254     // mc analysis 
255     if(fUseMCInfo){
256     //MC array need for maching
257       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
258       if (!mcArray) {
259         AliError("Could not find Monte-Carlo in AOD");
260         return;
261       }
262       // find associated MC particle for D* ->D0toKpi
263       Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
264       if(mcLabel>=0) isDStar = 1;
265     }
266     
267     Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
268     
269     // set the D0 search window bin by bin - useful to calculate side band bkg
270     if (ptbin==0){
271       if(fAnalysis==1){
272         fD0Window=0.035;
273         fPeakWindow=0.03;
274       }else{
275         fD0Window=0.020;
276         fPeakWindow=0.0018;
277       }
278     }
279     if (ptbin==1){
280       if(fAnalysis==1){
281         fD0Window=0.035;
282         fPeakWindow=0.03;
283       }else{
284         fD0Window=0.020;
285         fPeakWindow=0.0018;
286       }
287     }
288     if (ptbin==2){
289       if(fAnalysis==1){
290         fD0Window=0.035;
291         fPeakWindow=0.03;
292       }else{
293         fD0Window=0.020;
294         fPeakWindow=0.0018;
295       }
296     }
297     if (ptbin==3){
298       if(fAnalysis==1){
299         fD0Window=0.035;
300         fPeakWindow=0.03;
301       }else{
302         fD0Window=0.022;
303         fPeakWindow=0.0016;
304       }
305     }
306     if (ptbin==4){ 
307       if(fAnalysis==1){
308         fD0Window=0.035;
309         fPeakWindow=0.03;
310       }else{
311         fD0Window=0.026;
312         fPeakWindow=0.0014;
313       }
314     }
315     if (ptbin==5){
316       if(fAnalysis==1){
317         fD0Window=0.045;
318         fPeakWindow=0.03;
319       }else{
320         fD0Window=0.026;
321         fPeakWindow=0.0014;
322       }
323     } 
324    if (ptbin==6){
325       if(fAnalysis==1){
326         fD0Window=0.045;
327         fPeakWindow=0.03;
328       }else{
329         fD0Window=0.026;
330         fPeakWindow=0.006;
331       }
332     } 
333     if (ptbin==7){
334       if(fAnalysis==1){
335         fD0Window=0.055;
336         fPeakWindow=0.03;
337       }else{
338         fD0Window=0.026;
339         fPeakWindow=0.006;
340       }
341     }
342     if (ptbin>7){
343       if(fAnalysis==1){
344         fD0Window=0.074;
345         fPeakWindow=0.03;
346       }else{
347         fD0Window=0.026;
348         fPeakWindow=0.006;
349       }
350     }
351     fCEvents->Fill(9);
352     
353     nSelectedProd++;
354     nSelectedAna++;
355     
356     Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
357     
358     // here for lead-lead
359     Double_t invmassD0   = dstarD0pi->InvMassD0();  
360     if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) continue; 
361     
362     Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
363     Double_t invmassDelta = dstarD0pi->DeltaInvMass();
364     
365     if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
366     
367     Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
368     if(!isTkSelected) continue;
369     
370     if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;    
371     
372     Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
373     if (!isSelected) continue;
374
375     // fill PID
376     FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputPID);
377     SideBandBackground(dstarD0pi,fCuts,fOutputPID);
378     WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
379
380     // fill no pid
381     fCuts->SetUsePID(kFALSE);
382     FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputAll);
383     SideBandBackground(dstarD0pi,fCuts,fOutputAll);
384     WrongSignForDStar(dstarD0pi,fCuts,fOutputAll);
385     fCuts->SetUsePID(kTRUE);
386
387     // rare D search ------ 
388     if(fDoSearch){
389       TLorentzVector LorentzTrack1(0,0,0,0); // lorentz 4 vector
390       TLorentzVector LorentzTrack2(0,0,0,0); // lorentz 4 vector
391       
392       for (Int_t i=0; i<aodEvent->GetNTracks(); i++){ 
393         
394         AliAODTrack* aodTrack = aodEvent->GetTrack(i);
395         
396         if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
397         if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
398         if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
399         
400         //build the D1 mass
401         Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
402
403         LorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
404         LorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
405         
406         //D1 mass
407         Double_t d1mass = ((LorentzTrack1+LorentzTrack2).M());
408         //mass difference - at 0.4117 and 0.4566
409         fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
410       }
411     }
412
413     if(isDStar == 1) { 
414       fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
415     }
416     
417   }
418   
419   fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);  
420   fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE); 
421
422   AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
423   
424   PostData(1,fOutput);
425   PostData(2,fOutputAll);
426   PostData(3,fOutputPID);
427   PostData(5,fCounter);
428
429 }
430 //________________________________________ terminate ___________________________
431 void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
432 {    
433   // The Terminate() function is the last function to be called during
434   // a query. It always runs on the client, it can be used to present
435   // the results graphically or save the results to file.
436   
437   //Info("Terminate","");
438   AliAnalysisTaskSE::Terminate();
439   
440   fOutput = dynamic_cast<TList*> (GetOutputData(1));
441   if (!fOutput) {     
442     printf("ERROR: fOutput not available\n");
443     return;
444   }
445   
446   fCEvents        = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
447   fDeltaMassD1     = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
448   fTrueDiff2      = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
449
450   fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
451   if (!fOutputAll) {
452     printf("ERROR: fOutputAll not available\n");
453     return;
454   }
455   fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
456   if (!fOutputPID) {
457     printf("ERROR: fOutputPID not available\n");
458     return;
459   }
460  
461
462   return;
463 }
464 //___________________________________________________________________________
465 void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() { 
466  // output
467   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
468   
469   //slot #1  
470   //OpenFile(1);
471   fOutput = new TList();
472   fOutput->SetOwner();
473   fOutput->SetName("chist0");
474
475   fOutputAll = new TList();
476   fOutputAll->SetOwner();
477   fOutputAll->SetName("listAll");
478
479   fOutputPID = new TList();
480   fOutputPID->SetOwner();
481   fOutputPID->SetName("listPID");
482     
483   // define histograms
484   DefineHistograms();
485
486  //Counter for Normalization
487   TString normName="NormalizationCounter";
488   AliAnalysisDataContainer *cont = GetOutputSlot(4)->GetContainer();
489   if(cont)normName=(TString)cont->GetName();
490   fCounter = new AliNormalizationCounter(normName.Data());
491   fCounter->SetRejectPileUp(fCuts->GetOptPileUp());
492
493   PostData(1,fOutput);
494   PostData(2,fOutputAll);
495   PostData(3,fOutputPID);
496
497   return;
498 }
499 //___________________________________ hiostograms _______________________________________
500 void  AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
501
502   fCEvents = new TH1F("fCEvents","conter",11,0,11);
503   fCEvents->SetStats(kTRUE);
504   fCEvents->GetXaxis()->SetTitle("1");
505   fCEvents->GetYaxis()->SetTitle("counts");
506   fOutput->Add(fCEvents);
507
508   fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
509   fOutput->Add(fTrueDiff2);
510
511   fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
512   fOutput->Add(fDeltaMassD1);
513
514   const Int_t nhist=14;
515   TString nameMass=" ", nameSgn=" ", nameBkg=" ";
516
517   for(Int_t i=-2;i<nhist;i++){
518     nameMass="histDeltaMass_";
519     nameMass+=i+1;
520     nameSgn="histDeltaSgn_";
521     nameSgn+=i+1;
522     nameBkg="histDeltaBkg_";
523     nameBkg+=i+1; 
524     
525     if (i==-2) {
526       nameMass="histDeltaMass";
527       nameSgn="histDeltaSgn";
528       nameBkg="histDeltaBkg";
529     }
530     
531     TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
532     TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
533     TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
534     
535     nameMass="histD0Mass_";
536     nameMass+=i+1;
537     nameSgn="histD0Sgn_";
538     nameSgn+=i+1;
539     nameBkg="histD0Bkg_";
540     nameBkg+=i+1;
541     
542     if (i==-2) {
543       nameMass="histD0Mass";
544       nameSgn="histD0Sgn";
545       nameBkg="histD0Bkg";
546     }
547
548     TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
549     TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
550     TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
551
552     nameMass="histDstarMass_";
553     nameMass+=i+1;
554     nameSgn="histDstarSgn_";
555     nameSgn+=i+1;
556     nameBkg="histDstarBkg_";
557     nameBkg+=i+1;
558
559     if (i==-2) {
560       nameMass="histDstarMass";
561       nameSgn="histDstarSgn";
562       nameBkg="histDstarBkg";
563     }
564
565     TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
566     TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
567     TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
568
569     nameMass="histSideBandMass_";
570     nameMass+=i+1;
571     if (i==-2) { 
572       nameMass="histSideBandMass";
573     }
574     
575     TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
576
577     nameMass="histWrongSignMass_";
578     nameMass+=i+1;
579     if (i==-2) { 
580       nameMass="histWrongSignMass";
581     }
582     
583     TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
584
585
586     spectrumMass->Sumw2();
587     spectrumSgn->Sumw2();
588     spectrumBkg->Sumw2();
589     
590     spectrumMass->SetLineColor(6);
591     spectrumSgn->SetLineColor(2);
592     spectrumBkg->SetLineColor(4);
593     
594     spectrumMass->SetMarkerStyle(20);
595     spectrumSgn->SetMarkerStyle(20);
596     spectrumBkg->SetMarkerStyle(20);
597     spectrumMass->SetMarkerSize(0.6);
598     spectrumSgn->SetMarkerSize(0.6);
599     spectrumBkg->SetMarkerSize(0.6);
600     spectrumMass->SetMarkerColor(6);
601     spectrumSgn->SetMarkerColor(2);
602     spectrumBkg->SetMarkerColor(4);
603
604     spectrumD0Mass->Sumw2();
605     spectrumD0Sgn->Sumw2();
606     spectrumD0Bkg->Sumw2();
607
608     spectrumD0Mass->SetLineColor(6);
609     spectrumD0Sgn->SetLineColor(2);
610     spectrumD0Bkg->SetLineColor(4);
611
612     spectrumD0Mass->SetMarkerStyle(20);
613     spectrumD0Sgn->SetMarkerStyle(20);
614     spectrumD0Bkg->SetMarkerStyle(20);
615     spectrumD0Mass->SetMarkerSize(0.6);
616     spectrumD0Sgn->SetMarkerSize(0.6);
617     spectrumD0Bkg->SetMarkerSize(0.6);
618     spectrumD0Mass->SetMarkerColor(6);
619     spectrumD0Sgn->SetMarkerColor(2);
620     spectrumD0Bkg->SetMarkerColor(4);
621
622     spectrumDstarMass->Sumw2();
623     spectrumDstarSgn->Sumw2();
624     spectrumDstarBkg->Sumw2();
625
626     spectrumDstarMass->SetLineColor(6);
627     spectrumDstarSgn->SetLineColor(2);
628     spectrumDstarBkg->SetLineColor(4);
629
630     spectrumDstarMass->SetMarkerStyle(20);
631     spectrumDstarSgn->SetMarkerStyle(20);
632     spectrumDstarBkg->SetMarkerStyle(20);
633     spectrumDstarMass->SetMarkerSize(0.6);
634     spectrumDstarSgn->SetMarkerSize(0.6);
635     spectrumDstarBkg->SetMarkerSize(0.6);
636     spectrumDstarMass->SetMarkerColor(6);
637     spectrumDstarSgn->SetMarkerColor(2);
638     spectrumDstarBkg->SetMarkerColor(4);
639
640     spectrumSideBandMass->Sumw2();
641     spectrumSideBandMass->SetLineColor(4);
642     spectrumSideBandMass->SetMarkerStyle(20);
643     spectrumSideBandMass->SetMarkerSize(0.6);
644     spectrumSideBandMass->SetMarkerColor(4);
645
646     spectrumWrongSignMass->Sumw2();
647     spectrumWrongSignMass->SetLineColor(4);
648     spectrumWrongSignMass->SetMarkerStyle(20);
649     spectrumWrongSignMass->SetMarkerSize(0.6);
650     spectrumWrongSignMass->SetMarkerColor(4);
651
652     TH1F* allMass = (TH1F*)spectrumMass->Clone();
653     TH1F* allSgn  = (TH1F*)spectrumSgn->Clone();
654     TH1F* allBkg  = (TH1F*)spectrumBkg->Clone();
655
656     TH1F* pidMass = (TH1F*)spectrumMass->Clone();
657     TH1F* pidSgn  = (TH1F*)spectrumSgn->Clone();
658     TH1F* pidBkg  = (TH1F*)spectrumBkg->Clone();
659
660     fOutputAll->Add(allMass);
661     fOutputAll->Add(allSgn);
662     fOutputAll->Add(allBkg);
663
664     fOutputPID->Add(pidMass);
665     fOutputPID->Add(pidSgn);
666     fOutputPID->Add(pidBkg);
667
668     TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
669     TH1F* allD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
670     TH1F* allD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
671
672     TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
673     TH1F* pidD0Sgn  = (TH1F*)spectrumD0Sgn->Clone();
674     TH1F* pidD0Bkg  = (TH1F*)spectrumD0Bkg->Clone();
675
676     fOutputAll->Add(allD0Mass);
677     fOutputAll->Add(allD0Sgn);
678     fOutputAll->Add(allD0Bkg);
679
680     fOutputPID->Add(pidD0Mass);
681     fOutputPID->Add(pidD0Sgn);
682     fOutputPID->Add(pidD0Bkg);
683   
684     TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
685     TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
686     TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
687     
688     TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
689     TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
690     TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
691     
692     fOutputAll->Add(allDstarMass);
693     fOutputAll->Add(allDstarSgn);
694     fOutputAll->Add(allDstarBkg);
695
696     fOutputPID->Add(pidDstarMass);
697     fOutputPID->Add(pidDstarSgn);
698     fOutputPID->Add(pidDstarBkg);
699
700     TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
701     TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
702
703     fOutputAll->Add(allSideBandMass);
704     fOutputPID->Add(pidSideBandMass);
705    
706     TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
707     TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
708
709     fOutputAll->Add(allWrongSignMass);
710     fOutputPID->Add(pidWrongSignMass);
711    
712   }
713   
714   // pt spectra
715   nameMass="ptMass";
716   nameSgn="ptSgn";
717   nameBkg="ptBkg";
718   
719   TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
720   TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
721   TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
722   
723   ptspectrumMass->Sumw2();
724   ptspectrumSgn->Sumw2();
725   ptspectrumBkg->Sumw2();
726   
727   ptspectrumMass->SetLineColor(6);
728   ptspectrumSgn->SetLineColor(2);
729   ptspectrumBkg->SetLineColor(4);
730   
731   ptspectrumMass->SetMarkerStyle(20);
732   ptspectrumSgn->SetMarkerStyle(20);
733   ptspectrumBkg->SetMarkerStyle(20);
734   ptspectrumMass->SetMarkerSize(0.6);
735   ptspectrumSgn->SetMarkerSize(0.6);
736   ptspectrumBkg->SetMarkerSize(0.6);
737   ptspectrumMass->SetMarkerColor(6);
738   ptspectrumSgn->SetMarkerColor(2);
739   ptspectrumBkg->SetMarkerColor(4);
740   
741   TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
742   TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
743   TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
744   
745   TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
746   TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
747   TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
748     
749   fOutputAll->Add(ptallMass);
750   fOutputAll->Add(ptallSgn);
751   fOutputAll->Add(ptallBkg);
752   
753   fOutputPID->Add(ptpidMass);
754   fOutputPID->Add(ptpidSgn);
755   fOutputPID->Add(ptpidBkg);
756   
757   // eta spectra
758   nameMass="etaMass";
759   nameSgn="etaSgn";
760   nameBkg="etaBkg";
761   
762   TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
763   TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
764   TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
765   
766   etaspectrumMass->Sumw2();
767   etaspectrumSgn->Sumw2();
768   etaspectrumBkg->Sumw2();
769   
770   etaspectrumMass->SetLineColor(6);
771   etaspectrumSgn->SetLineColor(2);
772   etaspectrumBkg->SetLineColor(4);
773   
774   etaspectrumMass->SetMarkerStyle(20);
775   etaspectrumSgn->SetMarkerStyle(20);
776   etaspectrumBkg->SetMarkerStyle(20);
777   etaspectrumMass->SetMarkerSize(0.6);
778   etaspectrumSgn->SetMarkerSize(0.6);
779   etaspectrumBkg->SetMarkerSize(0.6);
780   etaspectrumMass->SetMarkerColor(6);
781   etaspectrumSgn->SetMarkerColor(2);
782   etaspectrumBkg->SetMarkerColor(4);
783   
784   TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
785   TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
786   TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
787   
788   TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
789   TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
790   TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
791     
792   fOutputAll->Add(etaallMass);
793   fOutputAll->Add(etaallSgn);
794   fOutputAll->Add(etaallBkg);
795   
796   fOutputPID->Add(etapidMass);
797   fOutputPID->Add(etapidSgn);
798   fOutputPID->Add(etapidBkg);
799   
800   return;
801 }
802 //________________________________________________________________________
803 void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
804   //
805   // Fill histos for D* spectrum
806   //
807
808   Int_t ptbin=cuts->PtBin(part->Pt());
809   Double_t invmassD0   = part->InvMassD0();  
810   
811   Double_t pt = part->Pt();
812   Double_t eta = part->Eta();  
813   
814   Double_t invmassDelta = part->DeltaInvMass();
815   Double_t invmassDstar = part->InvMassDstarKpipi();
816   
817   TString fillthis="";
818   Bool_t massInRange=kFALSE;
819   
820   Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();    
821   Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
822   
823   // delta M(Kpipi)-M(Kpi)
824   if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;  
825   
826   if(fUseMCInfo) {
827     if(isDStar==1) {
828       fillthis="histD0Sgn_";
829       fillthis+=ptbin;
830       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
831       fillthis="histD0Sgn";
832       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
833       fillthis="histDstarSgn_";
834       fillthis+=ptbin;
835       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
836       fillthis="histDstarSgn";
837       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
838       fillthis="histDeltaSgn_";
839       fillthis+=ptbin;
840       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
841       fillthis="histDeltaSgn";
842       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
843     if (massInRange) {
844         fillthis="ptSgn";
845         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
846         fillthis="etaSgn";
847         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
848       }
849     }
850     else {//background
851       fillthis="histD0Bkg_";
852       fillthis+=ptbin;
853       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
854       fillthis="histD0Bkg";
855       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
856       fillthis="histDstarBkg_";
857       fillthis+=ptbin;
858       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
859       fillthis="histDstarBkg";
860       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
861       fillthis="histDeltaBkg_";
862       fillthis+=ptbin;
863       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
864       fillthis="histDeltaBkg";
865       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
866      if (massInRange) {
867         fillthis="ptBkg";
868         ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
869         fillthis="etaBkg";
870         ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
871       }
872     }
873   }
874   //no MC info, just cut selection
875   fillthis="histD0Mass_";
876   fillthis+=ptbin;
877   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
878   fillthis="histD0Mass";
879   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
880   fillthis="histDstarMass_";
881   fillthis+=ptbin;
882   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
883   fillthis="histDstarMass";
884   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
885   fillthis="histDeltaMass_";
886   fillthis+=ptbin;
887   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
888   fillthis="histDeltaMass";
889   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
890   
891   if (massInRange) {
892     fillthis="ptMass";
893     ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
894     fillthis="etaMass";
895     ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
896   }
897  
898   return;
899 }
900 //______________________________ side band background for D*___________________________________
901 void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){
902
903   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
904   // (expected detector resolution) on the left and right frm the D0 mass. Each band
905   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
906
907   Int_t ptbin=cuts->PtBin(part->Pt());
908   
909   Bool_t massInRange=kFALSE;
910   
911   // select the side bands intervall
912   Double_t invmassD0    = part->InvMassD0();
913   if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
914     
915     // for pt and eta
916     Double_t invmassDelta = part->DeltaInvMass();
917     if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
918     
919     TString fillthis="";
920     fillthis="histSideBandMass_";
921     fillthis+=ptbin;
922     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
923     fillthis="histSideBandMass";
924     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
925     
926   }
927 }
928 //________________________________________________________________________________________________________________
929 void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part,  AliRDHFCutsDStartoKpipi *cuts, TList *listout){
930   //
931   // assign the wrong charge to the soft pion to create background
932   //
933   Int_t ptbin=cuts->PtBin(part->Pt());
934   
935   AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
936
937   Int_t okD0WrongSign,okD0barWrongSign;
938   Double_t wrongMassD0=0.;
939   
940   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
941    if (!isSelected){
942     return;
943   }
944
945   okD0WrongSign =  1;
946   okD0barWrongSign = 1;
947   
948   //if is D*+ than assume D0bar
949   if(part->Charge()>0 && (isSelected ==1)) { 
950     okD0WrongSign = 0;
951   }
952   if(part->Charge()<0 && (isSelected ==2)){
953     okD0barWrongSign = 0;
954   }
955   
956   // assign the wrong mass in case the cuts return both D0 and D0bar
957   if(part->Charge()>0 && (isSelected ==3)) { 
958     okD0WrongSign = 0;
959   } else if(part->Charge()<0 && (isSelected ==3)){
960     okD0barWrongSign = 0;
961   }
962   
963   //wrong D0 inv mass
964   if(okD0WrongSign!=0){
965     wrongMassD0 = theD0particle->InvMassD0();
966   }else if(okD0WrongSign==0){
967     wrongMassD0 = theD0particle->InvMassD0bar();
968   }
969   
970   if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
971     
972     // wrong D* inv mass   
973     Double_t e[3];
974     if (part->Charge()>0){
975       e[0]=theD0particle->EProng(0,321);
976       e[1]=theD0particle->EProng(1,211);
977     }else{
978       e[0]=theD0particle->EProng(0,211);
979       e[1]=theD0particle->EProng(1,321);
980     }
981     e[2]=part->EProng(0,211);
982     
983     Double_t esum = e[0]+e[1]+e[2];
984     Double_t pds = part->P();
985
986     Double_t   wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
987
988     TString fillthis="";
989     fillthis="histWrongSignMass_";
990     fillthis+=ptbin;
991     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
992     fillthis="histWrongSignMass";
993     ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
994     
995   }
996 }