Incrementing class version.
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskDStarCorrelations.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  * 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 //             Base class for DStar - Hadron Correlations Analysis
18 //
19 //-----------------------------------------------------------------------
20 //          
21 //
22 //                                                 Author S.Bjelogrlic
23 //                         Utrecht University 
24 //                      sandro.bjelogrlic@cern.ch
25 //
26 //-----------------------------------------------------------------------
27
28 /* $Id: AliAnalysisTaskDStarCorrelations.cxx 65139 2013-11-25 14:47:45Z fprino $ */
29
30 //#include <TDatabasePDG.h>
31 #include <TParticle.h>
32 #include <TVector3.h>
33 #include <TChain.h>
34 #include "TROOT.h"
35
36 #include "AliAnalysisTaskDStarCorrelations.h"
37 #include "AliRDHFCutsDStartoKpipi.h"
38 #include "AliHFAssociatedTrackCuts.h"
39 #include "AliAODRecoDecay.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAODRecoDecayHF2Prong.h"
42 #include "AliAODPidHF.h"
43 #include "AliVParticle.h"
44 #include "AliAnalysisManager.h"
45 #include "AliAODInputHandler.h"
46 #include "AliAODHandler.h"
47 #include "AliESDtrack.h"
48 #include "AliAODMCParticle.h"
49 #include "AliNormalizationCounter.h"
50 #include "AliReducedParticle.h"
51 #include "AliHFCorrelator.h"
52 #include "AliAODMCHeader.h"
53 #include "AliEventPoolManager.h"
54 #include "AliVertexingHFUtils.h"
55
56 using std::cout;
57 using std::endl;
58
59
60 ClassImp(AliAnalysisTaskDStarCorrelations)
61
62
63 //__________________________________________________________________________
64 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
65 AliAnalysisTaskSE(),
66 fhandler(0x0),
67 fmcArray(0x0),
68 fCounter(0x0),
69 fCorrelator(0x0),
70 fselect(0),
71 fmontecarlo(kFALSE),
72 fmixing(kFALSE),
73 fFullmode(kFALSE),
74 fSystem(pp),
75 fEfficiencyVariable(kNone),
76 fBkgMethod(kDZeroSB),
77 fReco(kTRUE),
78 fUseEfficiencyCorrection(kFALSE),
79 fUseDmesonEfficiencyCorrection(kFALSE),
80 fUseCentrality(kFALSE),
81 fUseHadronicChannelAtKineLevel(kFALSE),
82 fPhiBins(32),
83 fEvents(0),
84 fDebugLevel(0),
85 fDisplacement(0),
86 fDim(0),
87 fNofPtBins(0),
88 fMaxEtaDStar(0.9),
89 fDMesonSigmas(0),
90
91 fD0Window(0),
92
93 fOutput(0x0),
94 fDmesonOutput(0x0),
95 fTracksOutput(0x0),
96 fEMOutput(0x0),
97 fCorrelationOutput(0x0),
98 fOutputMC(0x0),
99
100 fCuts(0),
101 fAssocCuts(0),
102 fUtils(0),
103 fTracklets(0),
104 fDeffMapvsPt(0),
105 fDeffMapvsPtvsMult(0),
106 fDeffMapvsPtvsEta(0)
107 {
108     SetDim();
109     // default constructor
110 }
111
112 //__________________________________________________________________________
113 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :
114 AliAnalysisTaskSE(name),
115
116 fhandler(0x0),
117 fmcArray(0x0),
118 fCounter(0x0),
119 fCorrelator(0x0),
120 fselect(0),
121 fmontecarlo(kFALSE),
122 fmixing(kFALSE),
123 fFullmode(mode),
124 fSystem(syst),
125 fEfficiencyVariable(kNone),
126 fBkgMethod(kDZeroSB),
127 fReco(kTRUE),
128 fUseEfficiencyCorrection(kFALSE),
129 fUseDmesonEfficiencyCorrection(kFALSE),
130 fUseCentrality(kFALSE),
131 fUseHadronicChannelAtKineLevel(kFALSE),
132 fPhiBins(32),
133 fEvents(0),
134 fDebugLevel(0),
135 fDisplacement(0),
136 fDim(0),
137 fNofPtBins(0),
138 fMaxEtaDStar(0.9),
139 fDMesonSigmas(0),
140 fD0Window(0),
141
142 fOutput(0x0),
143 fDmesonOutput(0x0),
144 fTracksOutput(0x0),
145 fEMOutput(0x0),
146 fCorrelationOutput(0x0),
147 fOutputMC(0x0),
148
149 fCuts(0),
150 fAssocCuts(AsscCuts),
151 fUtils(0),
152 fTracklets(0),
153 fDeffMapvsPt(0),
154 fDeffMapvsPtvsMult(0),
155 fDeffMapvsPtvsEta(0)
156 {
157      Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
158   
159     SetDim();
160     if(fSystem == pp)  fUseCentrality = kFALSE; else fUseCentrality = kTRUE;
161     
162     if(fSystem == AA) fBkgMethod = kDStarSB; else fBkgMethod = kDZeroSB;
163   
164     fCuts=cuts;
165     fNofPtBins= fCuts->GetNPtBins();
166     //cout << "Enlarging the DZero window " << endl;
167     EnlargeDZeroMassWindow();
168    // cout << "Done" << endl;
169     
170    
171   DefineInput(0, TChain::Class());
172   
173   DefineOutput(1,TList::Class()); // histos from data and MC
174   DefineOutput(2,TList::Class()); // histos from MC
175     DefineOutput(3,TList::Class()); // histos from data and MC
176     DefineOutput(4,TList::Class()); // histos from MC
177     DefineOutput(5,TList::Class()); // histos from data and MC
178       DefineOutput(6,TList::Class()); // histos from data and MC
179     DefineOutput(7,AliNormalizationCounter::Class());   // normalization
180
181   DefineOutput(8,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
182   DefineOutput(9,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
183 }
184
185 //__________________________________________________________________________
186
187 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
188   //
189         // destructor
190         //
191         
192         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  
193         
194         if(fhandler) {delete fhandler; fhandler = 0;}
195         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
196         if(fmcArray) {delete fmcArray; fmcArray = 0;}
197         if(fCounter) {delete fCounter; fCounter = 0;}
198         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
199         if(fOutput) {delete fOutput; fOutput = 0;}
200         if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
201         if(fCuts) {delete fCuts; fCuts = 0;}
202         if(fAssocCuts) {delete fAssocCuts; fAssocCuts=0;}
203     if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}
204     if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}
205     if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}
206
207 }
208
209 //___________________________________________________________
210 void AliAnalysisTaskDStarCorrelations::Init(){
211         //
212         // Initialization
213         //
214         if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
215         
216         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
217         
218         
219     
220         // Post the D* cuts
221         PostData(8,copyfCuts);
222         
223         // Post the hadron cuts
224         PostData(9,fAssocCuts);
225     
226         return;
227 }
228
229
230 //_________________________________________________
231 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
232         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
233         
234         //slot #1
235         //OpenFile(0);
236         fOutput = new TList();
237         fOutput->SetOwner();
238         
239         fOutputMC = new TList();
240         fOutputMC->SetOwner();
241     
242     fDmesonOutput = new TList();
243         fDmesonOutput->SetOwner();
244     
245     fTracksOutput = new TList();
246         fTracksOutput->SetOwner();
247     
248     fEMOutput = new TList();
249         fEMOutput->SetOwner();
250     
251     fCorrelationOutput = new TList();
252         fCorrelationOutput->SetOwner();
253         
254         // define histograms
255         DefineHistoForAnalysis();
256     DefineThNSparseForAnalysis();
257     
258
259     
260
261         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(6)->GetContainer()->GetName()));
262         fCounter->Init();
263         
264     Double_t Pi = TMath::Pi();
265     //AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t useCentrality, AliRDHFCuts *cutObject)
266         fCorrelator = new AliHFCorrelator("Correlator",fAssocCuts,fUseCentrality,fCuts); // fAssocCuts is the hadron cut object, fSystem to switch between pp or PbPb
267         fCorrelator->SetDeltaPhiInterval(  -0.5*Pi, 1.5*Pi); // set correct phi interval
268         //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval
269         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
270         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
271         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
272         fCorrelator->SetUseMC(fmontecarlo);
273         fCorrelator->SetUseReco(fReco);
274     //  fCorrelator->SetKinkRemoval(kTRUE);
275         Bool_t pooldef = fCorrelator->DefineEventPool();
276         
277         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
278     
279     fUtils = new AliAnalysisUtils();
280     
281     
282         
283         PostData(1,fOutput); // set the outputs
284         PostData(2,fDmesonOutput); // set the outputs
285     PostData(3,fTracksOutput); // set the outputs
286     PostData(4,fEMOutput); // set the outputs
287     PostData(5,fCorrelationOutput); // set the outputs
288     PostData(6,fOutputMC); // set the outputs
289         PostData(7,fCounter); // set the outputs
290 }
291
292 //________________________________________  user exec ____________
293 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
294     // cout << "Task debug check 1 " << endl;
295      // ********************************************** LOAD THE EVENT ****************************************************
296     
297     if (!fInputEvent) {
298         Error("UserExec","NO EVENT FOUND!");
299                 return;
300     }
301     
302     AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
303     if(!aodEvent){
304         AliError("AOD event not found!");
305         return;
306     }
307     
308     fEvents++; // event counter
309     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
310     
311     fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); // store event in AliNormalizationCounter
312     
313     // load MC array
314     fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
315     if(fmontecarlo && !fmcArray){
316         AliError("Array of MC particles not found");
317         return;
318     }
319     
320     // ********************************************** START EVENT SELECTION ****************************************************
321     
322     Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
323     
324     if(!isEvSel) {
325         
326         if(fCuts->IsEventRejectedDueToPileup()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);                        cout << "Reject PILEUP" << endl;}
327         if(fCuts->IsEventRejectedDueToCentrality()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);                    cout << "Reject CENTRALITY" << endl;}
328         if(fCuts->IsEventRejectedDueToNotRecoVertex()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);                 cout << "Reject NOT RECO VTX" << endl;}
329         if(fCuts->IsEventRejectedDueToVertexContributors()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);            cout << "Reject VTX CONTRIB" << endl;}
330         if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);  cout << "Reject VTX no fid reg " << endl;}
331         if(fCuts->IsEventRejectedDueToTrigger()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);                       cout << "Reject TRIGGER" << endl;}
332         if(fCuts->IsEventRejectedDuePhysicsSelection()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);                cout << "Reject PHYS SEL" << endl;}
333         
334         return;
335     }
336     // added event selection for pA
337     
338     if(fSystem == pA){
339         
340         if(fUtils->IsFirstEventInChunk(aodEvent)) {
341             AliInfo("Rejecting the event - first in the chunk");
342             ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);
343             return;
344         }
345         if(!fUtils->IsVertexSelected2013pA(aodEvent)) {
346             ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);
347             AliInfo("Rejecting the event - bad vertex");
348             return;
349         }
350     }
351
352     fCorrelator->SetAODEvent(aodEvent); // set the event in the correlator class
353     
354     Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event within the pool definition
355         if(!correlatorON) {
356         AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
357         return; // if not the case, skips the event
358         }
359     
360     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); // count events that are passing the event selection
361     
362     
363     
364     // cout << "Task debug check 2 " << endl;
365      // ********************************************** CENTRALITY/MULTIPLICITY  ****************************************************
366
367
368     Double_t MultipOrCent = fCorrelator->GetCentrality(); // returns centrality or multiplicity
369     
370     
371     
372     // ********************************************** MC SETUP  ****************************************************
373     
374     if(fmontecarlo) {
375         
376         fCorrelator->SetMCArray(fmcArray); // export MC array into correlatior class
377         
378         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
379         if (fmontecarlo && !mcHeader) {
380             AliError("Could not find MC Header in AOD");
381             return;
382         }
383         
384       // select MC event type (PP, GS etc) - those are defined in the associated tracks cut object
385         Bool_t isMCeventgood = kFALSE;
386         
387         
388         Int_t eventType = mcHeader->GetEventType();
389         Int_t NMCevents = fAssocCuts->GetNofMCEventType();
390         
391         for(Int_t k=0; k<NMCevents; k++){
392             Int_t * MCEventType = fAssocCuts->GetMCEventType();
393             
394             if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
395             ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
396         }
397         
398         if(NMCevents && !isMCeventgood){
399             if(fDebugLevel)     std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
400             return;
401         }
402         
403         
404     }// end if fmontecarlo
405     
406     
407      // ********************************************** EVENT PROPERTIES  ****************************************************
408     // checks on vertex and multiplicity of the event
409         AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
410         Double_t zVtxPosition = vtx->GetZ(); // zvertex
411         
412  
413   
414     if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
415     
416     if(!fmixing) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);
417     if(!fmixing) ((TH1D*)fOutput->FindObject("MultiplicityOnlyCheck"))->Fill(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
418     
419     
420    // Int_t poolbin = fAssocCuts->GetPoolBin(MultipOrCent, zVtxPosition); // get the event pool bin - commented for the moment - to be checked
421     
422     
423     
424      // ********************************************** D * selection  ****************************************************
425     
426     
427     TClonesArray *arrayDStartoD0pi=0;
428         if(!aodEvent && AODEvent() && IsStandardAOD()) {
429         // In case there is an AOD handler writing a standard AOD, use the AOD
430         // event in memory rather than the input (ESD) event.
431         aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
432         // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
433         // have to taken from the AOD event hold by the AliAODExtension
434         AliAODHandler* aodHandler = (AliAODHandler*)
435             ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
436         if(aodHandler->GetExtensions()) {
437             AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
438             AliAODEvent *aodFromExt = ext->GetAOD();
439             arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
440         }
441         } else {
442         arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
443         }
444     
445     
446     
447     // D* related variables
448     
449     Double_t ptDStar;   // D* pt
450         Double_t phiDStar;  // D* phi
451         Double_t etaDStar;  // D* eta
452         Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; // flags for selection
453         Double_t invMassDZero; // D0 inv mass
454         Double_t deltainvMDStar; // D* - D0 invarian mass
455     
456     Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
457         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
458         
459         
460         //MC tagging for DStar
461         //D* and D0 prongs needed to MatchToMC method
462         Int_t pdgDgDStartoD0pi[2]={421,211};
463         Int_t pdgDgD0toKpi[2]={321,211};
464         
465         //Bool_t isDStarCand = kFALSE;
466     Bool_t isDfromB = kFALSE;
467         Bool_t isEventMixingFilledPeak = kFALSE;
468         Bool_t isEventMixingFilledSB = kFALSE;
469     Bool_t EventHasDStarCandidate = kFALSE;
470     Bool_t EventHasDZeroSideBandCandidate = kFALSE;
471     Bool_t EventHasDStarSideBandCandidate = kFALSE;
472     Bool_t isTrackForPeakFilled = kFALSE;
473     Bool_t isTrackForSBFilled = kFALSE;
474         //loop on D* candidates
475         
476         Int_t looponDCands = 0;
477         if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); // load N of D* candidates from AOD array
478         if(!fReco) looponDCands = fmcArray->GetEntriesFast(); // load N of D* candidates from MC array
479         
480         Int_t nOfDStarCandidates = 0;  // D candidates counter
481         Int_t nOfSBCandidates = 0;     // SB candidates counter
482         
483         Double_t DmesonEfficiency = 1.; // Efficiency of D meson
484         Double_t DmesonWeight = 1.;    // Applyed weight
485     Double_t efficiencyvariable = -999; // Variable to be used (defined by the AddTask)
486     Double_t dmDStarWindow = 0.0019/3; // DStar window
487     
488     // cout << "Task debug check 3 " << endl;
489     // loop over D meson candidates
490     for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
491     
492         // initialize variables
493         isInPeak = kFALSE;
494         isInDStarSideBand = kFALSE;
495         isInDZeroSideBand = kFALSE;
496         
497         EventHasDStarCandidate = kFALSE;
498         EventHasDZeroSideBandCandidate = kFALSE;
499         EventHasDStarSideBandCandidate = kFALSE;
500      
501         
502         isDStarMCtag = kFALSE;
503         isDfromB = kFALSE;
504         ptDStar = -123.4;
505         phiDStar = -999;
506         etaDStar = -56.;
507         invMassDZero = - 999;
508         deltainvMDStar = -998;
509         AliAODRecoCascadeHF* dstarD0pi;
510         AliAODRecoDecayHF2Prong* theD0particle;
511         AliAODMCParticle* DStarMC=0x0;
512         Short_t daughtercharge = -2;
513         Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
514         Int_t trackiddaugh1 = -1;
515         Int_t trackidsoftPi = -1;
516         Int_t ptbin = -1;
517         
518         // ============================================ using reconstruction on Data or MC
519         if(fReco){
520             // cout << "Task debug check 4 " << endl;
521             // get the D objects
522             dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
523             if(!dstarD0pi->GetSecondaryVtx()) continue;
524             theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
525             if (!theD0particle) continue;
526             
527             
528             // apply topological cuts
529             
530             // track quality cuts
531             Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
532             // region of interest + topological cuts + PID
533             Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
534             
535             //apply topological cuts
536             if(!isTkSelected) continue;
537             if(!isSelected) continue;
538             if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
539             
540             // get D candidate variables
541             ptDStar = dstarD0pi->Pt();
542             phiDStar = dstarD0pi->Phi();
543             etaDStar = dstarD0pi->Eta();
544             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
545             if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;
546             if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
547             if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
548             if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
549            
550             
551              if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
552             
553             phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
554             ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
555             
556             cout << "DStar pt = " << ptDStar << endl;
557              cout << "pt bin = " << ptbin << endl;
558             if(ptbin<1) continue;
559             
560              Double_t mD0Window= fD0Window[ptbin]/3;
561              invMassDZero = dstarD0pi->InvMassD0();
562              deltainvMDStar = dstarD0pi->DeltaInvMass();
563             
564             
565             // get the D meson efficiency
566             DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
567             
568             // check this!
569             if(fUseDmesonEfficiencyCorrection){
570                 if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
571                 else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI
572                     if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value
573                     else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low
574                 }
575             }
576             else DmesonWeight = 1.;
577             
578          
579             // using montecarlo on reconstruction
580             Int_t mcLabelDStar = -999;
581             if(fmontecarlo){
582                 // find associated MC particle for D* ->D0toKpi
583                 mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);
584                 if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
585                 if(!isDStarMCtag) continue;
586                 AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
587                 //check if DStar from B
588                 Int_t labelMother = MCDStar->GetMother();
589                 AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
590                 if(!mother) continue;
591                 Int_t motherPDG =TMath::Abs(mother->PdgCode());
592                 if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
593             }
594             
595             
596             // fill mass histograms
597             
598             // cout << "Task debug check 5 " << endl;
599             // fill D0 invariant mass
600             if(!fmixing) {
601                 cout << " histo name = " << Form("histDZeroMass_%d",ptbin) << endl;
602                 ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
603                 // cout << "Task debug check 5.1 " << endl;
604                 if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
605             } // end if !mixing
606             
607         
608             
609             
610             
611             // good D0 canidates
612             if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
613                 // fill D* invariant mass
614                 if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
615                    // fill D* invariant mass if weighting
616                     if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
617                     isInPeak = kTRUE;
618                     // fill other good candidate distributions - pt, phi eta
619                     if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)    {
620                         ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
621                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
622                         if(!fmixing)    ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
623                         nOfDStarCandidates++;
624                         EventHasDStarCandidate=kTRUE;
625                     } // end if in good D* mass window
626                     
627                     // count D* sideband candidates
628                     if(fBkgMethod == kDStarSB ){
629                         if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
630                             ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
631                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
632                             if(!fmixing)        ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
633                             nOfSBCandidates++;
634                             EventHasDStarSideBandCandidate = kTRUE;
635                             }
636                             
637                         } // end if using DStar sidebans
638                 
639                 
640             }// end good D0
641            
642             // cout << "Task debug check 6 " << endl;
643             // Sideband D0
644               if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
645                 isInDZeroSideBand = kTRUE;
646                  if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
647                      if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
648                      
649                      if(fBkgMethod == kDZeroSB){
650                          if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow)       {
651                          
652                              ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
653                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
654                              if(!fmixing)       ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
655                              nOfSBCandidates++;
656                              EventHasDZeroSideBandCandidate = kTRUE;
657                          }
658                      }
659                  }
660             }// end SideBand D0
661             // cout << "Task debug check 7 " << endl;
662             
663             if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
664             if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.151) continue; // skip all D* canidates outside the interesting mass ranges
665            // cout << "Good D* candidate" << endl;
666
667             // get charge of soft pion
668             daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
669             
670             // get ID of daughters used for reconstruction
671             trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
672             trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
673             trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
674         }// end if reconstruction
675         
676         
677         
678           // ============================================ using MC kinematics only
679         Int_t DStarLabel = -1;
680         
681         if(!fReco){ // use pure MC information
682             
683             // get the DStar Particle
684             DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
685             if (!DStarMC) {
686                 AliWarning("Careful: DStar MC Particle not found in tree, skipping");
687                 continue;
688             }
689             DStarLabel = DStarMC->GetLabel();
690             if(DStarLabel>0)isDStarMCtag = kTRUE;
691             
692             Int_t PDG =TMath::Abs(DStarMC->PdgCode());
693             if(PDG !=413) continue; // skip if it is not a DStar
694             
695             // check fiducial acceptance
696             if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
697             ptbin=fCuts->PtBin(DStarMC->Pt());
698             
699             //check if DStar from B
700             Int_t labelMother = DStarMC->GetMother();
701             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
702             if(!mother) continue;
703             Int_t motherPDG =TMath::Abs(mother->PdgCode());
704             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
705             
706             
707             Bool_t isDZero = kFALSE;
708             Bool_t isSoftPi = kFALSE;
709             
710             if(fUseHadronicChannelAtKineLevel){
711                 //check decay channel on MC ************************************************
712                 Int_t NDaugh = DStarMC->GetNDaughters();
713                 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
714                 
715                 for(Int_t i=0; i<NDaugh;i++){ // loop on DStar daughters
716                     Int_t daugh_label = DStarMC->GetDaughter(i);
717                     AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
718                     if(!mcDaughter) continue;
719                     Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
720                     if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
721                     
722                     if(daugh_pdg == 421) {
723                         Int_t NDaughD0 = mcDaughter->GetNDaughters();
724                         if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
725                         trackiddaugh0 = mcDaughter->GetDaughter(0);
726                         trackiddaugh1 = mcDaughter->GetDaughter(1);
727                         Bool_t isKaon = kFALSE;
728                         Bool_t isPion = kFALSE;
729                         
730                         for(Int_t k=0;k<NDaughD0;k++){
731                             Int_t labelD0daugh = mcDaughter->GetDaughter(k);
732                             AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
733                             if(!mcGrandDaughter) continue;
734                             Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
735                             if(granddaugh_pdg==321) isKaon = kTRUE;
736                             if(granddaugh_pdg==211) isPion = kTRUE;
737                         }
738                         if(!isKaon || !isPion) break; // skip if not correct decay channel of D0
739                         isDZero = kTRUE;
740                     } // end check on Dzero
741                     
742                     if(daugh_pdg == 211) {
743                         isSoftPi = kTRUE;
744                         daughtercharge = mcDaughter->Charge();
745                         trackidsoftPi = daugh_label;}
746                 } // end loop on D* daughters
747                 if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
748             } // end check decay channel
749             
750             ptDStar = DStarMC->Pt();
751             phiDStar = DStarMC->Phi();
752             etaDStar = DStarMC->Eta();
753             
754             if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
755             
756             
757         }// end if pure MC information
758         
759         
760         
761         
762         
763         // getting the number of triggers in the MCtag D* case
764         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar);
765         if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
766         if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
767         
768         
769         fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
770         fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
771         
772         
773         
774         
775          // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
776         
777         cout << "Correlating " << endl;
778         
779         Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
780         
781         // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
782         if(fmixing && !execPool) {
783             AliInfo("Mixed event analysis: pool is not ready");
784             if(!isEventMixingFilledPeak && isInPeak)  {
785                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
786                 isEventMixingFilledPeak = kTRUE;
787             }
788             if (!isEventMixingFilledSB && isInDZeroSideBand)  {
789                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
790                 isEventMixingFilledSB=kTRUE;
791             }
792             continue;
793         } // end if pool not ok
794         // if analysis is on mixed event and pool settings are  satisfied, fill relevant histograms and continue
795         if(fmixing&&execPool){
796             // pool is ready - run checks on bins filling
797             if(!isEventMixingFilledPeak && isInPeak)  {
798                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
799                 if(fFullmode) EventMixingChecks(aodEvent);
800                 isEventMixingFilledPeak = kTRUE;
801             }
802             
803             if(!isEventMixingFilledSB && isInDZeroSideBand) {
804                 ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
805                 isEventMixingFilledSB=kTRUE;
806             }
807         } // end if pool ok
808         
809         
810         
811         
812         Int_t NofEventsinPool = 1;
813         if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
814         
815         //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
816         
817         for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
818             
819             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); // process the associated tracks
820             if(!analyzetracks) {
821                 AliInfo("AliHFCorrelator::Cannot process the track array");
822                 continue;
823             }
824             
825             Double_t arraytofill[4];
826           // Double_t MCarraytofill[4]; // think about this
827             Double_t weight;
828             
829               Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
830             
831              //************************************************** LOOP ON TRACKS *****************************************************
832             
833              for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
834                  
835                  Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
836                  if(!runcorrelation) continue;
837                  
838                  Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
839                  Double_t DeltaEta = fCorrelator->GetDeltaEta();
840                  
841                  AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
842                  if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
843                  
844                  Double_t ptHad = hadron->Pt();
845                  Double_t phiHad = hadron->Phi();
846                     phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
847                  Double_t etaHad = hadron->Eta();
848                  Int_t label = hadron->GetLabel();
849                  Int_t trackid = hadron->GetID();
850                  Double_t efficiency = hadron->GetWeight();
851                  
852                  
853                  
854                   
855                  if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
856                      
857                         ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
858                         ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
859                      isTrackForPeakFilled =  kTRUE; // makes sure you do not fill twice in case of more candidates
860                  }
861                  
862                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
863                     ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
864                     ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
865                      isTrackForSBFilled = kTRUE;
866                  }
867                  
868                  if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
869                      ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
870                      ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
871                      isTrackForSBFilled = kTRUE;
872                  }
873                  
874                  
875                  weight = 1;
876                  if(fUseEfficiencyCorrection && efficiency){
877                      weight = DmesonWeight * (1./efficiency);
878                  }
879                  
880                  
881                  arraytofill[0] = DeltaPhi;
882                  arraytofill[1] = deltainvMDStar;
883                  arraytofill[2] = DeltaEta;
884                  arraytofill[3] = ptHad;
885                 // arraytofill[4] = poolbin;
886                  
887                  
888                  // skip the D daughters in the correlation
889               //   Bool_t isDdaughter = kFALSE;
890                  if(!fmixing && fReco){ // for reconstruction
891                      if(trackid == trackiddaugh0) continue;
892                      if(trackid == trackiddaugh1) continue;
893                      if(trackid == trackidsoftPi) continue;
894                  }
895                  
896                  if(!fmixing && !fReco){ // for kinematic MC
897                      AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
898                      if(!part) continue;
899                      if(IsDDaughter(DStarMC, part)) continue;
900                      cout << "Skipping DStar  daugheter " << endl;
901                  }
902                  
903                  
904                  
905                  // ============================================= FILL CORRELATION THNSparses ===============================
906                  
907                  // filling signal
908                  if(isInPeak){
909                      cout << "Filling signal " << endl;
910                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
911                      //cout ("CorrelationsDStarHadron_%d",ptbin)
912                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
913                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
914                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
915                  }
916                  
917                   // filling bkg
918                  if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
919                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
920                       cout << "Filling bkg from D* sidebands " << endl;
921                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
922                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
923                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
924                      
925                  } // end if bkg from DStar
926                  
927                  if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
928                    //  if(!fReco && TMath::Abs(etaHad)>0.8) continue;
929                      cout << "Filling bkg from Dzero sidebands " << endl;
930                      if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
931                      if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
932                      if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
933                      
934                  } // end if bkg from DZero
935                  
936                                                                         
937                  
938              } // end loop on associated tracks
939             
940         } // end loop on events in event pool
941         
942         
943     } // end loop on D mesons
944     
945     
946     
947     
948      Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
949     
950     
951     if(!updated) AliInfo("Pool was not updated");
952     
953     
954 } //end the exec
955
956 //________________________________________ terminate ___________________________
957 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
958 {    
959         // The Terminate() function is the last function to be called during
960         // a query. It always runs on the client, it can be used to present
961         // the results graphically or save the results to file.
962         
963         AliAnalysisTaskSE::Terminate();
964         
965         fOutput = dynamic_cast<TList*> (GetOutputData(1));
966         if (!fOutput) {     
967                 printf("ERROR: fOutput not available\n");
968                 return;
969         }
970
971         return;
972 }
973 //_____________________________________________________
974 Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
975     
976     //Daughter removal in MCKine case
977     Bool_t isDaughter = kFALSE;
978     Int_t labelD0 = d->GetLabel();
979     
980     Int_t mother = track->GetMother();
981     
982     //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
983     while (mother > 0){
984         AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
985         if (mcMoth){
986             if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
987             mother = mcMoth->GetMother(); //goes back by one
988         } else{
989             AliError("Failed casting the mother particle!");
990             break;
991         }
992     }
993     
994     return isDaughter;
995 }
996
997 //_____________________________________________________
998 void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
999     
1000     Double_t Pi = TMath::Pi();
1001         Int_t nbinscorr = fPhiBins;
1002         Double_t lowcorrbin = -0.5*Pi ;
1003     Double_t upcorrbin = 1.5*Pi ;
1004     
1005     
1006     // create ThNSparses
1007     
1008     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1009     
1010     
1011     //sparse bins
1012     
1013     //1 delta_phi
1014     //2 invariant mass D *
1015     //3 delta eta
1016     //4 track pt
1017     
1018     
1019     //Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
1020     
1021     
1022     // for reconstruction on Data and MC
1023     Int_t nbinsSparse[4]=         {nbinscorr ,     32 ,  32, 10};
1024     Double_t binLowLimitSparse[4]={lowcorrbin,0.14314 ,-1.6,  0};
1025     Double_t binUpLimitSparse[4]= {upcorrbin ,0.14794 , 1.6,  5};
1026     
1027     Int_t nbinsSparseDStarSB[4]=         {nbinscorr ,     40 ,  32, 10};
1028     Double_t binLowLimitSparseDStarSB[4]={lowcorrbin,0.14788 ,-1.6,  0};
1029     Double_t binUpLimitSparseDStarSB[4]= {upcorrbin ,0.1504 , 1.6,  5};
1030     
1031     TString signalSparseName = "";
1032     TString bkgSparseName = "";
1033     
1034     THnSparseF * CorrelationsSignal = NULL;
1035     THnSparseF * CorrelationsBkg = NULL;
1036     
1037     
1038     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1039     
1040   
1041     
1042     
1043     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1044         
1045         
1046         
1047         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1048
1049
1050         
1051         signalSparseName = "CorrelationsDStar";
1052         if(fselect==1) signalSparseName += "Hadron_";
1053         if(fselect==2) signalSparseName += "Kaon_";
1054         if(fselect==3) signalSparseName += "KZero_";
1055         
1056         bkgSparseName = "CorrelationsBkg";
1057         if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
1058         if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
1059         if(fselect==1) bkgSparseName += "Hadron_";
1060         if(fselect==2) bkgSparseName += "Kaon_";
1061         if(fselect==3) bkgSparseName += "KZero_";
1062         
1063         signalSparseName+=Form("%d",iBin);
1064         bkgSparseName+=Form("%d",iBin);
1065         cout << "ThNSparses name = " << signalSparseName << endl;
1066         
1067         // define thnsparses for signal candidates
1068         CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1069           CorrelationsSignal->Sumw2();
1070         fCorrelationOutput->Add(CorrelationsSignal);
1071         
1072         // define thnsparses for bkg candidates from DStar
1073         if(fBkgMethod == kDStarSB ){
1074         CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
1075         CorrelationsBkg->Sumw2();
1076         fCorrelationOutput->Add(CorrelationsBkg);
1077         }
1078         
1079         // define thnsparses for bkg candidates from DZero
1080         if(fBkgMethod == kDZeroSB ){
1081             CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass;  #Delta #eta;AssocTrk p_{T}",4,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1082             CorrelationsBkg->Sumw2();
1083             fCorrelationOutput->Add(CorrelationsBkg);
1084         }
1085         
1086     } // end loop on bins
1087     
1088     
1089     
1090 }
1091
1092 //__________________________________________________________________________________________________
1093 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
1094     
1095     Double_t Pi = TMath::Pi();
1096         Int_t nbinscorr = fPhiBins;
1097         Double_t lowcorrbin = -0.5*Pi ;
1098     Double_t upcorrbin = 1.5*Pi ;
1099     
1100     // event counter
1101     TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
1102     NofEvents->GetXaxis()->SetBinLabel(1," All events");
1103         NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1104         NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1105         NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1106         NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1107         NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1108         NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1109         NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1110     NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1111     NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1112     NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1113     fOutput->Add(NofEvents);
1114     
1115     //event properties
1116     TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1117     fOutput->Add(EventPropsCheck);
1118     
1119     //event properties
1120     TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",1000,0,1000);
1121     fOutput->Add(SPDMultiplicty);
1122     
1123     
1124     // ===================================================   D* histograms
1125     TH1F * D0mass = NULL;
1126     TH1F * DStarMass = NULL;
1127     TH1F * DStarFromSBMass = NULL;
1128     
1129     TH1F * D0massWeighted = NULL;
1130     TH1F * DStarMassWeighted = NULL;
1131     TH1F * DStarFromSBMassWeighted = NULL;
1132     
1133     TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "";
1134     
1135     Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
1136     Float_t * ptbinlims = fCuts->GetPtBinLimits();
1137     
1138     //GetMinPtCandidate()
1139
1140     
1141     for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
1142         
1143       
1144         
1145         if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
1146         
1147         
1148         std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
1149         
1150         nameDZeroMass = "histDZeroMass_";
1151         nameDStarMass = "histDStarMass_";
1152         nameDStarFromSBMass = "histDStarFromSBMass_";
1153         
1154         nameDZeroMass+=Form("%d",iBin);
1155         nameDStarMass+=Form("%d",iBin);
1156         nameDStarFromSBMass+=Form("%d",iBin);
1157         
1158         cout << "D zero histogram: " << nameDZeroMass << endl;
1159         
1160         D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1161         DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1162         DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1163         
1164         if(!fmixing){
1165         fDmesonOutput->Add(D0mass);
1166         fDmesonOutput->Add(DStarMass);
1167         fDmesonOutput->Add(DStarFromSBMass);
1168         }
1169         
1170         // if using D meson efficiency, define weighted histos
1171         if(fUseDmesonEfficiencyCorrection){
1172            
1173             nameDZeroMass = "histDZeroMassWeight_";
1174             nameDStarMass = "histDStarMassWeight_";
1175             nameDStarFromSBMass = "histDStarFromSBMassWeight_";
1176             
1177             nameDZeroMass+=Form("%d",iBin);
1178             nameDStarMass+=Form("%d",iBin);
1179             nameDStarFromSBMass+=Form("%d",iBin);
1180             
1181             D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invarians mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
1182             DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invarians mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1183             DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invarians mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
1184             
1185             if(!fmixing){
1186             fDmesonOutput->Add(D0massWeighted);
1187             fDmesonOutput->Add(DStarMassWeighted);
1188             fDmesonOutput->Add(DStarFromSBMassWeighted);
1189             }
1190         }
1191     }// end loop on pt bins
1192     
1193     
1194     // pt distributions
1195     TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
1196     fDmesonOutput->Add(RecoPtDStar);
1197         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
1198     fDmesonOutput->Add(RecoPtBkg);
1199     
1200     // phi distribution
1201     TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1202     TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,30,0,30);
1203     
1204     // eta distribution
1205     TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1206     TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,30,0,30);
1207     
1208     if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
1209     if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
1210     if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
1211     if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
1212     
1213     
1214     // single track related histograms
1215     // phi distribution
1216     TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1217     TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
1218     
1219     // eta distribution
1220     TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1221     TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
1222     
1223     if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
1224     if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
1225     if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
1226     if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
1227
1228     
1229     // Montecarlo for D*
1230     TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1231     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
1232     
1233     TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1234     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
1235         
1236         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1237         if(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
1238     
1239     
1240     // event mixing histograms
1241     TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1242     CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1243         CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1244     CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1245         CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1246         
1247     if(fmixing) fEMOutput->Add(CheckPoolReadiness);
1248     
1249     
1250    /* Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
1251     Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
1252     Int_t nPoolBins = NofCentBins*NofZVrtxBins;
1253         
1254     
1255     TH1D * PoolBinDistribution  = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1256     fEMOutput->Add(PoolBinDistribution);
1257     
1258     TH2D * EventDistributionPerPoolBin  = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins5,-0.5,nPoolBins-0.5);
1259     fEMOutput->Add(EventDistributionPerPoolBin);
1260     */
1261 }
1262
1263
1264 //__________________________________________________________________________________________________
1265 void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1266     
1267
1268   //Float_t* ptbins = fCuts->GetPtBinLimits();
1269   if(fD0Window) delete fD0Window;
1270   fD0Window = new Float_t[fNofPtBins];
1271     
1272   AliInfo("Enlarging the D0 mass windows from cut object\n"); 
1273   Int_t nvars = fCuts->GetNVars();
1274
1275   if(nvars<1){
1276     AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1277     return;
1278   }
1279   Float_t** rdcutsvalmine=new Float_t*[nvars];
1280   for(Int_t iv=0;iv<nvars;iv++){
1281     rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1282   }
1283     
1284     
1285     for (Int_t k=0;k<nvars;k++){
1286       for (Int_t j=0;j<fNofPtBins;j++){
1287             
1288         // enlarge D0 window
1289         if(k==0)        {
1290           fD0Window[j] =fCuts->GetCutValue(0,j);
1291           rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1292           cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1293         }
1294         else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1295                         
1296         // set same windows
1297         //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1298       }
1299     }
1300     
1301     fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1302     
1303     AliInfo("\n New windows set\n");     
1304     fCuts->PrintAll();
1305     
1306     
1307     for(Int_t iv=0;iv<nvars;iv++){
1308       delete rdcutsvalmine[iv];
1309     }
1310     delete [] rdcutsvalmine;
1311     
1312 }
1313
1314
1315 //____________________________  Run checks on event mixing ___________________________________________________
1316 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
1317         // check this function
1318     
1319         AliCentrality *centralityObj = 0;
1320         Int_t multiplicity = -1;
1321         Double_t MultipOrCent = -1;
1322         
1323         // get the pool for event mixing
1324         if(fSystem != AA){ // pp
1325                 multiplicity = AOD->GetNTracks();
1326                 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1327         }
1328         if(fSystem == AA){ // PbPb
1329                 
1330                 centralityObj = AOD->GetHeader()->GetCentralityP();
1331                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1332                 AliInfo(Form("Centrality is %f", MultipOrCent));
1333         }
1334         
1335         AliAODVertex *vtx = AOD->GetPrimaryVertex();
1336         Double_t zvertex = vtx->GetZ(); // zvertex
1337         
1338         
1339         
1340         
1341         AliEventPool * pool = fCorrelator->GetPool();
1342         
1343     
1344         
1345         
1346         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1347         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1348         
1349         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1350         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1351 }
1352         
1353
1354
1355
1356