]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx
CommitLineData
c683985a 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// Task for selecting D mesons to be used as an input for D-Jet correlations
18//
19//-----------------------------------------------------------------------
20// Authors:
21// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22// A.Grelli (Utrecht University) a.grelli@uu.nl
23// Xiaoming Zhang (LBNL) XMZhang@lbl.gov
24//-----------------------------------------------------------------------
25
26#include <TDatabasePDG.h>
27#include <TParticle.h>
28#include <TVector3.h>
29#include "TROOT.h"
30#include <TH3F.h>
31
32#include "AliRDHFCutsDStartoKpipi.h"
33#include "AliRDHFCutsD0toKpi.h"
34#include "AliAODMCHeader.h"
35#include "AliAODHandler.h"
36#include "AliAnalysisManager.h"
37#include "AliLog.h"
38#include "AliAODVertex.h"
39#include "AliAODRecoDecay.h"
40#include "AliAODRecoCascadeHF.h"
41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliESDtrack.h"
43#include "AliAODMCParticle.h"
44#include "AliAnalysisTaskSEDmesonsFilterCJ.h"
45
46ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
47
48//_______________________________________________________________________________
49
50AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
51AliAnalysisTaskSE(),
52fUseMCInfo(kFALSE),
53fUseReco(kTRUE),
54fCandidateType(0),
55fCandidateName(""),
56fPDGmother(0),
57fNProngs(0),
58fBranchName(""),
59fOutput(0),
60fCuts(0),
61fMinMass(0.),
62fMaxMass(0.),
bbb94467 63fCandidateArray(0),
e278882f 64fSideBandArray(0),
d69d07c3 65fhImpPar(),
66fhImpParB(),
67fhInvMassS(),
68fhInvMassB()
c683985a 69{
70 //
71 // Default constructor
72 //
73
74 for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
75 for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
76}
77
78//_______________________________________________________________________________
79
80AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
81AliAnalysisTaskSE(name),
82fUseMCInfo(kFALSE),
83fUseReco(kTRUE),
84fCandidateType(candtype),
85fCandidateName(""),
86fPDGmother(0),
87fNProngs(0),
88fBranchName(""),
89fOutput(0),
90fCuts(cuts),
91fMinMass(0.),
92fMaxMass(0.),
bbb94467 93fCandidateArray(0),
e278882f 94fSideBandArray(0),
d69d07c3 95fhImpPar(),
96fhImpParB(),
97fhInvMassS(),
98fhInvMassB()
c683985a 99{
100 //
101 // Constructor. Initialization of Inputs and Outputs
102 //
103
104 Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
105
106 for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
107 for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
108
109 const Int_t nptbins = fCuts->GetNPtBins();
110 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 };
111
112 switch (fCandidateType) {
113 case 0 :
114 fCandidateName = "D0";
115 fPDGmother = 421;
116 fNProngs = 2;
117 fPDGdaughters[0] = 211; // pi
118 fPDGdaughters[1] = 321; // K
119 fPDGdaughters[2] = 0; // empty
120 fPDGdaughters[3] = 0; // empty
121 fBranchName = "D0toKpi";
122 break;
123 case 1 :
124 fCandidateName = "Dstar";
125 fPDGmother = 413;
126 fNProngs = 3;
127 fPDGdaughters[1] = 211; // pi soft
128 fPDGdaughters[0] = 421; // D0
129 fPDGdaughters[2] = 211; // pi fromD0
130 fPDGdaughters[3] = 321; // K from D0
131 fBranchName = "Dstar";
132
133 if (nptbins<=13) {
134 for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
135 } else {
136 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
137 }
138 break;
139 default :
140 printf("%d not accepted!!\n",fCandidateType);
141 break;
142 }
143
144 if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
145 if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
146
147 DefineOutput(1, TList::Class()); // histos
148 DefineOutput(2, AliRDHFCuts::Class()); // my cuts
bbb94467 149 DefineOutput(3, TClonesArray::Class()); //array of candidates
150 DefineOutput(4, TClonesArray::Class()); //array of SB candidates
c683985a 151}
152
153//_______________________________________________________________________________
154
155AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
156{
157 //
158 // destructor
159 //
160
161 Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
162
163 if (fOutput) { delete fOutput; fOutput = 0; }
164 if (fCuts) { delete fCuts; fCuts = 0; }
165 if (fCandidateArray) { delete fCandidateArray; fCandidateArray = 0; }
bbb94467 166 delete fSideBandArray;
c683985a 167
168}
169
170//_______________________________________________________________________________
171
172void AliAnalysisTaskSEDmesonsFilterCJ::Init()
173{
174 //
175 // Initialization
176 //
177
178 if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
179
180 switch (fCandidateType) {
bbb94467 181 case 0:
182 {
c683985a 183 AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
184 copyfCutsDzero->SetName("AnalysisCutsDzero");
185 PostData(2, copyfCutsDzero); // Post the data
bbb94467 186 } break;
187 case 1:
188 {
189 AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
190 copyfCutsDstar->SetName("AnalysisCutsDStar");
191 PostData(2, copyfCutsDstar); // Post the cuts
192 } break;
193 default: return;
c683985a 194 }
195
196 return;
197}
198
199//_______________________________________________________________________________
200
201void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
202{
203 //
204 // output
205 //
206
207 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
208
209 fOutput = new TList(); fOutput->SetOwner();
210 DefineHistoForAnalysis(); // define histograms
211
bbb94467 212 if (fCandidateType==kD0toKpi){
213 fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
214 fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0);
215 }
c683985a 216
bbb94467 217 if (fCandidateType==kDstartoKpipi) {
218 fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
219 fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0);
c683985a 220 }
221
bbb94467 222 fCandidateArray->SetOwner();
223 fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
224
225 //this is used for the DStar side bands and MC!
226 fSideBandArray->SetOwner();
227 fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
228
c683985a 229 PostData(1, fOutput);
bbb94467 230 PostData(3, fCandidateArray);
231 PostData(4, fSideBandArray);
232
c683985a 233 return;
234}
235
236//_______________________________________________________________________________
237
238void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
239 //
240 // user exec
241 //
c683985a 242 // Load the event
243 AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
244
245 TClonesArray *arrayDStartoD0pi = 0;
246 if (!aodEvent && AODEvent() && IsStandardAOD()) {
247
248 // In case there is an AOD handler writing a standard AOD, use the AOD
249 // event in memory rather than the input (ESD) event.
250 aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
251
252 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
253 // have to taken from the AOD event hold by the AliAODExtension
254 AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
255 if(aodHandler->GetExtensions()) {
256 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
257 AliAODEvent *aodFromExt = ext->GetAOD();
258 arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
259 }
260 } else {
261 arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
262 }
263
264 if (!arrayDStartoD0pi) {
265 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
266 return;
267 } else {
268 AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
269 }
270
271 TClonesArray* mcArray = 0x0;
272 if (fUseMCInfo) {
273 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
274 if (!mcArray) {
275 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
276 return;
277 }
278 }
279
280 //Histograms
281 TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
282 TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
283 TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
284 TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
285
286 TH1F* hPtPion=0x0;
287 if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
288 hstat->Fill(0);
289
290 // fix for temporary bug in ESDfilter
291 // the AODs with null vertex pointer didn't pass the PhysSel
292 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
293
294 //Event selection
295 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
296 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
297 if (!iseventselected) return;
298 hstat->Fill(1);
299
300 const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
e278882f 301 if(!fUseMCInfo) hstat->Fill(2, nD);
c683985a 302
303 //D* and D0 prongs needed to MatchToMC method
304 Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
305 Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
306
307 //D0 from D0 bar
308 Int_t pdgdaughtersD0[2] = { 211, 321 }; // pi,K
309 Int_t pdgdaughtersD0bar[2] = { 321, 211 }; // K,pi
310
311 Int_t iCand =0;
312 Int_t iSBCand=0;
313 Int_t isSelected = 0;
314 AliAODRecoDecayHF *charmCand = 0;
bbb94467 315 AliAODRecoCascadeHF *dstar = 0;
c683985a 316 AliAODMCParticle *charmPart = 0;
317 Bool_t isMCBkg=kFALSE;
318
319 Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
320
321 Int_t mcLabel = -9999;
322 Int_t pdgMeson = 413;
323 if (fCandidateType==kD0toKpi) pdgMeson = 421;
324
bbb94467 325 //clear the TClonesArray from the previous event
326 fCandidateArray->Clear();
327 fSideBandArray->Clear();
328
c683985a 329 for (Int_t icharm=0; icharm<nD; icharm++) { //loop over D candidates
330 charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
331 if (!charmCand) continue;
332
8577f7bc 333 TString smcTruth="S";
334
bbb94467 335 if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
c683985a 336
337 if (fUseMCInfo) { // Look in MC, try to simulate the z
338 if (fCandidateType==kDstartoKpipi) {
bbb94467 339 mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
c683985a 340 }
341
342 if (fCandidateType==kD0toKpi)
343 mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
344
e278882f 345 if (mcLabel<=0) {
346 isMCBkg=kTRUE;
347 hstat->Fill(5);
348 }
349 else {
350 hstat->Fill(2);
351
352 }
c683985a 353 if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
8577f7bc 354
355 if (isMCBkg) smcTruth="B";
356
c683985a 357 }
358
359 Double_t ptD = charmCand->Pt();
360
361 // region of interest + cuts
362 if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
363
364 if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
365 //select by track cuts the side band candidates (don't want mass cut)
366 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
367 if (!isSelected) continue;
368 //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
369 Int_t bin = fCuts->PtBin(ptD);
370 if(bin<0 || bin>=fCuts->GetNPtBins()) {
371 AliError(Form("Pt %.3f out of bounds",ptD));
372 continue;
373 }
c683985a 374 //if data and Dstar from D0 side band
bbb94467 375 if (((dstar->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/|| ((dstar->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
c683985a 376
bbb94467 377 new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
c683985a 378 iSBCand++;
379 }
380 }
381 //candidate selected by cuts and PID
382 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
383 if (!isSelected) continue;
384
e278882f 385 Int_t nprongs=charmCand->GetNProngs();
386
c683985a 387 //for data and MC signal fill fCandidateArray
388 if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
389 // for data or MC with the requirement fUseReco fill with candidates
e278882f 390
c683985a 391 if(fUseReco) {
bbb94467 392 if (fCandidateType==kDstartoKpipi){
393 new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
d69d07c3 394 fhInvMassS->Fill(dstar->DeltaInvMass(),dstar->Pt());
e278882f 395
bbb94467 396 } else{
397 new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
398 //Printf("Filling reco");
e278882f 399 for(Int_t kd=0;kd<nprongs;kd++){
400 fhImpPar->Fill(charmCand->Getd0Prong(kd),charmCand->PtProng(kd));
401 }
402 Double_t masses[2];
403 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
404 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
405 fhInvMassS->Fill(masses[0]);
406 fhInvMassS->Fill(masses[1]);
407
bbb94467 408 }
409
c683985a 410 hstat->Fill(3);
e278882f 411
c683985a 412 }
413 // for MC with requirement particle level fill with AliAODMCParticle
414 else if (fUseMCInfo) {
415 new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
416 //Printf("Filling gen");
417 hstat->Fill(3);
418 }
419
420 iCand++;
421 }
422 //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
423 else if(fUseReco){
bbb94467 424 if (fCandidateType==kDstartoKpipi){
425 new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
e278882f 426 fhInvMassB->Fill(dstar->DeltaInvMass());
bbb94467 427 }
428 if (fCandidateType==kD0toKpi){
429 new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
e278882f 430 Double_t masses[2];
431 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
432 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
433 fhInvMassB->Fill(masses[0]);
434 fhInvMassB->Fill(masses[1]);
435
bbb94467 436 }
c683985a 437 iSBCand++;
438 }
439
440
441 Double_t masses[2];
442 if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
443
444 //softpion from D* decay
bbb94467 445 AliAODTrack *track2 = (AliAODTrack*)dstar->GetBachelor();
c683985a 446
447 // select D* in the D0 window.
448 // In the cut object window is loose to allow for side bands
449
450
451 // retrieve the corresponding pt bin for the candidate
452 // and set the expected D0 width (x3)
453 // static const Int_t n = fCuts->GetNPtBins();
454 Int_t bin = fCuts->PtBin(ptD);
455 if(bin<0 || bin>=fCuts->GetNPtBins()) {
456 AliError(Form("Pt %.3f out of bounds",ptD));
457 continue;
458 }
459
460 AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
461 //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
bbb94467 462 if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
463 masses[0] = dstar->DeltaInvMass(); //D*
c683985a 464 masses[1] = 0.; //dummy for D*
465
466 //D* delta mass
467 hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
468
469 // D* pt and soft pion pt for good candidates
470 Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
bbb94467 471 Double_t invmassDelta = dstar->DeltaInvMass();
c683985a 472 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
473 }
8577f7bc 474
475 if (fUseMCInfo){ //fill histograms of kinematics, using MC truth
476 //get histos
477 TH2F *halphaDD = (TH2F*)fOutput->FindObject(Form("halphaDD%s",smcTruth.Data()));
478 TH2F *halphaDpis = (TH2F*)fOutput->FindObject(Form("halphaDpis%s",smcTruth.Data()));
479 TH2F *halphaDpi = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
480 TH2F *halphaDK = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
481
482 TH2F *hdeltaRDD = (TH2F*)fOutput->FindObject(Form("hdeltaRDD%s",smcTruth.Data()));
483 TH2F *hdeltaRDpis = (TH2F*)fOutput->FindObject(Form("hdeltaRDpis%s",smcTruth.Data()));
484 TH2F *hdeltaRDpi = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
485 TH2F *hdeltaRDK = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
486
487 Double_t aD = dstar->Phi(),
488 apis= track2->Phi();
489
490 AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
491 Double_t aD0 = D0fromDstar->Phi();
492 Int_t isD0= D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
493 Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1),
494 api= isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
495 Double_t dRDD0 = DeltaR(dstar,D0fromDstar),
496 dRDpis = DeltaR(dstar,track2),
497 dRDpi = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(1) : (AliVParticle*)D0fromDstar->GetDaughter(0)),
498 dRDK = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(0) : (AliVParticle*)D0fromDstar->GetDaughter(1));
499
500 halphaDD-> Fill(aD-aD0,ptD);
501 halphaDpis->Fill(aD-apis,ptD);
502 halphaDpi-> Fill(aD-api,ptD);
503 halphaDK-> Fill(aD-aK,ptD);
504
505 hdeltaRDD-> Fill(dRDD0,ptD);
506 hdeltaRDpis->Fill(dRDpis,ptD);
507 hdeltaRDpi-> Fill(dRDpi,ptD);
508 hdeltaRDK-> Fill(dRDK,ptD);
509
e278882f 510 if (isMCBkg){
511 fhImpParB->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
512 fhImpParB->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
513 fhImpParB->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
514
515 } else {
516 fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
517 fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
518 fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
519
520 }
521
522
523 } else{
524 fhImpPar->Fill(charmCand->Getd0Prong(0),charmCand->PtProng(0)); //bachelor
525 AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
526 fhImpPar->Fill(D0fromDstar->Getd0Prong(0),D0fromDstar->PtProng(0));
527 fhImpPar->Fill(D0fromDstar->Getd0Prong(1),D0fromDstar->PtProng(1));
8577f7bc 528 }
529
c683985a 530 } //Dstar specific
531
532 if (fCandidateType==kD0toKpi) { //D0->Kpi
533
534 //needed quantities
535 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
536 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
537 hstat->Fill(3);
538
539 // mass vs pt
540 if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
541 if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
8577f7bc 542
543 if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
544
545 Double_t aD = charmCand->Phi();
546 Double_t adaugh[2]={charmCand->PhiProng(0),charmCand->PhiProng(1)};
547 AliAODTrack* p0=(AliAODTrack*)charmCand->GetDaughter(0);
548 AliAODTrack* p1=(AliAODTrack*)charmCand->GetDaughter(1);
549 Float_t dR0 = DeltaR(charmCand, p0), dR1 = DeltaR(charmCand, p1);
550 Bool_t isD0=kFALSE;
551 if(mcLabel==421) isD0=kTRUE;
552 if(mcLabel==-421) isD0=kFALSE;
553
554 if(isMCBkg) { //background
555 TH2F *halphaDpi = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
556 TH2F *halphaDK = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
557
558 TH2F *hdeltaRDpi = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
559 TH2F *hdeltaRDK = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
560
561
562 if (isSelected==1 || isSelected==3) { // selected as D0
563 halphaDK->Fill(aD-adaugh[0],ptD);
564 halphaDpi->Fill(aD-adaugh[1],ptD);
565
566 hdeltaRDK->Fill(dR0,ptD);
567 hdeltaRDpi->Fill(dR1,ptD);
568
569 }
570 if (isSelected>=2) { //selected as D0bar
571 halphaDpi->Fill(aD-adaugh[0],ptD);
572 halphaDK->Fill(aD-adaugh[1],ptD);
573
574 hdeltaRDpi->Fill(dR0,ptD);
575 hdeltaRDK->Fill(dR1,ptD);
576
577 }
578
579 }else{ //signal and reflections
580 TH2F *halphaDpiS = (TH2F*)fOutput->FindObject("halphaDpiS");
581 TH2F *halphaDKS = (TH2F*)fOutput->FindObject("halphaDKS");
582 TH2F *halphaDpiR = (TH2F*)fOutput->FindObject("halphaDpiR");
583 TH2F *halphaDKR = (TH2F*)fOutput->FindObject("halphaDKR");
584
585 TH2F *hdeltaRDpiS = (TH2F*)fOutput->FindObject("hdeltaRDpiS");
586 TH2F *hdeltaRDKS = (TH2F*)fOutput->FindObject("hdeltaRDKS");
587 TH2F *hdeltaRDpiR = (TH2F*)fOutput->FindObject("hdeltaRDpiR");
588 TH2F *hdeltaRDKR = (TH2F*)fOutput->FindObject("hdeltaRDKR");
589
590 if(isD0) { //D0
591 halphaDKS->Fill(aD-adaugh[0],ptD);
592 halphaDpiS->Fill(aD-adaugh[1],ptD);
593
594 hdeltaRDKS->Fill(dR0,ptD);
595 hdeltaRDpiS->Fill(dR1,ptD);
596 if(isSelected>=2){ //selected as D0bar
597 halphaDpiR->Fill(aD-adaugh[0],ptD);
598 halphaDKR->Fill(aD-adaugh[1],ptD);
599
600 hdeltaRDpiR->Fill(dR0,ptD);
601 hdeltaRDKR->Fill(dR1,ptD);
602 }
603 } else { //D0bar
604 halphaDKS->Fill(aD-adaugh[1],ptD);
605 halphaDpiS->Fill(aD-adaugh[0],ptD);
606
607 hdeltaRDKS->Fill(dR1,ptD);
608 hdeltaRDpiS->Fill(dR0,ptD);
609
610 if(isSelected>=2){ //selected as D0bar
611 halphaDpiR->Fill(aD-adaugh[1],ptD);
612 halphaDKR->Fill(aD-adaugh[0],ptD);
613
614 hdeltaRDpiR->Fill(dR1,ptD);
615 hdeltaRDKR->Fill(dR0,ptD);
616 }
617 }
618
619 } //end signal and reflections
620
621
622 }// end MC
623
624
c683985a 625 } //D0 specific
626
627 charmCand = 0;
628 isMCBkg=kFALSE;
629 } // end of D cand loop
630
631 hnCandEv->Fill(fCandidateArray->GetEntriesFast());
632 if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
633 Int_t nsbcand=fSideBandArray->GetEntriesFast();
634 hstat->Fill(4,nsbcand);
635 hnSBCandEv->Fill(nsbcand);
636 }
c3c1ab31 637 //Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
bbb94467 638 PostData(1, fOutput);
639 PostData(3, fCandidateArray);
640 PostData(4, fSideBandArray);
641
c683985a 642 return;
643}
644
645//_______________________________________________________________________________
646
647void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
648{
649 // The Terminate() function is the last function to be called during
650 // a query. It always runs on the client, it can be used to present
651 // the results graphically or save the results to file.
652
653 Info("Terminate"," terminate");
654 AliAnalysisTaskSE::Terminate();
655
656 fOutput = dynamic_cast<TList*>(GetOutputData(1));
657 if (!fOutput) {
658 printf("ERROR: fOutput not available\n");
659 return;
660 }
661
662 return;
663}
664
665//_______________________________________________________________________________
666
667void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
668{
669 //
670 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
671 //
672
673 Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
674
675 // compute the Delta mass for the D*
676 if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
677
678
679 fMinMass = mass - range;
680 fMaxMass = mass + range;
681
682 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
683 if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
684
685 return;
686}
687
688//_______________________________________________________________________________
689
690void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
691{
692 //
693 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
694 //
695
696 if (uplimit>lowlimit) {
697 fMinMass = lowlimit;
698 fMaxMass = uplimit;
699 } else {
700 printf("Error! Lower limit larger than upper limit!\n");
701 fMinMass = uplimit - uplimit*0.2;
702 fMaxMass = uplimit;
703 }
704
705 return;
706}
707
708//_______________________________________________________________________________
709
710Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
711{
712 //
713 // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
714 //
715
716 if (nptbins>30) {
717 AliInfo("Maximum number of bins allowed is 30!");
718 return kFALSE;
719 }
720
721 if (!width) return kFALSE;
722 for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
723
724 return kTRUE;
725}
726
727//_______________________________________________________________________________
728
729Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
730{
731 //
732 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
733 //
734
735 // Statistics
e278882f 736 TH1I* hstat = new TH1I("hstat","Statistics",6,-0.5,5.5);
c683985a 737 hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
738 hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
e278882f 739 if(fUseMCInfo){
740 if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N D");
741 else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
742 } else hstat->GetXaxis()->SetBinLabel(3, "N Cand");
c683985a 743 hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
744 if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
e278882f 745 if(fUseMCInfo) {
746 hstat->GetXaxis()->SetBinLabel(6, "N Background");
747 }
c683985a 748 hstat->SetNdivisions(1);
749 fOutput->Add(hstat);
750
751 TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
752 fOutput->Add(hnCandEv);
753
754 // Invariant mass related histograms
755 const Int_t nbinsmass = 200;
8577f7bc 756 const Int_t ptbinsD=100;
757 Float_t ptmin=0.,ptmax=50.;
758 TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
c683985a 759 hInvMass->SetStats(kTRUE);
760 hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
761 hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
762 fOutput->Add(hInvMass);
8577f7bc 763 if ((fCandidateType==kDstartoKpipi) || fUseMCInfo){
c683985a 764 TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
765 fOutput->Add(hnSBCandEv);
8577f7bc 766 }
767 if (fCandidateType==kDstartoKpipi) {
c683985a 768
769 TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
770 hPtPion->SetStats(kTRUE);
771 hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
772 hPtPion->GetYaxis()->SetTitle("entries");
773 fOutput->Add(hPtPion);
774 }
775
e278882f 776 fhImpPar = new TH2F("hImpPar", "Impact parameter of daughter tracks; Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
777 fOutput->Add(fhImpPar);
d69d07c3 778
779 fhInvMassS = new TH1F("hInvMassS", "D invariant mass distribution (filled with fCandidateArray)", nbinsmass, fMinMass, fMaxMass);
780 fhInvMassS->SetStats(kTRUE);
781 fhInvMassS->GetXaxis()->SetTitle("mass (GeV/c)");
782 fhInvMassS->GetYaxis()->SetTitle("p_{T} (GeV/c)");
783 fOutput->Add(fhInvMassS);
784
8577f7bc 785 const Int_t nbinsalpha=200;
786 Float_t minalpha=-TMath::Pi(), maxalpha=TMath::Pi();
787 const Int_t nbinsdeltaR= 200;
788 Float_t mindeltaR = 0., maxdeltaR = 10.;
789 if(fUseMCInfo){
790 if (fCandidateType==kDstartoKpipi){
791 TH2F* halphaDDS =new TH2F("halphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
792 TH2F* halphaDpisS=new TH2F("halphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
793 TH2F* halphaDpiS =new TH2F("halphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
794 TH2F* halphaDKS =new TH2F("halphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
795
796 TH2F* halphaDDB =new TH2F("halphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
797 TH2F* halphaDpisB=new TH2F("halphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
798 TH2F* halphaDpiB =new TH2F("halphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
799 TH2F* halphaDKB =new TH2F("halphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
800
801 TH2F* hdeltaRDDS =new TH2F("hdeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
802 TH2F* hdeltaRDpisS=new TH2F("hdeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
803 TH2F* hdeltaRDpiS =new TH2F("hdeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
804 TH2F* hdeltaRDKS =new TH2F("hdeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
805
806 TH2F* hdeltaRDDB =new TH2F("hdeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
807 TH2F* hdeltaRDpisB=new TH2F("hdeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
808 TH2F* hdeltaRDpiB =new TH2F("hdeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
809 TH2F* hdeltaRDKB =new TH2F("hdeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
810
811 fOutput->Add(halphaDDS);
812 fOutput->Add(halphaDpisS);
813 fOutput->Add(halphaDpiS);
814 fOutput->Add(halphaDKS);
815 fOutput->Add(halphaDDB);
816 fOutput->Add(halphaDpisB);
817 fOutput->Add(halphaDpiB);
818 fOutput->Add(halphaDKB);
819
820 fOutput->Add(hdeltaRDDS);
821 fOutput->Add(hdeltaRDpisS);
822 fOutput->Add(hdeltaRDpiS);
823 fOutput->Add(hdeltaRDKS);
824 fOutput->Add(hdeltaRDDB);
825 fOutput->Add(hdeltaRDpisB);
826 fOutput->Add(hdeltaRDpiB);
827 fOutput->Add(hdeltaRDKB);
828 }
829
830 if (fCandidateType==kD0toKpi){
831
832 TH2F* halphaDpiS=new TH2F("halphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
833 TH2F* halphaDKS =new TH2F("halphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
834 TH2F* halphaDpiR=new TH2F("halphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
835 TH2F* halphaDKR =new TH2F("halphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
836
837 TH2F* halphaDpiB=new TH2F("halphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
838 TH2F* halphaDKB =new TH2F("halphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
839
840
841 TH2F* hdeltaRDpiS=new TH2F("hdeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
842 TH2F* hdeltaRDKS =new TH2F("hdeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
843 TH2F* hdeltaRDpiR=new TH2F("hdeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
844 TH2F* hdeltaRDKR =new TH2F("hdeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
845
846 TH2F* hdeltaRDpiB=new TH2F("hdeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
847 TH2F* hdeltaRDKB =new TH2F("hdeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
848
849 fOutput->Add(halphaDpiS);
850 fOutput->Add(halphaDKS);
851 fOutput->Add(halphaDpiR);
852 fOutput->Add(halphaDKR);
853 fOutput->Add(halphaDpiB);
854 fOutput->Add(halphaDKB);
855
856 fOutput->Add(hdeltaRDpiS);
857 fOutput->Add(hdeltaRDKS);
858 fOutput->Add(hdeltaRDpiR);
859 fOutput->Add(hdeltaRDKR);
860 fOutput->Add(hdeltaRDpiB);
861 fOutput->Add(hdeltaRDKB);
862
863 }
e278882f 864 fhImpParB = new TH2F("hImpParB", "Impact parameter of daughter tracks (Background); Getd0Prong();#it{p}^{daugh}_{T} (GeV/c)",200, -0.1,0.1,ptbinsD, ptmin, ptmax); //same range of pt of the D, but pt daughter used
865 fOutput->Add(fhImpParB);
d69d07c3 866
e278882f 867 fhInvMassB = new TH1F("hInvMassB", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass);
868 fhInvMassB->SetStats(kTRUE);
869 fhInvMassB->GetXaxis()->SetTitle("mass (GeV/c)");
870 fhInvMassB->GetYaxis()->SetTitle("p_{T} (GeV/c)");
871 fOutput->Add(fhInvMassB);
872
8577f7bc 873 }
c683985a 874 return kTRUE;
875}
8577f7bc 876
877//_______________________________________________________________________________
878
879Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
880 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
881
882 if(!p1 || !p2) return -1;
883 Double_t phi1=p1->Phi(),eta1=p1->Eta();
884 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
885
886 Double_t dPhi=phi1-phi2;
887 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
888 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
889
890 Double_t dEta=eta1-eta2;
891 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
892 return deltaR;
893
894}