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