Fixes for std:: need with the trunk of Root (Jochen, Yves)
[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
53 using std::cout;
54 using std::endl;
55
56
57 ClassImp(AliAnalysisTaskDStarCorrelations)
58
59
60 //__________________________________________________________________________
61 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
62 AliAnalysisTaskSE(),
63 fhandler(0x0),
64 //fPoolMgr(0x0),      
65 fmcArray(0x0),
66 fCounter(0x0),
67 fCorrelator(0x0),
68 fselect(0),
69 fmontecarlo(kFALSE),
70 fmixing(kFALSE),
71 fSystem(kFALSE),
72 fEvents(0),
73 fDebug(0),
74 fDisplacement(0),
75
76 fOutput(0x0),            
77 fCuts(0),
78 fCuts2(0)
79 {
80 // default constructor  
81 }
82
83 //__________________________________________________________________________
84 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts) :
85 AliAnalysisTaskSE(name),
86
87 fhandler(0x0),
88 fmcArray(0x0),
89 fCounter(0x0),
90 fCorrelator(0x0),
91 fselect(0),
92 fmontecarlo(kFALSE),
93 fmixing(kFALSE),
94 fSystem(kFALSE),
95 fEvents(0),
96 fDebug(0),
97 fDisplacement(0),
98
99 fOutput(0x0),                
100 fCuts(0),
101 fCuts2(AsscCuts)
102 {
103         fCuts=cuts;
104         Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
105         DefineInput(0, TChain::Class());
106         DefineOutput(1,TList::Class()); // histos from data
107         DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
108         DefineOutput(3,AliNormalizationCounter::Class());   // normalization
109         DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my cuts
110 }
111
112 //__________________________________________________________________________
113
114 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
115         //
116         // destructor
117         //
118         
119         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  
120         
121         if(fhandler) {delete fhandler; fhandler = 0;}
122         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
123         if(fmcArray) {delete fmcArray; fmcArray = 0;}
124         if(fCounter) {delete fCounter; fCounter = 0;}
125         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
126         if(fOutput) {delete fOutput; fOutput = 0;}
127         if(fCuts) {delete fCuts; fCuts = 0;}
128         if(fCuts2) {delete fCuts2; fCuts2=0;}
129
130 }
131
132 //___________________________________________________________
133 void AliAnalysisTaskDStarCorrelations::Init(){
134         //
135         // Initialization
136         //
137         if(fDebug > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
138         
139         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
140         
141         
142         
143         
144         // Post the D* cuts
145         PostData(2,copyfCuts);
146         
147         // Post the hadron cuts
148         PostData(4,fCuts2);
149         
150
151         
152         return;
153 }
154
155
156 //_________________________________________________
157 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
158         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
159         
160         //slot #1  
161         //OpenFile(0);
162         fOutput = new TList();
163         fOutput->SetOwner();
164         
165         // define histograms
166         DefineHistoForAnalysis();
167         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(3)->GetContainer()->GetName()));
168         fCounter->Init();
169         
170     Double_t Pi = TMath::Pi();
171         fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fSystem); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb
172         fCorrelator->SetDeltaPhiInterval((-0.5-1./32)*Pi,(1.5-1./32)*Pi); // set correct phi interval
173         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
174         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
175         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
176         fCorrelator->SetUseMC(fmontecarlo);
177         Bool_t pooldef = fCorrelator->DefineEventPool();
178         
179         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
180
181         
182         PostData(1,fOutput); // set the outputs
183         PostData(3,fCounter); // set the outputs
184 }
185 //_________________________________________________
186 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
187
188         
189         
190         cout << " " << endl;
191         cout << "=================================================================================" << endl;
192         if(!fmixing){
193         if(fselect==1) cout << "TASK::Correlation with hadrons on SE "<< endl;
194         if(fselect==2) cout << "TASK::Correlation with kaons on SE "<< endl;
195         if(fselect==3) cout << "TASK::Correlation with kzeros on SE "<< endl;
196         }
197         if(fmixing){
198                 if(fselect==1) cout << "TASK::Correlation with hadrons on ME "<< endl;
199                 if(fselect==2) cout << "TASK::Correlation with kaons on ME "<< endl;
200                 if(fselect==3) cout << "TASK::Correlation with kzeros on ME "<< endl;
201         }
202         if (!fInputEvent) {
203                 Error("UserExec","NO EVENT FOUND!");
204                 return;
205         }
206         
207         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
208         if(!aodEvent){
209           AliError("AOD event not found!");
210           return;
211         }
212         
213         fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
214         
215         fEvents++; // event counter
216         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
217         fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
218         
219         if(fmontecarlo && !fmcArray){
220           AliError("Array of MC particles not found");
221           return;
222         }
223         Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
224         if(!isEvSel) return;
225         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
226         //
227         Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
228         if(!correlatorON) {
229                 AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
230                 return;
231         }
232         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);
233         
234         if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);
235         // D* reconstruction
236         TClonesArray *arrayDStartoD0pi=0;
237
238         
239         if(!aodEvent && AODEvent() && IsStandardAOD()) {
240                 // In case there is an AOD handler writing a standard AOD, use the AOD 
241                 // event in memory rather than the input (ESD) event.    
242                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
243                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
244                 // have to taken from the AOD event hold by the AliAODExtension
245                 AliAODHandler* aodHandler = (AliAODHandler*) 
246                 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
247                 if(aodHandler->GetExtensions()) {
248                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
249                         AliAODEvent *aodFromExt = ext->GetAOD();
250                         arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
251                 }
252         } else {
253                 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
254         }
255         
256         if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
257         
258         // initialize variables you will need for the D*
259         
260         Double_t ptDStar;//
261         Double_t phiDStar;//
262         Double_t etaDStar;//
263         Bool_t isInPeak, isInSideBand, isDStarMCtag;
264         Double_t invMassDZero;
265         Double_t deltainvMDStar;
266
267         
268         Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
269         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
270                 
271         
272         //MC tagging for DStar
273         //D* and D0 prongs needed to MatchToMC method
274         Int_t pdgDgDStartoD0pi[2]={421,211};
275         Int_t pdgDgD0toKpi[2]={321,211};
276
277         
278         //loop on D* candidates
279         for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
280                 isInPeak = kFALSE;
281                 isInSideBand = kFALSE;
282                 isDStarMCtag = kFALSE;
283                 ptDStar = -123.4;
284                 phiDStar = -999;
285                 etaDStar = -56.;
286                 invMassDZero = - 999;
287                 deltainvMDStar = -998;
288                 
289                 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
290                 if(!dstarD0pi->GetSecondaryVtx()) continue;
291                 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
292                 if (!theD0particle) continue;
293
294                 
295                 // track quality cuts
296                 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
297                 // region of interest + topological cuts + PID
298                 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
299                 //apply selections
300                 if(!isTkSelected) continue;
301                 if(!isSelected) continue;
302                 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;    
303                 Int_t mcLabelDStar = -999;
304                 if(fmontecarlo){
305                         // find associated MC particle for D* ->D0toKpi
306                         mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray,kFALSE);
307                         if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
308                 }
309                 
310                 ptDStar = dstarD0pi->Pt();
311                 phiDStar = dstarD0pi->Phi(); 
312                 etaDStar = dstarD0pi->Eta();
313                 
314                 phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the phi of the D meson in the correct range
315                 
316                 Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
317                 
318                 Double_t dmDStarWindow =0.0019;// 0.0019 = 3 sigma
319                 Double_t mD0Window=0.074;
320                 
321                 if (!fSystem){ // pp
322                         if (ptbin==1)  mD0Window = 0.026; //0.5-1
323                         if (ptbin==2)  mD0Window = 0.022; //1-2
324                         if (ptbin==3)  mD0Window = 0.024; //2-3
325                         if (ptbin==4)  mD0Window = 0.032;
326                         if (ptbin==5)  mD0Window = 0.032;
327                         if (ptbin==6)  mD0Window = 0.036;
328                         if (ptbin==7)  mD0Window = 0.036;
329                         if (ptbin==8)  mD0Window = 0.036;
330                         if (ptbin==9)  mD0Window = 0.058;
331                         if (ptbin==10) mD0Window = 0.058;
332                         if (ptbin>10)  mD0Window = 0.074;
333                 }
334                 if(fSystem){// PbPb
335                         if (ptbin==0)  mD0Window = 0.032; //1-1
336                         if (ptbin==1)  mD0Window = 0.032; //2-3
337                         if (ptbin==2)  mD0Window = 0.032; //3-4
338                         if (ptbin==3)  mD0Window = 0.032; //4-5
339                         if (ptbin==4)  mD0Window = 0.036; //5-6
340                         if (ptbin==5)  mD0Window = 0.036; //6-8
341                         if (ptbin==6)  mD0Window = 0.055; //8-12
342                         if (ptbin==7)  mD0Window = 0.074; //12-16
343                         if (ptbin==8)  mD0Window = 0.074; //16-24
344                         if (ptbin==9)  mD0Window = 0.074; //24-35 
345                  }
346                 
347                 invMassDZero = dstarD0pi->InvMassD0();
348                 ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
349                 
350                 deltainvMDStar = dstarD0pi->DeltaInvMass();
351                 
352         
353                 //good candidates
354                 if (TMath::Abs(invMassDZero-mPDGD0)<mD0Window){
355                         
356                         ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
357                         if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
358                                 
359                                 ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar);
360                                 isInPeak = kTRUE;
361                                 ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);
362                         }
363                 }// end if good candidates
364                 
365                 //sidebands
366                 if (TMath::Abs(invMassDZero-mPDGD0)>1.3*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<4.*mD0Window ){
367                         ((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
368                         ((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero);
369                         
370                         if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
371                                 ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar);
372                                 isInSideBand = kTRUE;
373                                 ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);
374                         }
375                         
376                 }//end if sidebands
377                 // getting the number of triggers in the MCtag D* case 
378                 
379
380         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);
381                 
382                 
383                 if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
384                 
385                 
386             fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
387                 
388                 Short_t daughtercharge = ((AliAODTrack*)theD0particle->GetDaughter(0))->Charge();
389                 fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
390                 
391                 
392                 Int_t trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
393                 Int_t trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
394                 Int_t trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
395                 
396                 Bool_t execPool = fCorrelator->ProcessEventPool();
397                 if(fmixing && !execPool) {
398                         AliInfo("Mixed event analysis: pool is not ready");
399                         continue;
400                 }
401                 
402                 Int_t NofEventsinPool = 1;
403                 if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
404                                 
405                 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
406                 
407                         Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
408                         
409                         if(!analyzetracks) {
410                                 AliInfo("AliHFCorrelator::Cannot process the track array");
411                                 continue;
412                         }
413                 
414                 
415                         Int_t NofTracks = fCorrelator->GetNofTracks();
416                 
417                         for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
418                         
419                                 Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
420                                 if(!runcorrelation) continue;
421                         
422                                 Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
423                                 Double_t DeltaEta = fCorrelator->GetDeltaEta();
424                         
425                                 AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
426                                 
427                                 Double_t ptHad = hadron->Pt();
428                                 Double_t phiHad = hadron->Phi();
429                                 Double_t etaHad = hadron->Eta(); 
430                                 Double_t label = hadron->GetLabel(); 
431                                 Double_t trackid = hadron->GetID();
432                                 
433                                 phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
434                                 
435                                 if(!fmixing){ // skip D* Daughetrs
436                                         if(trackid == trackiddaugh0) continue;
437                                         if(trackid == trackiddaugh1) continue;
438                                         if(trackid == trackidsoftPi) continue;
439                                 }
440                         
441                                 // from here on it is up to the user to decide what object to fill
442                                 
443                                 if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
444                                 
445                                         Int_t PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
446                                 
447                                         FillMCTagCorrelations(ptDStar,DeltaPhi,DeltaEta,ptHad,PartSource);
448                                 
449                                 }
450                         
451                                 if(isInPeak)  {
452                                 
453                                         if(fselect==1) ((TH3D*)fOutput->FindObject("DPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
454                                         if(fselect==2) ((TH3D*)fOutput->FindObject("DPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
455                                         if(fselect==3) ((TH3D*)fOutput->FindObject("DPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
456                                     ((TH2F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad);
457                                         //counterPeak++; // count tracks per peak per event
458                                 
459                                 }
460                         
461                                 if(isInSideBand) {
462                                 
463                                         if(fselect==1) ((TH3D*)fOutput->FindObject("bkgDPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
464                                         if(fselect==2) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
465                                         if(fselect==3) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
466                                 
467                                         
468                                 
469                                         //counterSB++;
470                                 }
471                         
472                         
473                         } // end loop on track candidates
474                 } // end loop on events in the pool
475                                 
476         }// end loop on D* candidates           
477         Bool_t updated = fCorrelator->PoolUpdate();
478         if(!updated) AliInfo("Pool was not updated");
479         
480                 //cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END OF THE EVENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
481         
482         
483                 
484 } //end the exec
485
486 //________________________________________ terminate ___________________________
487 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
488 {    
489         // The Terminate() function is the last function to be called during
490         // a query. It always runs on the client, it can be used to present
491         // the results graphically or save the results to file.
492         
493         AliAnalysisTaskSE::Terminate();
494         
495         fOutput = dynamic_cast<TList*> (GetOutputData(1));
496         if (!fOutput) {     
497                 printf("ERROR: fOutput not available\n");
498                 return;
499         }
500
501         return;
502 }
503
504
505 //_____________________________________________________
506 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
507         
508         Double_t Pi = TMath::Pi();
509         Int_t nbinscorr = 32;
510         Double_t lowcorrbin = -0.5*Pi - Pi/32; // shift the bin by half the width so that at 0 is it the bin center
511         Double_t upcorrbin = 1.5*Pi - Pi/32;    
512         
513         // ========================= histograms for both Data and MonteCarlo
514         
515         
516         TH1D * NofEvents = new TH1D("NofEvents","NofEvents",3,0,3);
517         fOutput->Add(NofEvents);
518                 
519         TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
520         fOutput->Add(D0InvMass);
521         
522         TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);
523         fOutput->Add(D0InvMassinSB);
524         
525         TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
526         fOutput->Add(DeltaInvMass);
527         
528         TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
529         fOutput->Add(bkgDeltaInvMass);
530         
531         TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",30,0,30);
532         fOutput->Add(RecoPtDStar);
533         
534         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",30,0,30);
535         fOutput->Add(RecoPtBkg);
536         
537         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",30,0,30);
538         fOutput->Add(MCtagPtDStar);
539         
540         TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
541         if(fselect==3) fOutput->Add(KZeroSpectra);
542         
543         TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);
544         if(fselect==3) fOutput->Add(KZeroSpectraifHF);
545         
546         TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","NofTracksInPeak",500,0.5,500.5);
547         fOutput->Add(NofTracksInPeak);
548         
549         TH1D * NofTracksInSB = new TH1D("NofTracksInSB","NofTracksInSB",500,0.5,500.5);
550         fOutput->Add(NofTracksInSB);
551         
552         TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
553         if(fmixing) fOutput->Add(EventMixingCheck);
554         
555
556         
557         TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
558         MCSources->GetXaxis()->SetBinLabel(1,"All ");
559         MCSources->GetXaxis()->SetBinLabel(2," from Heavy flavour");
560         MCSources->GetXaxis()->SetBinLabel(3," from c->D");
561         MCSources->GetXaxis()->SetBinLabel(4," from b->D");
562         MCSources->GetXaxis()->SetBinLabel(5," from b->B");
563         if(fmontecarlo) fOutput->Add(MCSources);
564         
565         TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
566         fOutput->Add(PhiEtaTrigger);
567         
568         TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
569         fOutput->Add(PhiEtaSideBand);
570         
571         TH2F * PhiEtaPart = new TH2F("PhiEtaPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
572         fOutput->Add(PhiEtaPart);
573
574
575         //correlations histograms
576         TString histoname1 = "DPhiDStar";
577         if(fselect==1) histoname1 += "Hadron";
578         if(fselect==2) histoname1 += "Kaon";
579         if(fselect==3) histoname1 += "KZero";
580
581         
582         TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
583         
584         TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
585         
586         //side band background histograms
587         TString histoname2 = "bkg";
588         histoname2 += histoname1;
589         TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
590         TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
591         
592         
593         fOutput->Add(DPhiDStar);
594
595         if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
596         
597         fOutput->Add(bkgDPhiDStar);
598
599         if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
600
601
602         // ========================= histos for analysis on MC
603         TString histoname3 = "MCTag";
604         histoname3 += histoname1;
605         TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
606          
607         TString histoname44 = "CharmDOrigin";
608         histoname44 += histoname1;
609         histoname44 += "MC";
610         
611         TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
612         
613         
614         TString histoname54 = "BeautyDOrigin";
615         histoname54 += histoname1;
616         histoname54 += "MC";
617         TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
618         
619         TString histoname55 = "BeautyBOrigin";
620         histoname55 += histoname1;
621         histoname55 += "MC";
622         TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
623         
624         if(fmontecarlo){
625         
626         fOutput->Add(MCTagDPhiDStar);
627         fOutput->Add(CharmDOriginDPhiDStar);
628         fOutput->Add(BeautyDOriginDPhiDStar);
629         fOutput->Add(BeautyBOriginDPhiDStar);
630         
631         }
632         
633
634         
635 }
636
637
638
639 //____________________________  Function for MC correlations ___________________________________________________
640 void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Int_t mcSource){
641
642
643         if(fselect==1) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
644         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKaon"))->Fill(DelPhi,ptTrig,DelEta);
645         if(fselect==3) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKZero"))->Fill(DelPhi,ptTrig,DelEta);
646
647
648
649 ((TH1F*)fOutput->FindObject("MCSources"))->Fill(0);
650
651 if (mcSource==44){ // is from charm ->D
652         if(fselect==1) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
653         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
654         if(fselect==3) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
655         
656         
657         ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
658         ((TH1F*)fOutput->FindObject("MCSources"))->Fill(2);
659         }
660
661 if (mcSource==54){ // is from beauty -> D
662         if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
663         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
664         if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
665         if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
666         if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(3);
667         }
668
669 if (mcSource==55){ // is from beauty ->B
670         if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
671         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
672         if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyOriginBDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
673         if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
674         if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(4);
675         }
676         return;
677 }
678
679
680
681
682