Usage of multi-vertexer to tag pileup events in D meson analysis (Zaida)
[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$ */
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 fReco(kTRUE),
77 fUseEfficiencyCorrection(kFALSE),
78 fUseDmesonEfficiencyCorrection(kFALSE),
79 fUseCentrality(kFALSE),
80 fUseHadronicChannelAtKineLevel(kFALSE),
81 fPhiBins(32),
82 fEvents(0),
83 fDebugLevel(0),
84 fDisplacement(0),
85 fDim(0),
86 fNofPtBins(0),
87 fMaxEtaDStar(0.9),
88 fDMesonSigmas(0),
89
90 fD0Window(0),
91
92 fOutput(0x0),
93 fOutputMC(0x0),
94 fCuts(0),
95 fCuts2(0),
96 fUtils(0),
97 fTracklets(0),
98 fDeffMapvsPt(0),
99 fDeffMapvsPtvsMult(0),
100 fDeffMapvsPtvsEta(0)
101 {
102     SetDim();
103     // default constructor
104 }
105
106 //__________________________________________________________________________
107 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :
108 AliAnalysisTaskSE(name),
109
110 fhandler(0x0),
111 fmcArray(0x0),
112 fCounter(0x0),
113 fCorrelator(0x0),
114 fselect(0),
115 fmontecarlo(kFALSE),
116 fmixing(kFALSE),
117 fFullmode(mode),
118 fSystem(syst),
119 fEfficiencyVariable(kNone),
120 fReco(kTRUE),
121 fUseEfficiencyCorrection(kFALSE),
122 fUseDmesonEfficiencyCorrection(kFALSE),
123 fUseCentrality(kFALSE),
124 fUseHadronicChannelAtKineLevel(kFALSE),
125 fPhiBins(32),
126 fEvents(0),
127 fDebugLevel(0),
128 fDisplacement(0),
129 fDim(0),
130 fNofPtBins(0),
131 fMaxEtaDStar(0.9),
132 fDMesonSigmas(0),
133 fD0Window(0),
134
135 fOutput(0x0),
136 fOutputMC(0x0),
137 fCuts(0),
138 fCuts2(AsscCuts),
139 fUtils(0),
140 fTracklets(0),
141 fDeffMapvsPt(0),
142 fDeffMapvsPtvsMult(0),
143 fDeffMapvsPtvsEta(0)
144 {
145      Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
146   
147     SetDim();
148     if(fSystem == AA)  fUseCentrality = kTRUE; else fUseCentrality = kFALSE;
149   
150     fCuts=cuts;
151     fNofPtBins= fCuts->GetNPtBins();
152     //cout << "Enlarging the DZero window " << endl;
153     EnlargeDZeroMassWindow();
154    // cout << "Done" << endl;
155     
156    
157   DefineInput(0, TChain::Class());
158   
159   DefineOutput(1,TList::Class()); // histos from data and MC
160   DefineOutput(2,TList::Class()); // histos from MC
161   DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
162   DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
163   DefineOutput(5,AliNormalizationCounter::Class());   // normalization
164 }
165
166 //__________________________________________________________________________
167
168 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
169   //
170         // destructor
171         //
172         
173         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  
174         
175         if(fhandler) {delete fhandler; fhandler = 0;}
176         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
177         if(fmcArray) {delete fmcArray; fmcArray = 0;}
178         if(fCounter) {delete fCounter; fCounter = 0;}
179         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
180         if(fOutput) {delete fOutput; fOutput = 0;}
181         if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
182         if(fCuts) {delete fCuts; fCuts = 0;}
183         if(fCuts2) {delete fCuts2; fCuts2=0;}
184     if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}
185     if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}
186     if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}
187
188 }
189
190 //___________________________________________________________
191 void AliAnalysisTaskDStarCorrelations::Init(){
192         //
193         // Initialization
194         //
195         if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
196         
197         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
198         
199         
200     
201         // Post the D* cuts
202         PostData(3,copyfCuts);
203         
204         // Post the hadron cuts
205         PostData(4,fCuts2);
206     
207         return;
208 }
209
210
211 //_________________________________________________
212 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
213         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
214         
215         //slot #1
216         //OpenFile(0);
217         fOutput = new TList();
218         fOutput->SetOwner();
219         
220         fOutputMC = new TList();
221         fOutputMC->SetOwner();
222         
223         // define histograms
224         DefineHistoForAnalysis();
225     DefineThNSparseForAnalysis();
226     
227
228         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
229         fCounter->Init();
230         
231     Double_t Pi = TMath::Pi();
232         fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb
233         fCorrelator->SetDeltaPhiInterval(  -0.5*Pi, 1.5*Pi); // set correct phi interval
234         //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval
235         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
236         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
237         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
238         fCorrelator->SetUseMC(fmontecarlo);
239         fCorrelator->SetUseReco(fReco);
240     //  fCorrelator->SetKinkRemoval(kTRUE);
241         Bool_t pooldef = fCorrelator->DefineEventPool();
242         
243         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
244     
245     fUtils = new AliAnalysisUtils();
246     
247     
248         
249         PostData(1,fOutput); // set the outputs
250         PostData(2,fOutputMC); // set the outputs
251         PostData(5,fCounter); // set the outputs
252 }
253
254 //_________________________________________________
255 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
256   
257   
258   if(fDebugLevel){
259     
260     if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;
261     if(!fReco) std::cout << "USING MC TRUTH" << std::endl;
262     std::cout << " " << std::endl;
263     std::cout << "=================================================================================" << std::endl;
264     if(!fmixing){
265       if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;
266       if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;
267       if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;
268     }
269     if(fmixing){
270       if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;
271       if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;
272       if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;
273     }
274     
275   }// end if debug
276   
277   
278   
279   if (!fInputEvent) {
280     Error("UserExec","NO EVENT FOUND!");
281                 return;
282   }
283   
284   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
285   if(!aodEvent){
286     AliError("AOD event not found!");
287     return;
288   }
289   
290     fTracklets = aodEvent->GetTracklets();
291     
292     fEvents++; // event counter
293     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
294     
295     fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);
296     
297     // load MC array
298     fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
299     if(fmontecarlo && !fmcArray){
300       AliError("Array of MC particles not found");
301       return;
302     }
303     
304     
305     
306     
307     // ********************************************** START EVENT SELECTION ****************************************************
308     
309     Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
310     
311     if(!isEvSel) {
312       
313       if(fCuts->IsEventRejectedDueToPileup()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);
314       if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);
315       if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);
316       if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);
317       if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);
318       if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);
319       if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);
320       
321       return;
322     }
323     
324     // added event selection for pA
325     
326     if(fSystem == pA){
327       
328       if(fUtils->IsFirstEventInChunk(aodEvent)) {
329         AliInfo("Rejecting the event - first in the chunk");
330         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);
331         return;
332         }
333       if(!fUtils->IsVertexSelected2013pA(aodEvent)) {
334         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);
335         AliInfo("Rejecting the event - bad vertex");
336         return;
337       }
338     }
339     // ******************************** END event selections **************************************************
340     
341     AliCentrality *centralityObj = 0;
342     Double_t MultipOrCent = -1;
343     
344     if(fUseCentrality){
345       /* if(fSystem == AA ){    */      centralityObj = aodEvent->GetHeader()->GetCentralityP();
346       MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
347       //AliInfo(Form("Centrality is %f", MultipOrCent));
348     }
349     
350     if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.);
351     
352         
353     fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
354     
355     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
356     
357     Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings
358         if(!correlatorON) {
359           AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
360           return;
361         }
362         
363         if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);
364                 
365         // check the event type
366         // load MC header
367         if(fmontecarlo){
368           AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
369           if (fmontecarlo && !mcHeader) {
370             AliError("Could not find MC Header in AOD");
371             return;
372           }
373           
374           Bool_t isMCeventgood = kFALSE;
375           
376           
377           Int_t eventType = mcHeader->GetEventType();
378           Int_t NMCevents = fCuts2->GetNofMCEventType();
379           
380           for(Int_t k=0; k<NMCevents; k++){
381             Int_t * MCEventType = fCuts2->GetMCEventType();
382             
383             if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
384             ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
385           }
386           
387           if(NMCevents && !isMCeventgood){
388             if(fDebugLevel)     std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
389             return;
390           }
391           
392         } // end if montecarlo
393         
394         
395         // checks on vertex and multiplicity of the event
396         AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
397         Double_t zVtxPosition = vtx->GetZ(); // zvertex
398         
399         
400         if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);
401         
402         
403         
404         // D* reconstruction
405         TClonesArray *arrayDStartoD0pi=0;
406         if(!aodEvent && AODEvent() && IsStandardAOD()) {
407           // In case there is an AOD handler writing a standard AOD, use the AOD
408           // event in memory rather than the input (ESD) event.
409           aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
410           // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
411           // have to taken from the AOD event hold by the AliAODExtension
412           AliAODHandler* aodHandler = (AliAODHandler*)
413             ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
414           if(aodHandler->GetExtensions()) {
415             AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
416             AliAODEvent *aodFromExt = ext->GetAOD();
417             arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
418           }
419         } else {
420           arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
421         }
422         
423         if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
424         
425     
426     // get the poolbin
427     
428     Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition);
429   
430     
431         // initialize variables you will need for the D*
432         
433         Double_t ptDStar;//
434         Double_t phiDStar;//
435         Double_t etaDStar;//
436         Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag;
437         Double_t invMassDZero;
438         Double_t deltainvMDStar;
439     
440         
441         Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
442         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
443         
444         
445         //MC tagging for DStar
446         //D* and D0 prongs needed to MatchToMC method
447         Int_t pdgDgDStartoD0pi[2]={421,211};
448         Int_t pdgDgD0toKpi[2]={321,211};
449         
450         Bool_t isDStarCand = kFALSE;
451     Bool_t isDfromB = kFALSE;
452         Bool_t isEventMixingFilledPeak = kFALSE;
453         Bool_t isEventMixingFilledSB = kFALSE;
454     Bool_t EventHasDStarCandidate = kFALSE;
455     Bool_t EventHasDZeroSideBandCandidate = kFALSE;
456     Bool_t EventHasDStarSideBandCandidate = kFALSE;
457         //loop on D* candidates
458         
459         Int_t looponDCands = 0;
460         if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast();
461         if(!fReco) looponDCands = fmcArray->GetEntriesFast();
462         
463         Int_t nOfDStarCandidates = 0;
464         Int_t nOfSBCandidates = 0;
465         
466         Double_t DmesonEfficiency = 1.;
467         Double_t DmesonWeight = 1.;
468     Double_t efficiencyvariable = -999;
469     
470         
471         
472         for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
473           isInPeak = kFALSE;
474           isInDStarSideBand = kFALSE;
475           isInDZeroSideBand = kFALSE;
476           isDStarMCtag = kFALSE;
477           isDfromB = kFALSE;
478           ptDStar = -123.4;
479           phiDStar = -999;
480           etaDStar = -56.;
481           invMassDZero = - 999;
482           deltainvMDStar = -998;
483           AliAODRecoCascadeHF* dstarD0pi;
484           AliAODRecoDecayHF2Prong* theD0particle;
485           AliAODMCParticle* DStarMC=0x0;
486       Short_t daughtercharge = -2;
487           Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
488           Int_t trackiddaugh1 = -1;
489           Int_t trackidsoftPi = -1;
490           
491           // start the if reconstructed candidates from here ************************************************
492           
493           if(fReco){//// if reconstruction is applied
494             dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
495             if(!dstarD0pi->GetSecondaryVtx()) continue;
496             theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
497             if (!theD0particle) continue;
498             
499             
500             // track quality cuts
501             Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
502             // region of interest + topological cuts + PID
503             Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
504             
505             //apply track selections
506             if(!isTkSelected) continue;
507             if(!isSelected) continue;
508         if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
509           
510           
511           ptDStar = dstarD0pi->Pt();
512           phiDStar = dstarD0pi->Phi();
513           etaDStar = dstarD0pi->Eta();
514           if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
515           if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;
516           if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;
517           if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();
518           if(fEfficiencyVariable == kNone) efficiencyvariable = 0;
519           
520         // get the D meson efficiency
521         DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
522             
523          
524
525             if(fUseDmesonEfficiencyCorrection){
526               if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
527               else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI
528                 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
529                 else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low
530               }
531             }
532             else DmesonWeight = 1.; 
533             
534             // continue;
535           
536             Int_t mcLabelDStar = -999;
537           if(fmontecarlo){
538               // find associated MC particle for D* ->D0toKpi
539               mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);
540               if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
541               if(!isDStarMCtag) continue;
542             AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
543             //check if DStar from B
544             Int_t labelMother = MCDStar->GetMother();
545             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
546             if(!mother) continue;
547             Int_t motherPDG =TMath::Abs(mother->PdgCode());
548             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
549             }
550                         
551             
552             phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar);
553             
554             // set the phi of the D meson in the correct range
555             
556             Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
557             
558             
559             Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma
560             //  Double_t mD0Window=0.074/3;
561             
562             Double_t mD0Window= fD0Window[ptbin]/3;
563             //cout << "Check with new window " << fD0Window[ptbin]/3 << endl;
564             
565             
566             
567             invMassDZero = dstarD0pi->InvMassD0();
568             if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
569             
570             deltainvMDStar = dstarD0pi->DeltaInvMass();
571             
572             //good D0 candidates
573             if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
574               
575               if(!fmixing)      ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);
576               // good D*?
577               if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow){
578                 
579                 if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight);
580                 if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);
581                 isInPeak = kTRUE;
582         EventHasDStarCandidate = kTRUE;
583                 nOfDStarCandidates++;
584               } // end Good D*
585               
586                 //  D* sideband?
587               if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)<fDMesonSigmas[3]*dmDStarWindow)){
588                 isInDStarSideBand = kTRUE;
589               EventHasDStarSideBandCandidate = kTRUE;
590               } // end D* sideband
591               
592             }// end good D0 candidates
593             
594             //D0 sidebands
595             if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
596               if(!fmixing)((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);
597               if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight);
598               
599               if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0] *dmDStarWindow){ // is in DStar peak region?
600                 if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight);
601                 isInDZeroSideBand = kTRUE;
602         EventHasDZeroSideBandCandidate = kTRUE;
603                 nOfSBCandidates++;
604                 if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);
605               }
606               
607             }//end if sidebands
608             
609            
610             
611             
612             if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
613             
614             
615             // check properties of the events containing the D*
616
617      
618             
619             isDStarCand = kTRUE;
620             
621             // charge of the daughter od the
622             daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
623             
624             
625             trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
626             trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
627             trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
628             
629             // end here the reco
630             
631             
632           }// end of if for applied reconstruction to D*
633           
634           Int_t DStarLabel = -1;
635           
636           if(!fReco){ // use pure MC information
637           
638         // get the DStar Particle
639             DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
640             if (!DStarMC) {
641               AliWarning("Careful: DStar MC Particle not found in tree, skipping");
642               continue;
643             }
644             DStarLabel = DStarMC->GetLabel();
645         if(DStarLabel>0)isDStarMCtag = kTRUE;
646             
647             Int_t PDG =TMath::Abs(DStarMC->PdgCode());
648             if(PDG !=413) continue; // skip if it is not a DStar
649             // check fiducial acceptance
650         if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
651             
652             //check if DStar from B
653             Int_t labelMother = DStarMC->GetMother();
654             AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
655                  if(!mother) continue;
656             Int_t motherPDG =TMath::Abs(mother->PdgCode());
657             if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
658           
659           Bool_t isDZero = kFALSE;
660           Bool_t isSoftPi = kFALSE;
661           
662           if(fUseHadronicChannelAtKineLevel){
663           //check decay channel on MC ************************************************
664                 Int_t NDaugh = DStarMC->GetNDaughters();
665                 if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs
666                 
667                 for(Int_t i=0; i<NDaugh;i++){ // loop on daughters
668                         Int_t daugh_label = DStarMC->GetDaughter(i);
669                         AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));
670                         if(!mcDaughter) continue;
671                         Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());
672                         if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;
673               
674                         if(daugh_pdg == 421) {isDZero = kTRUE;
675                             Int_t NDaughD0 = mcDaughter->GetNDaughters();
676                             if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs
677                             trackiddaugh0 = mcDaughter->GetDaughter(0);
678                             trackiddaugh1 = mcDaughter->GetDaughter(1);
679                             Bool_t isKaon = kFALSE;
680                             Bool_t isPion = kFALSE;
681                 
682                             for(Int_t k=0;k<NDaughD0;k++){
683                                 Int_t labelD0daugh = mcDaughter->GetDaughter(k);
684                                 AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));
685                                 if(!mcGrandDaughter) continue;
686                                 Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());
687                                 if(granddaugh_pdg==321) isKaon = kTRUE;
688                                 if(granddaugh_pdg==211) isPion = kTRUE;
689                             }
690                             if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0
691                         }
692               
693                         if(daugh_pdg == 211) {
694                             isSoftPi = kTRUE;
695                             daughtercharge = mcDaughter->Charge();
696                             trackidsoftPi = daugh_label;}
697                     }
698               if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
699           } // end check decay channel
700             
701             ptDStar = DStarMC->Pt();
702             phiDStar = DStarMC->Phi();
703             etaDStar = DStarMC->Eta();
704           
705           if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
706             
707           } // end use pure MC information
708           
709     
710         // getting the number of triggers in the MCtag D* case
711         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);
712           if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);
713           if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);
714           
715           
716           fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
717           fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
718           
719           
720           // ************************************************ CORRELATION ANALYSIS STARTS HERE
721           
722           
723           Bool_t execPool = fCorrelator->ProcessEventPool();
724           
725           
726           if(fmixing && !execPool) {
727             AliInfo("Mixed event analysis: pool is not ready");
728             if(!isEventMixingFilledPeak && isInPeak)  {
729               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1);
730               isEventMixingFilledPeak = kTRUE;
731             }
732             if (!isEventMixingFilledSB && isInDZeroSideBand)  {
733               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3);
734               isEventMixingFilledSB=kTRUE;
735             }
736             
737             continue;
738           }
739           
740           // check event topology
741           if(fmixing&&execPool){
742             // pool is ready - run checks on bins filling
743             if(!isEventMixingFilledPeak && isInPeak)  {
744               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0);
745               if(fFullmode) EventMixingChecks(aodEvent);
746               isEventMixingFilledPeak = kTRUE;
747             }
748             
749             if(!isEventMixingFilledSB && isInDZeroSideBand) {
750               ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2);
751               isEventMixingFilledSB=kTRUE;
752             }
753           }
754           
755           Int_t NofEventsinPool = 1;
756           if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
757           
758           
759           for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
760             
761             Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
762             if(!analyzetracks) {
763               AliInfo("AliHFCorrelator::Cannot process the track array");
764               continue;
765             }
766             
767             //initialization of variables for correlations with leading particles
768             Double_t DeltaPhiLeading = -999.;
769           Double_t DeltaEtaLeading = -999.;
770         
771             
772             Int_t NofTracks = fCorrelator->GetNofTracks();
773             
774             
775             if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks);
776             if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks);
777             
778             
779             
780             Double_t arraytofill[5];
781             Double_t MCarraytofill[7];
782             
783             
784             Double_t weight;
785             
786           for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
787               Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
788               if(!runcorrelation) continue;
789               
790               Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
791               Double_t DeltaEta = fCorrelator->GetDeltaEta();
792               
793               AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
794               if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
795               
796               Double_t ptHad = hadron->Pt();
797               Double_t phiHad = hadron->Phi();
798               Double_t etaHad = hadron->Eta();
799               Int_t label = hadron->GetLabel();
800               Int_t trackid = hadron->GetID();
801               Double_t efficiency = hadron->GetWeight();
802               
803               weight = 1;
804               if(fUseEfficiencyCorrection && efficiency){
805                   weight = DmesonWeight * (1./efficiency);
806               }
807         
808               phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
809               
810               
811               if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency);
812               
813               arraytofill[0] = DeltaPhi;
814               arraytofill[1] = DeltaEta;
815               arraytofill[2] = ptDStar;
816               arraytofill[3] = ptHad;
817               arraytofill[4] = poolbin;
818                 
819               
820               MCarraytofill[0] = DeltaPhi;
821               MCarraytofill[1] = DeltaEta;
822               MCarraytofill[2] = ptDStar;
823               MCarraytofill[3] = ptHad;
824               MCarraytofill[4] = poolbin;
825               
826              
827               if(fmontecarlo){
828                 if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks);
829                 if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks);
830                 if(label<0) continue; // skip track with wrong label
831               }
832               
833                Bool_t isDdaughter = kFALSE;
834             // skip the D daughters in the correlation
835               if(!fmixing && fReco){
836                   if(trackid == trackiddaugh0) continue;
837                   if(trackid == trackiddaugh1) continue;
838                   if(trackid == trackidsoftPi) continue;
839               }
840               
841               if(!fmixing && !fReco){
842                   AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
843                   if(!part) continue;
844                   if(IsDDaughter(DStarMC, part)) continue;
845                   cout << "Skipping DStar  daugheter " << endl;
846               }
847               if(!fmixing && !fReco && fmontecarlo){  // skip D* Daughetrs if it is Pure MCDStar
848                     Int_t hadronlabel = label;
849                     for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers
850                         if(DStarLabel<0){ break;}
851                         if(hadronlabel<0) { break;}
852                         AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fmcArray->At(hadronlabel));
853                         if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;}
854                         hadronlabel = mcParticle->GetMother();
855                         if(hadronlabel == DStarLabel) isDdaughter = kTRUE;
856                     }
857                 
858                   if(isDdaughter && fDebugLevel){
859                       std::cout << "It is the D* daughter with label " << label << std::endl;
860                       std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl;
861                       std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl;
862                       std::cout << "Soft pi label = " << trackidsoftPi << std::endl;
863                   }
864                     
865                                         if(isDdaughter) continue; // skip if track is from DStar
866                                 }
867                 
868                                 // ================ FILL CORRELATION HISTOGRAMS ===============================
869                                 
870                 // monte carlo case (mc tagged D*)
871               if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo
872                 
873                 Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
874                               
875                 MCarraytofill[5] = 0;
876                 if(PartSource[0]) MCarraytofill[5] = 1;
877                 if(PartSource[1]) MCarraytofill[5] = 2;
878                 if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3;
879                     if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4;
880                     if(PartSource[3]) MCarraytofill[5] = 5;
881                     if(!isDfromB) MCarraytofill[6] = 0;
882                     if(isDfromB) MCarraytofill[6] = 1;
883                     if(!fReco && TMath::Abs(etaHad)>0.8) {
884                       delete [] PartSource;
885                       continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance
886                     }
887                     ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill);
888                     
889                     delete[] PartSource;
890                                 }
891               
892                 // Good DStar canidates
893                                 if(isInPeak)  {
894                     
895                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance
896                     if(fselect==1)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight);
897                     if(fselect==2)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight);
898                     if(fselect==3)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight);
899                                         
900                                     ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent);
901                                         if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad);
902                                     
903                                 }
904               
905                 // Sidebands from D0 candidate
906                                 if(isInDZeroSideBand) {
907                                         
908                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance
909                     if(fselect==1)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);
910                     if(fselect==2)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);
911                     if(fselect==3)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);
912                     
913                                         if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad);
914                     
915                 }
916               
917               // Sidebands from D* candidate
918                 if(isInDStarSideBand) {
919                                         
920                                         if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance
921                     if(fselect==1 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);
922                     if(fselect==2 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);
923                     if(fselect==3 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);
924
925                                 }
926                 
927                 
928                         } // end loop on track candidates
929             
930                         
931                         
932             // fill the leading particle histograms
933             
934             if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
935                         if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
936             
937                 } // end loop on events in the pool
938         
939         }// end loop on D* candidates
940         
941     
942  
943         // check events with D* or SB canidates
944     if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition);
945      if(fFullmode && EventHasDZeroSideBandCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition);
946     
947     if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition);
948     
949     
950         if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent);
951     if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent);
952         
953     // update event pool
954     Bool_t updated = fCorrelator->PoolUpdate();
955         
956     //  if(updated) EventMixingChecks(aodEvent);
957         if(!updated) AliInfo("Pool was not updated");
958         
959     
960 } //end the exec
961
962 //________________________________________ terminate ___________________________
963 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
964 {    
965         // The Terminate() function is the last function to be called during
966         // a query. It always runs on the client, it can be used to present
967         // the results graphically or save the results to file.
968         
969         AliAnalysisTaskSE::Terminate();
970         
971         fOutput = dynamic_cast<TList*> (GetOutputData(1));
972         if (!fOutput) {     
973                 printf("ERROR: fOutput not available\n");
974                 return;
975         }
976
977         return;
978 }
979 //_____________________________________________________
980 Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {
981     
982     //Daughter removal in MCKine case
983     Bool_t isDaughter = kFALSE;
984     Int_t labelD0 = d->GetLabel();
985     
986     Int_t mother = track->GetMother();
987     
988     //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
989     while (mother > 0){
990         AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
991         if (mcMoth){
992             if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
993             mother = mcMoth->GetMother(); //goes back by one
994         } else{
995             AliError("Failed casting the mother particle!");
996             break;
997         }
998     }
999     
1000     return isDaughter;
1001 }
1002
1003 //_____________________________________________________
1004 void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
1005     
1006     //cout << "DEFINING THNSPARSES "<< endl;
1007     
1008     Double_t Pi = TMath::Pi();
1009         Int_t nbinscorr = fPhiBins;
1010         Double_t lowcorrbin = -0.5*Pi;
1011         Double_t upcorrbin = 1.5*Pi;
1012     // define the THnSparseF
1013     
1014     //sparse bins
1015     
1016     //1 delta_phi
1017     //2 delta_eta
1018     //3 D* pt
1019     //4 multiplicity
1020     //5 track pt
1021     //6 zVtx position
1022     
1023     Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins());
1024     
1025     
1026     Int_t nbinsSparse[5]=         {nbinscorr,   32,30,250,nbinsPool};
1027     Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0,  0,-0.5};
1028     Double_t binUpLimitSparse[5]= {upcorrbin,  1.6,30,25,nbinsPool-0.5};
1029   
1030     Int_t MCnbinsSparse[7]=         {nbinscorr,   32,30,250,nbinsPool,10,2};    
1031     Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0,  0,-0.5,-0.5,-0.5};      //  
1032     Double_t MCbinUpLimitSparse[7]= {upcorrbin,  1.6,30,25,nbinsPool-0.5,9.5,1.5};
1033     
1034     TString sparsename = "CorrelationsDStar";
1035     if(fselect==1) sparsename += "Hadron";
1036         if(fselect==2) sparsename += "Kaon";
1037         if(fselect==3) sparsename += "KZero";
1038     
1039     TString D0Bkgsparsename = "DZeroBkg";
1040     D0Bkgsparsename += sparsename;
1041     
1042     TString DStarBkgsparsename = "DStarBkg";
1043     DStarBkgsparsename += sparsename;
1044     
1045     TString MCSparseName = "MCDStar";
1046     MCSparseName += sparsename;
1047     // signal correlations
1048     THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1049     
1050     // bkg correlations from D0 sidebands
1051     THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1052     
1053     // bkg correlations from D* sidebands
1054     THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
1055     
1056     
1057     THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse);
1058     
1059     MCCorrelations->GetAxis(5)->SetBinLabel(1," All ");
1060         MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour");
1061         MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D");
1062         MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D");
1063         MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B");
1064         MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour");
1065         MCCorrelations->GetAxis(5)->SetBinLabel(7," from c");
1066         MCCorrelations->GetAxis(5)->SetBinLabel(8," from b");
1067     
1068     MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c");
1069     MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b");
1070     
1071     Correlations->Sumw2();
1072     DZeroBkgCorrelations->Sumw2();
1073     DStarBkgCorrelations->Sumw2();
1074     
1075     fOutput->Add(Correlations);
1076     fOutput->Add(DZeroBkgCorrelations);
1077     if(fFullmode) fOutput->Add(DStarBkgCorrelations);
1078     if(fmontecarlo) fOutputMC->Add(MCCorrelations);
1079     
1080     
1081 }
1082 //__________________________________________________________________________________________________
1083 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
1084         
1085         Double_t Pi = TMath::Pi();
1086         Int_t nbinscorr = fPhiBins;
1087         Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center
1088         Double_t upcorrbin = 1.5*Pi ;
1089         
1090         // ========================= histograms for both Data and MonteCarlo
1091         
1092         
1093         TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);
1094     NofEvents->GetXaxis()->SetBinLabel(1," All events");
1095         NofEvents->GetXaxis()->SetBinLabel(2," Selected events");
1096         NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");
1097         NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");
1098         NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");
1099         NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");
1100         NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");
1101         NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");
1102     NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");
1103     NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");
1104     NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");
1105     fOutput->Add(NofEvents);
1106         
1107         
1108         
1109         
1110         TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
1111         if(!fmixing && fFullmode) fOutput->Add(D0InvMass);
1112         
1113         TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);
1114         if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB);
1115         
1116         //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
1117         //if(!fmixing) fOutput->Add(DeltaInvMass);
1118     TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution; D* p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);
1119         if(!fmixing) fOutput->Add(DeltaInvMass);
1120         
1121         TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution; SB p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);
1122         if(!fmixing) fOutput->Add(bkgDeltaInvMass);
1123     
1124     DeltaInvMass->Sumw2();
1125     bkgDeltaInvMass->Sumw2();
1126         
1127         TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);
1128         if(!fmixing) fOutput->Add(RecoPtDStar);
1129         
1130         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);
1131         if(!fmixing) fOutput->Add(RecoPtBkg);
1132     
1133     TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);
1134     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);
1135     
1136     TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);
1137     if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);
1138         
1139         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
1140         if(!fmixing) fOutput->Add(MCtagPtDStar);
1141         
1142         TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
1143     if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra);
1144         
1145         TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);
1146         if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF);
1147         
1148         TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5);
1149         if(fFullmode) fOutput->Add(NofTracksInPeak);
1150         
1151         TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5);
1152         if(fFullmode) fOutput->Add(NofTracksInSB);
1153         
1154         TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);
1155         if(fFullmode)fOutput->Add(TracksInPeakSpectra);
1156         
1157         TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);
1158         if(fFullmode)fOutput->Add(TracksInSBSpectra);
1159         
1160         
1161         //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
1162         //if(fmixing) fOutput->Add(EventMixingCheck);
1163     
1164     
1165     TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1166         if(fFullmode)fOutput->Add(EventPropsCheck);
1167     
1168     TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1169         if(fFullmode)fOutput->Add(EventPropsCheckifDStar);
1170     
1171     TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1172         if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB);
1173     
1174     TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);
1175         if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB);
1176     
1177         
1178     TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005);
1179         if(fFullmode)fOutput->Add(WeightChecks);
1180     
1181         
1182         
1183         TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);
1184         fOutput->Add(PhiEtaTrigger);
1185         
1186         TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);
1187         fOutput->Add(PhiEtaSideBand);
1188         
1189         TH3F * PhiEtaPart = new TH3F("PhiEtaPart","#phi distribution of the associated particle; #phi; #eta; multiplicity",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9,100,0,1000);
1190         fOutput->Add(PhiEtaPart);
1191     
1192     TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);
1193         if(fFullmode)fOutput->Add(DStarCandidates);
1194     
1195     TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);
1196         if(fFullmode)fOutput->Add(SBCandidates);
1197         
1198         
1199         //correlations histograms
1200         TString histoname1 = "DPhiDStar";
1201         if(fselect==1) histoname1 += "Hadron";
1202         if(fselect==2) histoname1 += "Kaon";
1203         if(fselect==3) histoname1 += "KZero";
1204         
1205         /*
1206      TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1207      
1208      TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1209      
1210      //side band background histograms
1211      TString histoname2 = "bkg";
1212      histoname2 += histoname1;
1213      TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1214      TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1215      
1216      
1217      fOutput->Add(DPhiDStar);
1218      
1219      if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
1220      
1221      fOutput->Add(bkgDPhiDStar);
1222      
1223      if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
1224      
1225      */
1226         // leading particle
1227         TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1228         TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1229         
1230         if(fFullmode)fOutput->Add(leadingcand);
1231         if(fFullmode)fOutput->Add(leadingsidebands);
1232         
1233         // ========================= histos for analysis on MC only
1234         
1235         TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1236         if(fmontecarlo) fOutputMC->Add(EventTypeMC);
1237         
1238         TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
1239         MCSources->GetXaxis()->SetBinLabel(1," All ");
1240         MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");
1241         MCSources->GetXaxis()->SetBinLabel(3," from c->D");
1242         MCSources->GetXaxis()->SetBinLabel(4," from b->D");
1243         MCSources->GetXaxis()->SetBinLabel(5," from b->B");
1244         MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
1245         MCSources->GetXaxis()->SetBinLabel(7," from c");
1246         MCSources->GetXaxis()->SetBinLabel(8," from b");
1247         
1248         if(fmontecarlo) fOutputMC->Add(MCSources);
1249     
1250     // leading particle from mc source
1251     TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);
1252         LeadingMCSources->GetXaxis()->SetBinLabel(1," All ");
1253         LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");
1254         LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");
1255         LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");
1256         LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");
1257         LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
1258         LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");
1259         LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");
1260         
1261         if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources);
1262         
1263     // all hadrons
1264         TString histoname3 = "MCTag";
1265         histoname3 += histoname1;
1266         TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1267         
1268         TString histoname44 = "CharmDOrigin";
1269         histoname44 += histoname1;
1270         histoname44 += "MC";
1271         
1272         TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1273         
1274         
1275         TString histoname54 = "BeautyDOrigin";
1276         histoname54 += histoname1;
1277         histoname54 += "MC";
1278         TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1279         
1280         TString histoname55 = "BeautyBOrigin";
1281         histoname55 += histoname1;
1282         histoname55 += "MC";
1283         TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1284         
1285         TString histoname4 = "CharmQuarkOrigin";
1286         histoname4 += histoname1;
1287         histoname4 += "MC";
1288         TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1289         
1290         TString histoname5 = "BeautyQuarkOrigin";
1291         histoname5 += histoname1;
1292         histoname5 += "MC";
1293         TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1294         
1295         TString histoname6 = "NonHFOrigin";
1296         histoname6 += histoname1;
1297         histoname6 += "MC";
1298         TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1299     
1300     if(fmontecarlo && fFullmode){
1301         
1302         fOutputMC->Add(MCTagDPhiDStar);
1303         fOutputMC->Add(CharmDOriginDPhiDStar);
1304         fOutputMC->Add(BeautyDOriginDPhiDStar);
1305         fOutputMC->Add(BeautyBOriginDPhiDStar);
1306         fOutputMC->Add(CharmQuarkOriginDPhiDStar);
1307         fOutputMC->Add(BeautyQuarkOriginDPhiDStar);
1308                 fOutputMC->Add(NonHFOriginDPhiDStar);
1309         
1310         }
1311     
1312     // ========================= histos for analysis on MC
1313     // all leading hadron
1314         TString Leadinghistoname3 = "LeadingMCTag";
1315         Leadinghistoname3 += histoname1;
1316         TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1317     
1318         TString Leadinghistoname44 = "LeadingCharmDOrigin";
1319         Leadinghistoname44 += histoname1;
1320         Leadinghistoname44 += "MC";
1321         
1322         TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1323         
1324         
1325         TString Leadinghistoname54 = "LeadingBeautyDOrigin";
1326         Leadinghistoname54 += histoname1;
1327         Leadinghistoname54 += "MC";
1328         TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1329         
1330         TString Leadinghistoname55 = "LeadingBeautyBOrigin";
1331         Leadinghistoname55 += histoname1;
1332         Leadinghistoname55 += "MC";
1333         TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1334         
1335         TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";
1336         Leadinghistoname4 += histoname1;
1337         Leadinghistoname4 += "MC";
1338         TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1339         
1340         TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";
1341         Leadinghistoname5 += histoname1;
1342         Leadinghistoname5 += "MC";
1343         TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);
1344     
1345     
1346         
1347         
1348         if(fmontecarlo && fFullmode){
1349                 
1350                 fOutputMC->Add(LeadingMCTagDPhiDStar);
1351                 fOutputMC->Add(LeadingCharmDOriginDPhiDStar);
1352                 fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);
1353                 fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);
1354                 fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);
1355                 fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);
1356                 
1357         }
1358         
1359         TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5);
1360         MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");
1361         MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");
1362         MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");
1363         MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");
1364         MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");
1365         MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");
1366         if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);
1367         
1368         TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5);
1369         if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels);
1370         
1371         // ============================= EVENT MIXING CHECKS ======================================
1372         
1373         Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();
1374         Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();
1375         Int_t NofCentBins = fCuts2->GetNCentPoolBins();
1376         Double_t * CentBins = fCuts2->GetCentPoolBins();
1377         Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();
1378         Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();
1379     
1380     
1381     
1382         Int_t k =0;
1383         
1384         if(fSystem == AA) k = 100; // PbPb centrality
1385         if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity
1386         
1387         
1388         //Double_t minvalue = CentBins[0];
1389         //Double_t maxvalue = CentBins[NofCentBins+1];
1390         //Double_t Zminvalue = ZVrtxBins[0];
1391         //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];
1392         
1393         Double_t minvalue, maxvalue;
1394     Double_t Zminvalue, Zmaxvalue;
1395     
1396     Zminvalue = -15.;
1397     Zmaxvalue = 15;
1398     if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb
1399     if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity
1400     
1401         //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};
1402     // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};
1403     //  Double_t * events = Nevents;
1404     Double_t eventsv[] ={0,1000000};
1405     //Double_t * events = new Double_t[2];
1406    // events[0] = 0;
1407 //      events[1] = 1000000;
1408     Double_t *events = eventsv;
1409     Int_t Nevents = 1000000;
1410     //  TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events);
1411     
1412     TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]);
1413     
1414         EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
1415         EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
1416         EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");
1417         if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin);
1418         
1419         Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;
1420         //Int_t Diff = MaxNofTracks-MinNofTracks;
1421         
1422     //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};
1423     //  Double_t  * trackN = Ntracks;
1424         
1425         TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks);
1426         NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
1427         NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
1428         NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");
1429         
1430         if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin);
1431         
1432         TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);
1433         NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");
1434         NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");
1435         if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls);
1436         
1437     
1438         
1439         TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);
1440         EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");
1441         EventProps->GetYaxis()->SetTitle("Z vertex [cm]");
1442         if(fmixing && fFullmode) fOutput->Add(EventProps);
1443     
1444     TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);
1445     CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");
1446         CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");
1447     CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");
1448         CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");
1449         
1450     if(fmixing) fOutput->Add(CheckPoolReadiness);
1451     
1452         
1453 }
1454
1455 //__________________________________________________________________________________________________
1456 void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){
1457     
1458
1459   //Float_t* ptbins = fCuts->GetPtBinLimits();
1460   if(fD0Window) delete fD0Window;
1461   fD0Window = new Float_t[fNofPtBins];
1462     
1463   AliInfo("Enlarging the D0 mass windows from cut object\n"); 
1464   Int_t nvars = fCuts->GetNVars();
1465
1466   if(nvars<1){
1467     AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");
1468     return;
1469   }
1470   Float_t** rdcutsvalmine=new Float_t*[nvars];
1471   for(Int_t iv=0;iv<nvars;iv++){
1472     rdcutsvalmine[iv]=new Float_t[fNofPtBins];
1473   }
1474     
1475     
1476     for (Int_t k=0;k<nvars;k++){
1477       for (Int_t j=0;j<fNofPtBins;j++){
1478             
1479         // enlarge D0 window
1480         if(k==0)        {
1481           fD0Window[j] =fCuts->GetCutValue(0,j);
1482           rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);
1483           cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;
1484         }
1485         else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);
1486                         
1487         // set same windows
1488         //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);
1489       }
1490     }
1491     
1492     fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);
1493     
1494     AliInfo("\n New windows set\n");     
1495     fCuts->PrintAll();
1496     
1497     
1498     for(Int_t iv=0;iv<nvars;iv++){
1499       delete rdcutsvalmine[iv];
1500     }
1501     delete [] rdcutsvalmine;
1502     
1503 }
1504
1505
1506 //____________________________  Run checks on event mixing ___________________________________________________
1507 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
1508         
1509     
1510         AliCentrality *centralityObj = 0;
1511         Int_t multiplicity = -1;
1512         Double_t MultipOrCent = -1;
1513         
1514         // get the pool for event mixing
1515         if(fSystem != AA){ // pp
1516                 multiplicity = AOD->GetNTracks();
1517                 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1518         }
1519         if(fSystem == AA){ // PbPb
1520                 
1521                 centralityObj = AOD->GetHeader()->GetCentralityP();
1522                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1523                 AliInfo(Form("Centrality is %f", MultipOrCent));
1524         }
1525         
1526         AliAODVertex *vtx = AOD->GetPrimaryVertex();
1527         Double_t zvertex = vtx->GetZ(); // zvertex
1528         
1529         
1530         
1531         
1532         AliEventPool * pool = fCorrelator->GetPool();
1533         
1534     
1535         
1536         
1537         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1538         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1539         
1540         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool
1541         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool
1542 }
1543         
1544
1545
1546
1547