| 92d9c9cb |
1 | /**************************************************************************\r |
| 2 | * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r |
| 3 | * *\r |
| 4 | * Author: The ALICE Off-line Project. *\r |
| 5 | * Contributors are mentioned in the code where appropriate. *\r |
| 6 | * *\r |
| 7 | * Permission to use, copy, modify and distribute this software and its *\r |
| 8 | * documentation strictly for non-commercial purposes is hereby granted *\r |
| 9 | * without fee, provided that the above copyright notice appears in all *\r |
| 10 | * copies and that both the copyright notice and this permission notice *\r |
| 11 | * appear in the supporting documentation. The authors make no claims *\r |
| 12 | * about the suitability of this software for any purpose. It is *\r |
| 13 | * provided "as is" without express or implied warranty. *\r |
| 14 | **************************************************************************/\r |
| 15 | //\r |
| 16 | //\r |
| 17 | // Analysis Taks for reconstructed particle correlation \r |
| 18 | // (first implementation done for D mesons) with jets \r |
| 19 | // (use the so called Emcal framework)\r |
| 20 | //\r |
| 21 | //-----------------------------------------------------------------------\r |
| 22 | // Authors:\r |
| 23 | // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch\r |
| 24 | // A. Grelli (Utrecht University) a.grelli@uu.nl\r |
| 25 | // X. Zhang (LBNL) XMZhang@lbl.gov\r |
| 26 | //-----------------------------------------------------------------------\r |
| 27 | \r |
| 28 | #include <TDatabasePDG.h>\r |
| 29 | #include <TParticle.h>\r |
| 30 | #include "TROOT.h"\r |
| 31 | #include <TH3F.h>\r |
| b770538f |
32 | #include <THnSparse.h>\r |
| 92d9c9cb |
33 | \r |
| 34 | #include "AliAnalysisTaskFlavourJetCorrelations.h"\r |
| 35 | #include "AliAODMCHeader.h"\r |
| 36 | #include "AliAODHandler.h"\r |
| 37 | #include "AliAnalysisManager.h"\r |
| 38 | #include "AliLog.h"\r |
| 39 | #include "AliEmcalJet.h"\r |
| 40 | #include "AliJetContainer.h"\r |
| 41 | #include "AliAODRecoDecay.h"\r |
| 42 | #include "AliAODRecoCascadeHF.h"\r |
| 43 | #include "AliAODRecoDecayHF2Prong.h"\r |
| 44 | #include "AliESDtrack.h"\r |
| 45 | #include "AliAODMCParticle.h"\r |
| 46 | #include "AliPicoTrack.h"\r |
| 47 | #include "AliRDHFCutsD0toKpi.h"\r |
| 48 | #include "AliRDHFCutsDStartoKpipi.h"\r |
| 49 | \r |
| 50 | ClassImp(AliAnalysisTaskFlavourJetCorrelations)\r |
| 51 | \r |
| 52 | \r |
| 53 | //_______________________________________________________________________________\r |
| 54 | \r |
| 55 | AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :\r |
| 56 | AliAnalysisTaskEmcalJet("",kFALSE),\r |
| 57 | fUseMCInfo(kTRUE), \r |
| 58 | fUseReco(kTRUE),\r |
| 59 | fCandidateType(),\r |
| 60 | fPDGmother(),\r |
| 61 | fNProngs(),\r |
| 62 | fPDGdaughters(),\r |
| 63 | fBranchName(),\r |
| 64 | fmyOutput(0),\r |
| 65 | fCuts(0),\r |
| 66 | fMinMass(),\r |
| 67 | fMaxMass(), \r |
| 68 | fJetArrName(0),\r |
| 69 | fCandArrName(0),\r |
| 70 | fLeadingJetOnly(kFALSE),\r |
| 71 | fJetRadius(0)\r |
| 72 | {\r |
| 73 | //\r |
| 74 | // Default ctor\r |
| 75 | //\r |
| 76 | //SetMakeGeneralHistograms(kTRUE);\r |
| 77 | \r |
| 78 | }\r |
| 79 | \r |
| 80 | //_______________________________________________________________________________\r |
| 81 | \r |
| 82 | AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :\r |
| 83 | AliAnalysisTaskEmcalJet(name,kFALSE),\r |
| 84 | fUseMCInfo(kTRUE),\r |
| 85 | fUseReco(kTRUE), \r |
| 86 | fCandidateType(),\r |
| 87 | fPDGmother(),\r |
| 88 | fNProngs(),\r |
| 89 | fPDGdaughters(),\r |
| 90 | fBranchName(),\r |
| 91 | fmyOutput(0),\r |
| 92 | fCuts(0),\r |
| 93 | fMinMass(),\r |
| 94 | fMaxMass(), \r |
| 95 | fJetArrName(0),\r |
| 96 | fCandArrName(0),\r |
| 97 | fLeadingJetOnly(kFALSE),\r |
| 98 | fJetRadius(0)\r |
| 99 | {\r |
| 100 | //\r |
| 101 | // Constructor. Initialization of Inputs and Outputs\r |
| 102 | //\r |
| 103 | \r |
| 104 | Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");\r |
| 105 | fCuts=cuts;\r |
| 106 | fCandidateType=candtype;\r |
| 107 | const Int_t nptbins=fCuts->GetNPtBins();\r |
| 108 | Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};\r |
| 109 | \r |
| 110 | switch(fCandidateType){\r |
| 111 | case 0:\r |
| 112 | fPDGmother=421;\r |
| 113 | fNProngs=2;\r |
| 114 | fPDGdaughters[0]=211;//pi \r |
| 115 | fPDGdaughters[1]=321;//K\r |
| 116 | fPDGdaughters[2]=0; //empty\r |
| 117 | fPDGdaughters[3]=0; //empty\r |
| 118 | fBranchName="D0toKpi";\r |
| 119 | fCandArrName="D0";\r |
| 120 | break;\r |
| 121 | case 1: \r |
| 122 | fPDGmother=413;\r |
| 123 | fNProngs=3;\r |
| 124 | fPDGdaughters[1]=211;//pi soft\r |
| 125 | fPDGdaughters[0]=421;//D0\r |
| 126 | fPDGdaughters[2]=211;//pi fromD0\r |
| 127 | fPDGdaughters[3]=321; // K from D0\r |
| 128 | fBranchName="Dstar";\r |
| 129 | fCandArrName="Dstar";\r |
| 130 | \r |
| 131 | if(nptbins<=13){\r |
| 132 | for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];\r |
| 133 | } else {\r |
| 134 | AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));\r |
| 135 | }\r |
| 136 | break;\r |
| 137 | default:\r |
| 138 | printf("%d not accepted!!\n",fCandidateType);\r |
| 139 | break;\r |
| 140 | }\r |
| 141 | \r |
| 142 | if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);\r |
| 143 | if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);\r |
| 144 | \r |
| 145 | DefineOutput(1,TList::Class()); // histos\r |
| 146 | DefineOutput(2,AliRDHFCuts::Class()); // my cuts\r |
| 147 | \r |
| 148 | }\r |
| 149 | \r |
| 150 | //_______________________________________________________________________________\r |
| 151 | \r |
| 152 | AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {\r |
| 153 | //\r |
| 154 | // destructor\r |
| 155 | //\r |
| 156 | \r |
| 157 | Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor"); \r |
| 158 | \r |
| 159 | delete fmyOutput;\r |
| 160 | delete fCuts;\r |
| 161 | \r |
| 162 | }\r |
| 163 | \r |
| 164 | //_______________________________________________________________________________\r |
| 165 | \r |
| 166 | void AliAnalysisTaskFlavourJetCorrelations::Init(){\r |
| 167 | //\r |
| 168 | // Initialization\r |
| 169 | //\r |
| 170 | \r |
| 171 | if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");\r |
| 172 | switch(fCandidateType){\r |
| 173 | case 0:\r |
| 174 | {\r |
| 175 | AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));\r |
| 176 | copyfCuts->SetName("AnalysisCutsDzero");\r |
| 177 | // Post the data\r |
| 178 | PostData(2,copyfCuts);\r |
| 179 | }\r |
| 180 | break;\r |
| 181 | case 1:\r |
| 182 | {\r |
| 183 | AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));\r |
| 184 | copyfCuts->SetName("AnalysisCutsDStar");\r |
| 185 | // Post the cuts\r |
| 186 | PostData(2,copyfCuts);\r |
| 187 | }\r |
| 188 | break;\r |
| 189 | default:\r |
| 190 | return;\r |
| 191 | }\r |
| 192 | \r |
| 193 | return;\r |
| 194 | }\r |
| 195 | \r |
| 196 | //_______________________________________________________________________________\r |
| 197 | \r |
| 198 | void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { \r |
| 199 | // output \r |
| 200 | Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r |
| 201 | fmyOutput = new TList();\r |
| 202 | fmyOutput->SetOwner();\r |
| 203 | fmyOutput->SetName("pippo");\r |
| 204 | // define histograms\r |
| 205 | DefineHistoForAnalysis();\r |
| 206 | PostData(1,fmyOutput);\r |
| 207 | \r |
| 208 | return;\r |
| 209 | }\r |
| 210 | \r |
| 211 | //_______________________________________________________________________________\r |
| 212 | \r |
| 213 | void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)\r |
| 214 | {\r |
| 215 | // user exec\r |
| 216 | if (!fInitialized){\r |
| 217 | AliAnalysisTaskEmcalJet::ExecOnce();\r |
| 218 | }\r |
| 219 | \r |
| 220 | // Load the event\r |
| 221 | AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r |
| 222 | \r |
| 223 | TClonesArray *arrayDStartoD0pi=0;\r |
| 224 | TClonesArray *trackArr = 0;\r |
| 225 | TClonesArray *candidatesArr = 0;\r |
| 226 | TClonesArray *sbcandArr = 0;\r |
| 227 | \r |
| 228 | if (!aodEvent && AODEvent() && IsStandardAOD()) {\r |
| 229 | \r |
| 230 | // In case there is an AOD handler writing a standard AOD, use the AOD \r |
| 231 | // event in memory rather than the input (ESD) event. \r |
| 232 | aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r |
| 233 | \r |
| 234 | // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r |
| 235 | // have to taken from the AOD event hold by the AliAODExtension\r |
| 236 | AliAODHandler* aodHandler = (AliAODHandler*) \r |
| 237 | ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r |
| 238 | if(aodHandler->GetExtensions()) {\r |
| 239 | AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r |
| 240 | AliAODEvent *aodFromExt = ext->GetAOD();\r |
| 241 | arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());\r |
| 242 | }\r |
| 243 | } else {\r |
| 244 | arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());\r |
| 245 | }\r |
| 246 | \r |
| 247 | if (!arrayDStartoD0pi) {\r |
| 248 | AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));\r |
| 249 | // return;\r |
| 250 | } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast())); \r |
| 251 | \r |
| 252 | TClonesArray* mcArray = 0x0;\r |
| 253 | if (fUseMCInfo) {\r |
| 254 | mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r |
| 255 | if (!mcArray) {\r |
| 256 | printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");\r |
| 257 | return;\r |
| 258 | }\r |
| 259 | }\r |
| 260 | \r |
| 261 | //retrieve jets\r |
| 262 | trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));\r |
| 263 | //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));\r |
| 264 | //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));\r |
| 265 | fJetRadius=GetJetContainer(0)->GetJetRadius();\r |
| 266 | \r |
| 267 | if(!trackArr){\r |
| 268 | AliInfo(Form("Could not find the track array\n"));\r |
| 269 | return;\r |
| 270 | }\r |
| 271 | \r |
| 272 | \r |
| 273 | TString arrname="fCandidateArray";\r |
| 274 | TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r |
| 275 | candidatesArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));\r |
| 276 | if (!candidatesArr) {\r |
| 277 | Printf("%s array not found",candarrname.Data());\r |
| 278 | InputEvent()->GetList()->ls();\r |
| 279 | return;\r |
| 280 | }\r |
| 281 | //Printf("ncandidates found %d",candidatesArr->GetEntriesFast());\r |
| 282 | \r |
| 283 | TString arrSBname="fSideBandArray";\r |
| 284 | TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r |
| 285 | sbcandArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));\r |
| 286 | if (fCandidateType==1 && !sbcandArr) {\r |
| 287 | Printf("%s array not found",sbarrname.Data());\r |
| 288 | InputEvent()->GetList()->ls();\r |
| 289 | return;\r |
| 290 | }\r |
| 291 | //Printf("nSBCand or Bkg found %d",sbcandArr->GetEntriesFast());\r |
| 292 | \r |
| 293 | \r |
| 294 | //Histograms\r |
| 295 | TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");\r |
| 296 | TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");\r |
| 297 | TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");\r |
| 298 | TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");\r |
| 299 | TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");\r |
| 300 | TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");\r |
| 301 | TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");\r |
| 302 | TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");\r |
| 303 | TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");\r |
| 304 | TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");\r |
| 305 | TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");\r |
| 306 | TH1F* hdeltaRJetTracks=((TH1F*)fmyOutput->FindObject("hdeltaRJetTracks"));\r |
| 307 | \r |
| 308 | hstat->Fill(0);\r |
| 309 | \r |
| 310 | // fix for temporary bug in ESDfilter \r |
| 311 | // the AODs with null vertex pointer didn't pass the PhysSel\r |
| 312 | if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r |
| 313 | \r |
| 314 | //Event selection\r |
| 315 | Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);\r |
| 316 | TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();\r |
| 317 | if(!iseventselected) return;\r |
| 318 | \r |
| 319 | hstat->Fill(1);\r |
| 320 | \r |
| 321 | //trigger on jets\r |
| 322 | \r |
| 323 | Int_t njets=GetJetContainer()->GetNJets();\r |
| 324 | if(njets == 0) return;\r |
| 325 | \r |
| 326 | //retrieve charm candidates selected\r |
| 327 | Int_t candidates = candidatesArr->GetEntriesFast();\r |
| 328 | \r |
| 329 | // we start with jets\r |
| 330 | Double_t ejet = 0;\r |
| 331 | Double_t phiJet = 0;\r |
| 332 | Double_t etaJet = 0;\r |
| 333 | Double_t ptjet = 0;\r |
| 334 | Double_t leadingJet =0;\r |
| 335 | \r |
| 336 | Int_t ntrarr=trackArr->GetEntriesFast();\r |
| 337 | hNtrArr->Fill(ntrarr);\r |
| 338 | \r |
| 339 | for(Int_t i=0;i<ntrarr;i++){\r |
| 340 | AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));\r |
| 341 | //reusing histograms\r |
| 342 | hPtJetTrks->Fill(jtrack->Pt());\r |
| 343 | hPhiJetTrks->Fill(jtrack->Phi());\r |
| 344 | hEtaJetTrks->Fill(jtrack->Eta());\r |
| 345 | hEjetTrks->Fill(jtrack->E());\r |
| 346 | }\r |
| 347 | \r |
| 348 | \r |
| 349 | //option to use only the leading jet\r |
| 350 | if(fLeadingJetOnly){\r |
| 351 | for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) { \r |
| 352 | AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);\r |
| 353 | ptjet = jetL->Pt();\r |
| 354 | if(ptjet>leadingJet ) leadingJet = ptjet;\r |
| 355 | \r |
| 356 | }\r |
| 357 | }\r |
| 358 | \r |
| 359 | Int_t cntjet=0;\r |
| 360 | //loop over jets and consider only the leading jet in the event\r |
| 361 | for (Int_t iJets = 0; iJets<njets; iJets++) { \r |
| 362 | //Printf("Jet N %d",iJets);\r |
| 363 | AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);\r |
| 364 | //Printf("Corr task Accept Jet");\r |
| 365 | if(!AcceptJet(jet)) {\r |
| 366 | hstat->Fill(5);\r |
| 367 | continue;\r |
| 368 | }\r |
| 369 | \r |
| 370 | //jets variables\r |
| 371 | ejet = jet->E();\r |
| 372 | phiJet = jet->Phi();\r |
| 373 | etaJet = jet->Eta();\r |
| 374 | ptjet = jet->Pt();\r |
| 375 | \r |
| 376 | // choose the leading jet\r |
| 377 | if(fLeadingJetOnly && (ejet<leadingJet)) continue;\r |
| 378 | //used for normalization\r |
| 379 | hstat->Fill(3);\r |
| 380 | cntjet++;\r |
| 381 | // fill energy, eta and phi of the jet\r |
| 382 | hEjet ->Fill(ejet);\r |
| 383 | hPhiJet ->Fill(phiJet);\r |
| 384 | hEtaJet ->Fill(etaJet);\r |
| 385 | hPtJet ->Fill(ptjet);\r |
| 386 | \r |
| 387 | //loop on jet particles\r |
| 388 | Int_t ntrjet= jet->GetNumberOfTracks(); \r |
| 389 | for(Int_t itrk=0;itrk<ntrjet;itrk++){\r |
| 390 | \r |
| 391 | AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr); \r |
| 392 | hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));\r |
| 393 | \r |
| 394 | }//end loop on jet tracks\r |
| 395 | \r |
| 396 | //Printf("N candidates %d ", candidates);\r |
| 397 | for(Int_t ic = 0; ic < candidates; ic++) {\r |
| 398 | \r |
| 399 | // D* candidates\r |
| 400 | AliVParticle* charm=0x0;\r |
| 401 | charm=(AliVParticle*)candidatesArr->At(ic);\r |
| 402 | if(!charm) continue;\r |
| 403 | hstat->Fill(2);\r |
| 404 | \r |
| 405 | FlagFlavour(charm, jet);\r |
| 406 | if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);\r |
| 407 | \r |
| 408 | FillHistogramsRecoJetCorr(charm, jet);\r |
| 409 | \r |
| 410 | }\r |
| 411 | //retrieve side band background candidates for Dstar\r |
| 412 | Int_t nsbcand = 0; if(sbcandArr) nsbcand=sbcandArr->GetEntriesFast();\r |
| 413 | \r |
| 414 | for(Int_t ib=0;ib<nsbcand;ib++){\r |
| 415 | if(fCandidateType==kDstartoKpipi && !fUseMCInfo){\r |
| 416 | AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)sbcandArr->At(ib);\r |
| 417 | if(!sbcand) continue;\r |
| 418 | SideBandBackground(sbcand,jet);\r |
| 419 | }\r |
| 420 | if(fUseMCInfo){\r |
| 421 | AliAODRecoDecayHF* charmbg = 0x0;\r |
| 422 | charmbg=(AliAODRecoDecayHF*)candidatesArr->At(ib);\r |
| 423 | if(!charmbg) continue;\r |
| 424 | MCBackground(charmbg,jet);\r |
| 425 | }\r |
| 426 | }\r |
| 427 | } // end of jet loop\r |
| 428 | \r |
| 429 | hNJetPerEv->Fill(cntjet);\r |
| 430 | \r |
| 431 | PostData(1,fmyOutput);\r |
| 432 | \r |
| 433 | }\r |
| 434 | \r |
| 435 | //_______________________________________________________________________________\r |
| 436 | \r |
| 437 | void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)\r |
| 438 | { \r |
| 439 | // The Terminate() function is the last function to be called during\r |
| 440 | // a query. It always runs on the client, it can be used to present\r |
| 441 | // the results graphically or save the results to file.\r |
| 442 | \r |
| 443 | Info("Terminate"," terminate");\r |
| 444 | AliAnalysisTaskSE::Terminate();\r |
| 445 | \r |
| 446 | fmyOutput = dynamic_cast<TList*> (GetOutputData(1));\r |
| 447 | if (!fmyOutput) { \r |
| 448 | printf("ERROR: fmyOutput not available\n");\r |
| 449 | return;\r |
| 450 | }\r |
| 451 | }\r |
| 452 | \r |
| 453 | //_______________________________________________________________________________\r |
| 454 | \r |
| 455 | void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){\r |
| 456 | Float_t mass=0;\r |
| 457 | Int_t abspdg=TMath::Abs(pdg);\r |
| 458 | \r |
| 459 | mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();\r |
| 460 | // compute the Delta mass for the D*\r |
| 461 | if(fCandidateType==kDstartoKpipi){\r |
| 462 | Float_t mass1=0;\r |
| 463 | mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r |
| 464 | mass = mass-mass1;\r |
| 465 | }\r |
| 466 | \r |
| 467 | fMinMass = mass-range;\r |
| 468 | fMaxMass = mass+range;\r |
| 469 | \r |
| 470 | AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));\r |
| 471 | if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");\r |
| 472 | }\r |
| 473 | \r |
| 474 | //_______________________________________________________________________________\r |
| 475 | \r |
| 476 | void AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){\r |
| 477 | if(uplimit>lowlimit)\r |
| 478 | {\r |
| 479 | fMinMass = lowlimit;\r |
| 480 | fMaxMass = uplimit;\r |
| 481 | }\r |
| 482 | else{\r |
| 483 | printf("Error! Lower limit larger than upper limit!\n");\r |
| 484 | fMinMass = uplimit - uplimit*0.2;\r |
| 485 | fMaxMass = uplimit;\r |
| 486 | }\r |
| 487 | }\r |
| 488 | \r |
| 489 | //_______________________________________________________________________________\r |
| 490 | \r |
| 491 | Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){\r |
| 492 | if(nptbins>30) {\r |
| 493 | AliInfo("Maximum number of bins allowed is 30!");\r |
| 494 | return kFALSE;\r |
| 495 | }\r |
| 496 | if(!width) return kFALSE;\r |
| 497 | for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];\r |
| 498 | return kTRUE;\r |
| 499 | }\r |
| 500 | \r |
| 501 | //_______________________________________________________________________________\r |
| 502 | \r |
| 503 | Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{\r |
| 504 | if(!part || !jet){\r |
| 505 | printf("Particle or jet do not exist!\n");\r |
| 506 | return -99;\r |
| 507 | }\r |
| 508 | Double_t p[3],pj[3];\r |
| 509 | Bool_t okpp=part->PxPyPz(p);\r |
| 510 | Bool_t okpj=jet->PxPyPz(pj);\r |
| 511 | if(!okpp || !okpj){\r |
| 512 | printf("Problems getting momenta\n");\r |
| 513 | return -999;\r |
| 514 | }\r |
| 515 | Double_t pjet=jet->P();\r |
| 516 | Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);\r |
| 517 | return z;\r |
| 518 | }\r |
| 519 | \r |
| 520 | //_______________________________________________________________________________\r |
| 521 | \r |
| 522 | Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){\r |
| 523 | \r |
| 524 | // Statistics \r |
| 525 | TH1I* hstat=new TH1I("hstat","Statistics",6,-0.5,5.5);\r |
| 526 | hstat->GetXaxis()->SetBinLabel(1,"N ev anal");\r |
| 527 | hstat->GetXaxis()->SetBinLabel(2,"N ev sel");\r |
| 528 | hstat->GetXaxis()->SetBinLabel(3,"N cand sel cuts");\r |
| 529 | hstat->GetXaxis()->SetBinLabel(4,"N jets");\r |
| 530 | hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");\r |
| 531 | hstat->GetXaxis()->SetBinLabel(6,"N jet rej");\r |
| 532 | hstat->SetNdivisions(1);\r |
| 533 | fmyOutput->Add(hstat);\r |
| 534 | \r |
| 535 | const Int_t nbinsmass=200;\r |
| b770538f |
536 | const Int_t nbinsptjet=500;\r |
| 537 | const Int_t nbinsptD=100;\r |
| 538 | const Int_t nbinsz=100;\r |
| 539 | const Int_t nbinsphi=200;\r |
| 92d9c9cb |
540 | \r |
| b770538f |
541 | const Float_t ptjetlims[2]={0.,200.};\r |
| 542 | const Float_t ptDlims[2]={0.,50.};\r |
| 543 | const Float_t zlims[2]={0.,1.2};\r |
| 544 | const Float_t philims[2]={0.,6.3};\r |
| 92d9c9cb |
545 | \r |
| 546 | if(fCandidateType==kDstartoKpipi) \r |
| 547 | {\r |
| 548 | \r |
| b770538f |
549 | TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);\r |
| 92d9c9cb |
550 | hDiffSideBand->SetStats(kTRUE);\r |
| 551 | hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");\r |
| 552 | hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r |
| b770538f |
553 | hDiffSideBand->Sumw2();\r |
| 92d9c9cb |
554 | fmyOutput->Add(hDiffSideBand); \r |
| 555 | \r |
| 556 | \r |
| 557 | TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);\r |
| 558 | hPtPion->SetStats(kTRUE);\r |
| 559 | hPtPion->GetXaxis()->SetTitle("GeV/c");\r |
| 560 | hPtPion->GetYaxis()->SetTitle("Entries");\r |
| b770538f |
561 | hPtPion->Sumw2();\r |
| 92d9c9cb |
562 | fmyOutput->Add(hPtPion);\r |
| 563 | \r |
| 564 | }\r |
| 565 | \r |
| 566 | // jet related fistograms\r |
| 567 | \r |
| 568 | TH1F* hEjetTrks = new TH1F("hEjetTrks", "Jet tracks energy distribution;Energy (GeV)",500,0,200);\r |
| b770538f |
569 | hEjetTrks->Sumw2();\r |
| 570 | TH1F* hPhiJetTrks = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);\r |
| 571 | hPhiJetTrks->Sumw2();\r |
| 92d9c9cb |
572 | TH1F* hEtaJetTrks = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta", 100,-1.5,1.5);\r |
| b770538f |
573 | hEtaJetTrks->Sumw2();\r |
| 574 | TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);\r |
| 575 | hPtJetTrks->Sumw2();\r |
| 92d9c9cb |
576 | \r |
| 577 | TH1F* hEjet = new TH1F("hEjet", "Jet energy distribution;Energy (GeV)",500,0,200);\r |
| b770538f |
578 | hEjet->Sumw2();\r |
| 579 | TH1F* hPhiJet = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)", nbinsphi,philims[0],philims[1]);\r |
| 580 | hPhiJet->Sumw2();\r |
| 92d9c9cb |
581 | TH1F* hEtaJet = new TH1F("hEtaJet","Jet #eta distribution; #eta", 100,-1.5,1.5);\r |
| b770538f |
582 | hEtaJet->Sumw2();\r |
| 583 | TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);\r |
| 584 | hPtJet->Sumw2();\r |
| 92d9c9cb |
585 | \r |
| b770538f |
586 | TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);\r |
| 587 | hPtJetWithD->Sumw2();\r |
| 92d9c9cb |
588 | //for the MC this histogram is filled with the real background\r |
| b770538f |
589 | TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);\r |
| 590 | hPtJetWithDsb->Sumw2();\r |
| 92d9c9cb |
591 | \r |
| 592 | TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);\r |
| b770538f |
593 | hdeltaRJetTracks->Sumw2();\r |
| 92d9c9cb |
594 | \r |
| 595 | TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);\r |
| b770538f |
596 | hNtrArr->Sumw2();\r |
| 92d9c9cb |
597 | TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);\r |
| b770538f |
598 | hNJetPerEv->Sumw2();\r |
| 92d9c9cb |
599 | \r |
| 600 | fmyOutput->Add(hEjetTrks);\r |
| 601 | fmyOutput->Add(hPhiJetTrks);\r |
| 602 | fmyOutput->Add(hEtaJetTrks);\r |
| 603 | fmyOutput->Add(hPtJetTrks);\r |
| 604 | fmyOutput->Add(hEjet);\r |
| 605 | fmyOutput->Add(hPhiJet);\r |
| 606 | fmyOutput->Add(hEtaJet);\r |
| 607 | fmyOutput->Add(hPtJet);\r |
| 608 | fmyOutput->Add(hPtJetWithD);\r |
| 609 | fmyOutput->Add(hPtJetWithDsb);\r |
| 610 | fmyOutput->Add(hdeltaRJetTracks);\r |
| 611 | fmyOutput->Add(hNtrArr);\r |
| 612 | fmyOutput->Add(hNJetPerEv);\r |
| 613 | \r |
| 614 | TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);\r |
| b770538f |
615 | hDeltaRD->Sumw2();\r |
| 92d9c9cb |
616 | fmyOutput->Add(hDeltaRD);\r |
| 617 | \r |
| 618 | //background (side bands for the Dstar and like sign for D0)\r |
| b770538f |
619 | TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);\r |
| 92d9c9cb |
620 | hInvMassptD->SetStats(kTRUE);\r |
| 621 | hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");\r |
| 622 | hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r |
| b770538f |
623 | hInvMassptD->Sumw2();\r |
| 92d9c9cb |
624 | \r |
| 625 | fmyOutput->Add(hInvMassptD);\r |
| 626 | \r |
| 627 | if(fUseMCInfo){\r |
| b770538f |
628 | TH2F* hInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);\r |
| 92d9c9cb |
629 | hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));\r |
| 630 | hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r |
| b770538f |
631 | hInvMassptDbg->Sumw2();\r |
| 92d9c9cb |
632 | fmyOutput->Add(hInvMassptDbg);\r |
| 633 | \r |
| 634 | }\r |
| b770538f |
635 | const Int_t nAxis=5;\r |
| 636 | const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass};\r |
| 637 | const Double_t minSparse[nAxis]={zlims[0],philims[0],ptjetlims[0],ptDlims[0],fMinMass};\r |
| 638 | const Double_t maxSparse[nAxis]={zlims[1],philims[1],ptjetlims[1],ptDlims[1],fMinMass};\r |
| 639 | THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, and mass", nAxis, nbinsSparse, minSparse, maxSparse);\r |
| 640 | hsDphiz->Sumw2();\r |
| 641 | \r |
| 642 | fmyOutput->Add(hsDphiz);\r |
| 643 | \r |
| 92d9c9cb |
644 | PostData(1, fmyOutput);\r |
| 645 | \r |
| 646 | return kTRUE; \r |
| 647 | }\r |
| 648 | \r |
| 649 | //_______________________________________________________________________________\r |
| 650 | \r |
| 651 | void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){\r |
| 652 | \r |
| 653 | Double_t ptD=candidate->Pt();\r |
| 654 | Double_t ptjet=jet->Pt();\r |
| 655 | Double_t deltaR=DeltaR(candidate,jet);\r |
| 656 | Double_t deltaphi = jet->Phi()-candidate->Phi();\r |
| 657 | if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());\r |
| 658 | if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());\r |
| 659 | Double_t z=Z(candidate,jet);\r |
| 660 | \r |
| 661 | TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");\r |
| 662 | hDeltaRD->Fill(deltaR); \r |
| 663 | if(fUseReco){\r |
| 664 | if(fCandidateType==kD0toKpi) {\r |
| 665 | AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;\r |
| 666 | FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());\r |
| 667 | \r |
| 668 | }\r |
| 669 | \r |
| 670 | if(fCandidateType==kDstartoKpipi) {\r |
| 671 | AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;\r |
| 672 | FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);\r |
| 673 | \r |
| 674 | }\r |
| 675 | } else{ //generated level\r |
| 676 | //AliInfo("Non reco");\r |
| b770538f |
677 | FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);\r |
| 92d9c9cb |
678 | }\r |
| 679 | }\r |
| 680 | \r |
| 681 | //_______________________________________________________________________________\r |
| 682 | \r |
| 683 | void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){\r |
| 684 | \r |
| 685 | //dPhi and z not used at the moment,but will be (re)added\r |
| 686 | \r |
| 687 | Double_t masses[2]={0.,0.};\r |
| 688 | Int_t pdgdaughtersD0[2]={211,321};//pi,K \r |
| 689 | Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r |
| 690 | \r |
| 691 | masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r |
| 692 | masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r |
| 693 | \r |
| 694 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r |
| b770538f |
695 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r |
| 696 | Double_t point[5]={z,dPhi,ptj,ptD,masses[0]};\r |
| 92d9c9cb |
697 | \r |
| 698 | Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);\r |
| 699 | if(isselected==1 || isselected==3) {\r |
| 700 | \r |
| 701 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);\r |
| 702 | \r |
| 703 | FillMassHistograms(masses[0], ptD, deltaR);\r |
| b770538f |
704 | hsDphiz->Fill(point,1.);\r |
| 92d9c9cb |
705 | }\r |
| 706 | if(isselected>=2) {\r |
| 707 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);\r |
| 708 | \r |
| 709 | FillMassHistograms(masses[1], ptD, deltaR);\r |
| b770538f |
710 | point[4]=masses[1];\r |
| 711 | hsDphiz->Fill(point,1.);\r |
| 92d9c9cb |
712 | }\r |
| 713 | \r |
| 714 | }\r |
| 715 | \r |
| 716 | //_______________________________________________________________________________\r |
| 717 | \r |
| 718 | void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){\r |
| 719 | //dPhi and z not used at the moment,but will be (re)added\r |
| 720 | \r |
| 721 | AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();\r |
| 722 | Double_t deltamass= dstar->DeltaInvMass(); \r |
| 723 | //Double_t massD0= dstar->InvMassD0();\r |
| 724 | \r |
| 725 | \r |
| 726 | TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");\r |
| 727 | hPtPion->Fill(softpi->Pt());\r |
| 728 | \r |
| 729 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r |
| b770538f |
730 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r |
| 731 | Double_t point[5]={z,dPhi,ptj,ptD,deltamass};\r |
| 732 | \r |
| 92d9c9cb |
733 | if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);\r |
| 734 | \r |
| 735 | FillMassHistograms(deltamass, ptD, deltaR);\r |
| b770538f |
736 | hsDphiz->Fill(point,1.);\r |
| 92d9c9cb |
737 | \r |
| 738 | }\r |
| 739 | \r |
| 740 | //_______________________________________________________________________________\r |
| 741 | \r |
| b770538f |
742 | void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){\r |
| 92d9c9cb |
743 | \r |
| 744 | Double_t pdgmass=0;\r |
| 745 | TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r |
| b770538f |
746 | THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r |
| 747 | Double_t point[5]={z,dPhi,ptjet,ptD,pdgmass};\r |
| 748 | \r |
| 92d9c9cb |
749 | if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r |
| 750 | if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();\r |
| b770538f |
751 | point[4]=pdgmass;\r |
| 752 | \r |
| 753 | if(deltaR<fJetRadius) {\r |
| 754 | hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet\r |
| 755 | hsDphiz->Fill(point,1.);\r |
| 756 | }\r |
| 92d9c9cb |
757 | \r |
| 758 | }\r |
| 759 | \r |
| 760 | //_______________________________________________________________________________\r |
| 761 | \r |
| 762 | void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){\r |
| 763 | \r |
| 764 | if(deltaR<fJetRadius) {\r |
| 765 | TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");\r |
| 766 | hInvMassptD->Fill(mass,ptD);\r |
| 767 | }\r |
| 768 | }\r |
| 769 | \r |
| 770 | void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){\r |
| 771 | Double_t deltaR=DeltaR(charm, jet);\r |
| 772 | AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;\r |
| 773 | if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;\r |
| 774 | if (deltaR<fJetRadius) jet->AddFlavourTag(tag);\r |
| 775 | \r |
| 776 | return;\r |
| 777 | \r |
| 778 | }\r |
| 779 | \r |
| 780 | //_______________________________________________________________________________\r |
| 781 | \r |
| 782 | void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){\r |
| 783 | \r |
| 784 | // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas \r |
| 785 | // (expected detector resolution) on the left and right frm the D0 mass. Each band\r |
| 786 | // has a width of ~5 sigmas. Two band needed for opening angle considerations \r |
| 787 | TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");\r |
| 788 | TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r |
| 789 | \r |
| 790 | Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r |
| 791 | \r |
| 792 | Int_t bin = fCuts->PtBin(candDstar->Pt());\r |
| 793 | Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];\r |
| 794 | Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];\r |
| 795 | \r |
| 796 | Double_t invM=candDstar->InvMassD0(), deltaM=candDstar->DeltaInvMass(); \r |
| 797 | //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);\r |
| 798 | Double_t ptD=candDstar->Pt();\r |
| 799 | Double_t ptjet=jet->Pt();\r |
| 800 | Double_t dPhi=jet->Phi()-candDstar->Phi();\r |
| 801 | Double_t deltaR=DeltaR(candDstar,jet);\r |
| 802 | if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();\r |
| 803 | if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();\r |
| 804 | \r |
| 805 | //fill the background histogram with the side bands or when looking at MC with the real background\r |
| 806 | if((invM>sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)){ \r |
| 807 | hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background \r |
| 808 | //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);\r |
| 809 | \r |
| 810 | if(deltaR<fJetRadius){ // evaluate in the near side \r |
| 811 | //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);\r |
| 812 | hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r |
| 813 | }\r |
| 814 | }\r |
| 815 | }\r |
| 816 | \r |
| 817 | //_______________________________________________________________________________\r |
| 818 | \r |
| 819 | void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){\r |
| 820 | \r |
| 821 | //need mass, deltaR, pt jet, ptD\r |
| 822 | //two cases: D0 and Dstar\r |
| 823 | \r |
| 824 | Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());\r |
| 825 | if(!isselected) return;\r |
| 826 | \r |
| 827 | Double_t ptD=candbg->Pt();\r |
| 828 | Double_t ptjet=jet->Pt();\r |
| 829 | Double_t deltaR=DeltaR(candbg,jet);\r |
| 830 | \r |
| 831 | TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");\r |
| 832 | TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r |
| 833 | \r |
| 834 | \r |
| 835 | if(fCandidateType==kDstartoKpipi){\r |
| 836 | AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;\r |
| 837 | Double_t deltaM=dstarbg->DeltaInvMass();\r |
| 838 | hInvMassptDbg->Fill(deltaM,ptD);\r |
| 839 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r |
| 840 | }\r |
| 841 | \r |
| 842 | if(fCandidateType==kD0toKpi){\r |
| 843 | Double_t masses[2]={0.,0.};\r |
| 844 | Int_t pdgdaughtersD0[2]={211,321};//pi,K \r |
| 845 | Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r |
| 846 | \r |
| 847 | masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r |
| 848 | masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r |
| 849 | \r |
| 850 | \r |
| 851 | if(isselected==1 || isselected==3) {\r |
| 852 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);\r |
| 853 | hInvMassptDbg->Fill(masses[0],ptD);\r |
| 854 | }\r |
| 855 | if(isselected>=2) {\r |
| 856 | if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);\r |
| 857 | hInvMassptDbg->Fill(masses[1],ptD);\r |
| 858 | }\r |
| 859 | \r |
| 860 | \r |
| 861 | }\r |
| 862 | }\r |
| 863 | \r |
| 864 | //_______________________________________________________________________________\r |
| 865 | \r |
| 866 | Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {\r |
| 867 | //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)\r |
| 868 | \r |
| 869 | if(!p1 || !p2) return -1;\r |
| 870 | Double_t phi1=p1->Phi(),eta1=p1->Eta();\r |
| 871 | Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;\r |
| 872 | \r |
| 873 | Double_t dPhi=phi1-phi2;\r |
| 874 | if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());\r |
| 875 | if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());\r |
| 876 | \r |
| 877 | Double_t dEta=eta1-eta2;\r |
| 878 | Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );\r |
| 879 | return deltaR;\r |
| 880 | \r |
| 881 | }\r |
git://git.uio.no / u / mrichter / AliRoot.git / blame