]>
Commit | Line | Data |
---|---|---|
10eaad41 | 1 | /************************************************************************** |
2 | * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr) * | |
3 | * Contributors are mentioned in the code where appropriate. * | |
4 | * * | |
5 | * Permission to use, copy, modify and distribute this software and its * | |
6 | * documentation strictly for non-commercial purposes is hereby granted * | |
7 | * without fee, provided that the above copyright notice appears in all * | |
8 | * copies and that both the copyright notice and this permission notice * | |
9 | * appear in the supporting documentation. The authors make no claims * | |
10 | * about the suitability of this software for any purpose. It is * | |
11 | * provided "as is" without express or implied warranty. * | |
12 | **************************************************************************/ | |
13 | ||
14 | //---------------------------------------------------------------------------------------------------------------- | |
f9afc48d | 15 | // class AliResonanceKink |
10eaad41 | 16 | // Example of an analysis task for reconstructing resonances having at least one kaon-kink in their decay |
17 | // products. It provides basic plots as well as plots helping to calculate the corrections. | |
18 | // Usage: To analyse a resonance having a kaon in its decay products, one has to modify the integer | |
f9afc48d | 19 | // variables resonancePDG, fdaughter1pdg and fdaughter2pdg accordingly as well as daughter1pdgMass and daughter2pdgMass. |
10eaad41 | 20 | // Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed |
21 | //----------------------------------------------------------------------------------------------------------------- | |
22 | ||
23 | #include "TChain.h" | |
24 | #include "TTree.h" | |
25 | #include "TH2D.h" | |
26 | #include "TParticle.h" | |
92adf4f6 | 27 | #include "TDatabasePDG.h" |
28 | #include "TParticlePDG.h" | |
10eaad41 | 29 | #include "TF1.h" |
f9afc48d | 30 | #include "TList.h" |
31 | #include "TString.h" | |
10eaad41 | 32 | #include "AliMCEventHandler.h" |
33 | #include "AliMCEvent.h" | |
f9afc48d | 34 | #include "AliResonanceKink.h" |
10eaad41 | 35 | #include "AliESDkink.h" |
36 | #include "AliStack.h" | |
f9afc48d | 37 | #include "AliESDtrack.h" |
38 | #include "AliESDEvent.h" | |
92adf4f6 | 39 | #include "AliExternalTrackParam.h" |
7ccf0419 | 40 | #include "AliMCParticle.h" |
10eaad41 | 41 | |
f9afc48d | 42 | ClassImp(AliResonanceKink) |
10eaad41 | 43 | |
44 | //________________________________________________________________________ | |
f9afc48d | 45 | AliResonanceKink::AliResonanceKink() |
a1cdcf46 | 46 | : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0), |
92adf4f6 | 47 | fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(0), fLowX(0), fHighX(0), fdaughter1pdg(0), fdaughter2pdg(0), fresonancePDGcode(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0), |
48 | fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE) | |
7ccf0419 | 49 | , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0) |
10eaad41 | 50 | { |
51 | // Constructor | |
52 | } | |
53 | ||
54 | //________________________________________________________________________ | |
f9afc48d | 55 | AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG) |
a1cdcf46 | 56 | : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0), |
92adf4f6 | 57 | fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), |
58 | fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE) | |
7ccf0419 | 59 | , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0) |
10eaad41 | 60 | { |
f9afc48d | 61 | // Constructor |
10eaad41 | 62 | |
10eaad41 | 63 | fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0); |
64 | ||
f9afc48d | 65 | fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX); |
66 | fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX); | |
67 | fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020) | |
10eaad41 | 68 | |
69 | fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0); | |
70 | fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1); | |
71 | fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
72 | fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0); | |
73 | fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); | |
74 | fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
75 | fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0); | |
76 | fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1); | |
77 | fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
78 | fhdr=new TH1D("fhdr"," ", 100,0.0,5.0); | |
79 | fhdz=new TH1D("fhdz"," ", 100,0.0,5.0); | |
f9afc48d | 80 | |
81 | f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0); | |
82 | f1->SetParameter(0,0.493677); | |
83 | f1->SetParameter(1,0.9127037); | |
84 | f1->SetParameter(2,TMath::Pi()); | |
85 | ||
86 | f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0); | |
87 | f2->SetParameter(0,0.13957018); | |
88 | f2->SetParameter(1,0.2731374); | |
89 | f2->SetParameter(2,TMath::Pi()); | |
90 | ||
10eaad41 | 91 | fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0); |
92adf4f6 | 92 | |
7ccf0419 | 93 | fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0); |
94 | fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
95 | fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
96 | fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
10eaad41 | 97 | } |
98 | ||
99 | //________________________________________________________________________ | |
f9afc48d | 100 | AliResonanceKink:: ~AliResonanceKink() |
10eaad41 | 101 | { |
f9afc48d | 102 | // Destructor |
103 | if(fOpeningAngle) delete fOpeningAngle; | |
104 | if(fInvariantMass) delete fInvariantMass; | |
105 | if(fInvMassTrue) delete fInvMassTrue; | |
106 | if(fPhiBothKinks) delete fPhiBothKinks; | |
107 | if(fRecPt) delete fRecPt; | |
108 | if(fRecEta) delete fRecEta; | |
109 | if(fRecEtaPt) delete fRecEtaPt; | |
110 | if(fSimPt) delete fSimPt; | |
111 | if(fSimEta) delete fSimEta; | |
112 | if(fSimEtaPt) delete fSimEtaPt; | |
113 | if(fSimPtKink) delete fSimPtKink; | |
114 | if(fSimEtaKink) delete fSimEtaKink; | |
115 | if(fSimEtaPtKink) delete fSimEtaPtKink; | |
116 | if(fhdr) delete fhdr; | |
117 | if(fhdz) delete fhdz; | |
7ccf0419 | 118 | if(fvtxz) delete fvtxz; |
119 | if(fInvmassPt) delete fInvmassPt; | |
120 | if(fInvmassPtTrue) delete fInvmassPtTrue; | |
121 | if(fMCInvmassPt) delete fMCInvmassPt; | |
122 | if(fMCInvmassPtTrue) delete fMCInvmassPtTrue; | |
f9afc48d | 123 | } |
124 | //________________________________________________________________________ | |
125 | TList* AliResonanceKink::GetHistogramList() | |
126 | { | |
127 | // Adding histograms to the list | |
128 | fListOfHistos=new TList(); | |
10eaad41 | 129 | |
f9afc48d | 130 | fListOfHistos->Add(fOpeningAngle); |
131 | fListOfHistos->Add(fInvariantMass); | |
132 | fListOfHistos->Add(fInvMassTrue); | |
133 | fListOfHistos->Add(fPhiBothKinks); | |
134 | fListOfHistos->Add(fRecPt); | |
135 | fListOfHistos->Add(fRecEta); | |
136 | fListOfHistos->Add(fRecEtaPt); | |
137 | fListOfHistos->Add(fSimPt); | |
138 | fListOfHistos->Add(fSimEta); | |
139 | fListOfHistos->Add(fSimEtaPt); | |
140 | fListOfHistos->Add(fSimPtKink); | |
141 | fListOfHistos->Add(fSimEtaKink); | |
142 | fListOfHistos->Add(fSimEtaPtKink); | |
143 | fListOfHistos->Add(fhdr); | |
144 | fListOfHistos->Add(fhdz); | |
145 | fListOfHistos->Add(fvtxz); | |
7ccf0419 | 146 | fListOfHistos->Add(fInvmassPt); |
147 | fListOfHistos->Add(fInvmassPtTrue); | |
148 | fListOfHistos->Add(fMCInvmassPt); | |
149 | fListOfHistos->Add(fMCInvmassPtTrue); | |
150 | ||
f9afc48d | 151 | return fListOfHistos; |
152 | } | |
153 | ||
154 | //________________________________________________________________________ | |
155 | void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx) | |
156 | { | |
157 | // Initialisation of the output histograms | |
158 | fNbins=nbins; | |
159 | fLowX=nlowx; | |
160 | fHighX=nhighx; | |
161 | ||
162 | fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0); | |
163 | ||
164 | fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX); | |
165 | fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX); | |
166 | fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020) | |
167 | ||
168 | fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0); | |
169 | fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1); | |
170 | fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
171 | fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0); | |
172 | fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); | |
173 | fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
174 | fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0); | |
175 | fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1); | |
176 | fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1); | |
177 | fhdr=new TH1D("fhdr"," ", 100,0.0,5.0); | |
178 | fhdz=new TH1D("fhdz"," ", 100,0.0,5.0); | |
10eaad41 | 179 | |
f9afc48d | 180 | f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0); |
181 | f1->SetParameter(0,0.493677); | |
182 | f1->SetParameter(1,0.9127037); | |
183 | f1->SetParameter(2,TMath::Pi()); | |
184 | ||
185 | f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0); | |
186 | f2->SetParameter(0,0.13957018); | |
187 | f2->SetParameter(1,0.2731374); | |
188 | f2->SetParameter(2,TMath::Pi()); | |
f27d6e67 | 189 | |
f9afc48d | 190 | fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0); |
7ccf0419 | 191 | |
192 | fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
193 | fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
194 | fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
195 | fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0); | |
196 | ||
f9afc48d | 197 | } |
198 | ||
199 | //________________________________________________________________________ | |
200 | void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) | |
201 | { | |
202 | // Main loop | |
203 | // Called for each event | |
204 | Int_t resonancePDGcode, antiresonancePDGcode; | |
205 | ||
206 | if (fdaughter1pdg==kdaughterKaon) { | |
207 | resonancePDGcode=fresonancePDGcode; | |
208 | antiresonancePDGcode=-fresonancePDGcode; | |
10eaad41 | 209 | } |
f9afc48d | 210 | if (fdaughter1pdg!=kdaughterKaon) { |
211 | resonancePDGcode=-fresonancePDGcode; | |
212 | antiresonancePDGcode=fresonancePDGcode; | |
213 | } | |
214 | if (fdaughter1pdg==fdaughter2pdg) { | |
215 | resonancePDGcode=fresonancePDGcode; | |
216 | antiresonancePDGcode=fresonancePDGcode; | |
217 | } | |
10eaad41 | 218 | |
f9afc48d | 219 | Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass(); |
220 | Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass(); | |
221 | ||
222 | if (!esd) { | |
10eaad41 | 223 | Printf("ERROR: fESD not available"); |
224 | return; | |
225 | } | |
226 | ||
f9afc48d | 227 | if (!mcEvent) { |
228 | Printf("ERROR: mcEvent not available"); | |
229 | return; | |
230 | } | |
231 | ||
10eaad41 | 232 | AliStack* stack=mcEvent->Stack(); |
233 | ||
234 | if(fAnalysisType == "MC") { | |
7ccf0419 | 235 | |
236 | const AliESDVertex* vertex = GetEventVertex(esd); | |
237 | if(!vertex) return; | |
238 | Double_t vtx[3]; | |
239 | vertex->GetXYZ(vtx); | |
240 | fvtxz->Fill(vtx[2]); | |
241 | ||
10eaad41 | 242 | for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc) |
243 | { | |
244 | TParticle* particle = stack->Particle(iMc); | |
245 | ||
246 | if (!particle) | |
247 | { | |
1ac5c71d | 248 | if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc); |
10eaad41 | 249 | continue; |
250 | } | |
251 | ||
7ccf0419 | 252 | if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) { |
04c3c355 | 253 | |
254 | if(particle->GetNDaughters()>2) continue; | |
7ccf0419 | 255 | |
10eaad41 | 256 | Int_t firstD=particle->GetFirstDaughter(); |
257 | Int_t lastD=particle->GetLastDaughter(); | |
7ccf0419 | 258 | |
259 | if((firstD<0)||(firstD>=stack->GetNtrack())) continue; | |
260 | if((lastD<0)||(lastD>=stack->GetNtrack())) continue; | |
10eaad41 | 261 | |
eaef9478 | 262 | TParticle *daughterParticle1 = 0; |
263 | TParticle *daughterParticle2 = 0; | |
264 | AliMCParticle *mcDaughter1 = 0; | |
265 | AliMCParticle *mcDaughter2 = 0; | |
7ccf0419 | 266 | |
267 | if(fdaughter1pdg==kdaughterKaon) { | |
268 | daughterParticle1=stack->Particle(firstD); | |
269 | daughterParticle2=stack->Particle(lastD); | |
270 | mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD); | |
271 | mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD); | |
10eaad41 | 272 | } |
79b67623 | 273 | else |
7ccf0419 | 274 | if(fdaughter2pdg==kdaughterKaon) { |
275 | daughterParticle1=stack->Particle(lastD); | |
276 | daughterParticle2=stack->Particle(firstD); | |
277 | mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD); | |
278 | mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(firstD); | |
279 | } //to ensure that the first daughter is always the kaon | |
79b67623 | 280 | |
281 | if(TMath::Abs(daughterParticle1->GetPdgCode())!=321) continue; | |
282 | ||
eaef9478 | 283 | TParticle *daughters1Daughter = 0; |
284 | TParticle *daughters2Daughter = 0; | |
285 | Int_t mcProcessDaughters1Daughter = -999; | |
286 | Int_t mcProcessDaughters2Daughter = -999; | |
287 | AliMCParticle *mcDaughters1Daughter = 0; | |
7ccf0419 | 288 | |
289 | if(mcDaughter1->Charge()==0) continue; | |
290 | if(mcDaughter2->Charge()==0) continue; //accept resonance decays in two charged daughters | |
291 | ||
eaef9478 | 292 | Int_t nDecayKaonDaughter = -99; |
7ccf0419 | 293 | for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) { |
294 | if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) { | |
295 | daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia); | |
296 | mcProcessDaughters1Daughter=daughters1Daughter->GetUniqueID(); | |
297 | if(mcProcessDaughters1Daughter==4) { | |
298 | nDecayKaonDaughter=daughterParticle1->GetFirstDaughter()+ia; | |
299 | break; | |
300 | } | |
301 | } | |
302 | } | |
303 | ||
eaef9478 | 304 | Int_t nProcessDaughter = -99; |
7ccf0419 | 305 | for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) { |
306 | if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) { | |
307 | daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib); | |
308 | mcProcessDaughters2Daughter=daughters2Daughter->GetUniqueID(); | |
309 | if((mcProcessDaughters2Daughter==4)||(mcProcessDaughters2Daughter==13)) { | |
310 | nProcessDaughter=mcProcessDaughters2Daughter; | |
311 | break; | |
312 | } | |
313 | } | |
314 | } | |
315 | ||
eaef9478 | 316 | Int_t numberOfCharged = 0; |
7ccf0419 | 317 | if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) { |
318 | for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) { | |
319 | if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic)); | |
320 | else continue; | |
321 | if(mcDaughters1Daughter->Charge()!=0) numberOfCharged=numberOfCharged+1; | |
322 | } | |
323 | } | |
324 | ||
325 | if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) { | |
326 | ||
10eaad41 | 327 | fSimEta->Fill(particle->Eta()); |
10eaad41 | 328 | fSimEtaPt->Fill(particle->Pt(), particle->Eta()); |
7ccf0419 | 329 | fSimPt->Fill(particle->Pt()); |
330 | fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt()); | |
331 | if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter | |
332 | ||
333 | ||
334 | if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) { | |
335 | if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable | |
336 | fSimPtKink->Fill(particle->Pt()); | |
337 | fSimEtaKink->Fill(particle->Eta()); | |
338 | fSimEtaPtKink->Fill(particle->Pt(), particle->Eta()); | |
339 | fMCInvmassPt->Fill(particle->GetMass(), particle->Pt()); | |
340 | } | |
341 | } | |
342 | ||
343 | } // for the generated spectra | |
344 | ||
345 | } //for the particular resonance | |
346 | } //MC loop | |
10eaad41 | 347 | |
348 | } // end fAnalysisType==MC | |
349 | else | |
350 | ||
351 | if(fAnalysisType == "ESD") { | |
f9afc48d | 352 | const AliESDVertex* vertex = GetEventVertex(esd); |
10eaad41 | 353 | if(!vertex) return; |
354 | Double_t vtx[3]; | |
355 | vertex->GetXYZ(vtx); | |
356 | fvtxz->Fill(vtx[2]); | |
10eaad41 | 357 | Double_t ptrackpos[3], ptrackneg[3]; |
358 | ||
359 | TLorentzVector p4pos, anp4pos; | |
360 | TLorentzVector p4neg, anp4neg; | |
361 | TLorentzVector p4comb, anp4comb; | |
362 | ||
f9afc48d | 363 | for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) { |
364 | AliESDtrack* trackpos = esd->GetTrack(iTracks); | |
10eaad41 | 365 | if (!trackpos) { |
1ac5c71d | 366 | if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks); |
10eaad41 | 367 | continue; |
368 | } | |
369 | if (trackpos->GetSign() < 0) continue; | |
370 | ||
92adf4f6 | 371 | AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam(); |
372 | if(!tpcTrackpos) continue; | |
373 | ptrackpos[0]=tpcTrackpos->Px(); | |
374 | ptrackpos[1]=tpcTrackpos->Py(); | |
375 | ptrackpos[2]=tpcTrackpos->Pz(); | |
10eaad41 | 376 | |
92adf4f6 | 377 | Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos); |
378 | if(firstLevelAcceptPosTrack==kFALSE) continue; | |
10eaad41 | 379 | |
380 | TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]); | |
92adf4f6 | 381 | |
10eaad41 | 382 | TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel())); |
383 | if (!partpos) continue; | |
384 | Int_t pdgpos = partpos->GetPdgCode(); | |
385 | Int_t mumlabelpos=partpos->GetFirstMother(); | |
386 | mumlabelpos = TMath::Abs(mumlabelpos); | |
387 | TParticle * mumpos=stack->Particle(mumlabelpos); | |
388 | if (!mumpos) continue; | |
389 | Int_t mumpdgpos = mumpos->GetPdgCode(); | |
390 | ||
391 | Int_t indexKinkPos=trackpos->GetKinkIndex(0); | |
92adf4f6 | 392 | Bool_t posKaonKinkFlag=0; |
393 | if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom); | |
10eaad41 | 394 | |
92adf4f6 | 395 | if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass); |
10eaad41 | 396 | |
397 | if(indexKinkPos==0) { | |
92adf4f6 | 398 | |
399 | Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos); | |
400 | if(secondLevelAcceptPosTrack==kFALSE) continue; | |
401 | ||
f9afc48d | 402 | p4pos.SetVectM(posTrackMom, daughter2pdgMass); |
10eaad41 | 403 | |
404 | } | |
405 | ||
f9afc48d | 406 | for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) { |
10eaad41 | 407 | if(iTracks==j) continue; |
f9afc48d | 408 | AliESDtrack* trackneg=esd->GetTrack(j); |
10eaad41 | 409 | if (trackneg->GetSign() > 0) continue; |
410 | ||
92adf4f6 | 411 | AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam(); |
412 | if(!tpcTrackneg) continue; | |
413 | ptrackneg[0]=tpcTrackneg->Px(); | |
414 | ptrackneg[1]=tpcTrackneg->Py(); | |
415 | ptrackneg[2]=tpcTrackneg->Pz(); | |
10eaad41 | 416 | |
92adf4f6 | 417 | Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg); |
418 | if(firstLevelAcceptNegTrack==kFALSE) continue; | |
10eaad41 | 419 | |
420 | TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]); | |
10eaad41 | 421 | |
422 | TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel())); | |
423 | if (!partneg) continue; | |
424 | Int_t pdgneg = partneg->GetPdgCode(); | |
425 | Int_t mumlabelneg=partneg->GetFirstMother(); | |
426 | mumlabelneg = TMath::Abs(mumlabelneg); | |
427 | TParticle * mumneg=stack->Particle(mumlabelneg); | |
428 | if (!mumneg) continue; | |
429 | Int_t mumpdgneg = mumneg->GetPdgCode(); | |
430 | ||
431 | Int_t indexKinkNeg=trackneg->GetKinkIndex(0); | |
92adf4f6 | 432 | Bool_t negKaonKinkFlag=0; |
433 | if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom); | |
10eaad41 | 434 | |
f9afc48d | 435 | if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass); |
10eaad41 | 436 | |
437 | if(indexKinkNeg==0) { | |
92adf4f6 | 438 | |
439 | Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg); | |
440 | if(secondLevelAcceptNegTrack==kFALSE) continue; | |
441 | ||
f9afc48d | 442 | anp4neg.SetVectM(negTrackMom, daughter2pdgMass); |
10eaad41 | 443 | |
444 | } | |
445 | ||
446 | Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag()); | |
447 | ||
92adf4f6 | 448 | if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) { |
10eaad41 | 449 | p4comb=anp4pos; |
450 | p4comb+=p4neg; | |
451 | if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M()); | |
452 | } | |
453 | ||
454 | if(negKaonKinkFlag==1) { | |
455 | p4comb=p4pos; | |
456 | p4comb+=p4neg; | |
457 | fInvariantMass->Fill(p4comb.M()); | |
7ccf0419 | 458 | fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt()); |
f9afc48d | 459 | if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg) |
460 | &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) { | |
10eaad41 | 461 | fOpeningAngle->Fill(openingAngle); |
462 | fInvMassTrue->Fill(p4comb.M()); | |
7ccf0419 | 463 | fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt()); |
10eaad41 | 464 | if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) { |
465 | fRecPt->Fill(p4comb.Vect().Pt()); | |
466 | fRecEta->Fill(p4comb.Vect().Eta()); | |
467 | fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta()); | |
468 | ||
469 | } | |
470 | ||
471 | } | |
472 | ||
473 | } | |
474 | ||
92adf4f6 | 475 | if(posKaonKinkFlag==1) { |
10eaad41 | 476 | anp4comb=anp4pos; |
477 | anp4comb+=anp4neg; | |
478 | fInvariantMass->Fill(anp4comb.M()); | |
7ccf0419 | 479 | fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt()); |
f9afc48d | 480 | if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg) |
481 | &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) { | |
10eaad41 | 482 | fOpeningAngle->Fill(openingAngle); |
483 | fInvMassTrue->Fill(p4comb.M()); | |
7ccf0419 | 484 | fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt()); |
10eaad41 | 485 | if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) { |
486 | fRecPt->Fill(anp4comb.Vect().Pt()); | |
487 | fRecEta->Fill(anp4comb.Vect().Eta()); | |
488 | fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta()); | |
489 | } | |
490 | ||
491 | } | |
492 | ||
493 | } | |
494 | ||
495 | } //inner track loop | |
496 | ||
497 | } //outer track loop | |
498 | ||
499 | } // end fAnalysisType == ESD | |
500 | ||
10eaad41 | 501 | } |
502 | ||
10eaad41 | 503 | //____________________________________________________________________// |
f9afc48d | 504 | Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const { |
10eaad41 | 505 | // Calculates the number of sigma to the vertex. |
506 | ||
507 | Float_t b[2]; | |
508 | Float_t bRes[2]; | |
509 | Float_t bCov[3]; | |
510 | ||
92adf4f6 | 511 | esdTrack->GetImpactParametersTPC(b,bCov); |
10eaad41 | 512 | |
513 | if (bCov[0]<=0 || bCov[2]<=0) { | |
10eaad41 | 514 | bCov[0]=0; bCov[2]=0; |
515 | } | |
f9afc48d | 516 | |
10eaad41 | 517 | bRes[0] = TMath::Sqrt(bCov[0]); |
518 | bRes[1] = TMath::Sqrt(bCov[2]); | |
519 | ||
520 | if (bRes[0] == 0 || bRes[1] ==0) return -1; | |
521 | ||
522 | Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2)); | |
523 | ||
524 | if (TMath::Exp(-d * d / 2) < 1e-10) return 1000; | |
525 | ||
526 | d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2); | |
527 | ||
528 | return d; | |
529 | } | |
530 | ||
f9afc48d | 531 | //________________________________________________________________________ |
532 | const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const | |
10eaad41 | 533 | { |
534 | // Get the vertex | |
535 | ||
536 | const AliESDVertex* vertex = esd->GetPrimaryVertex(); | |
537 | ||
538 | if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex; | |
539 | else | |
540 | { | |
541 | vertex = esd->GetPrimaryVertexSPD(); | |
542 | if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex; | |
543 | else | |
544 | return 0; | |
545 | } | |
546 | } | |
547 | ||
92adf4f6 | 548 | //________________________________________________________________________ |
549 | ||
550 | Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd, | |
551 | const AliESDVertex *localvertex, AliESDtrack* localtrack) { | |
552 | // Apply the selections for each kink | |
553 | ||
554 | Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0; | |
555 | Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance. | |
556 | Double_t dca3D = 0.0; | |
557 | ||
558 | AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam(); | |
559 | if(!tpcTrack) { | |
560 | gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; | |
561 | dca[0] = -100.; dca[1] = -100.; dca3D = -100.; | |
562 | cov[0] = -100.; cov[1] = -100.; cov[2] = -100.; | |
563 | } | |
564 | else { | |
565 | gPt = tpcTrack->Pt(); | |
566 | gPx = tpcTrack->Px(); | |
567 | gPy = tpcTrack->Py(); | |
568 | gPz = tpcTrack->Pz(); | |
569 | tpcTrack->PropagateToDCA(localvertex, | |
570 | localesd->GetMagneticField(),100.,dca,cov); | |
571 | } | |
572 | ||
573 | if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) { | |
1ac5c71d | 574 | if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex); |
92adf4f6 | 575 | return kFALSE; |
576 | } | |
577 | ||
578 | if(TMath::Abs(dca[0]) > fMaxDCAxy) { | |
1ac5c71d | 579 | if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)", TMath::Abs(dca[0]), fMaxDCAxy); |
92adf4f6 | 580 | return kFALSE; |
581 | } | |
582 | ||
583 | if(TMath::Abs(dca[1]) > fMaxDCAzaxis) { | |
1ac5c71d | 584 | if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)", TMath::Abs(dca[1]), fMaxDCAzaxis); |
92adf4f6 | 585 | return kFALSE; |
586 | } | |
587 | ||
588 | if(gPt < fMinPtTrackCut) { | |
1ac5c71d | 589 | if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut); |
92adf4f6 | 590 | return kFALSE; |
591 | } | |
592 | ||
593 | return kTRUE; | |
594 | } | |
595 | ||
596 | //________________________________________________________________________ | |
597 | Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd, const AliESDVertex *localvertex, AliESDtrack *localtrack) { | |
598 | // Apply the selections for each track | |
f9afc48d | 599 | |
92adf4f6 | 600 | Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0; |
601 | Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance. | |
602 | Double_t dca3D = 0.0; | |
603 | ||
604 | AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam(); | |
605 | if(!tpcTrack) { | |
606 | gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; | |
607 | dca[0] = -100.; dca[1] = -100.; dca3D = -100.; | |
608 | cov[0] = -100.; cov[1] = -100.; cov[2] = -100.; | |
609 | } | |
610 | else { | |
611 | gPt = tpcTrack->Pt(); | |
612 | gPx = tpcTrack->Px(); | |
613 | gPy = tpcTrack->Py(); | |
614 | gPz = tpcTrack->Pz(); | |
615 | tpcTrack->PropagateToDCA(localvertex, | |
616 | localesd->GetMagneticField(),100.,dca,cov); | |
617 | } | |
618 | ||
619 | Int_t fcls[200]; | |
620 | Int_t nClustersTPC=localtrack->GetTPCclusters(fcls); | |
621 | Float_t chi2perTPCcluster=-1.0; | |
622 | if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC); | |
623 | ||
624 | Double_t extCov[15]; | |
625 | localtrack->GetExternalCovariance(extCov); | |
626 | ||
627 | if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) { | |
1ac5c71d | 628 | if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC"); |
92adf4f6 | 629 | return kFALSE; |
630 | } | |
631 | ||
632 | if(nClustersTPC < fMinTPCclusters) { | |
1ac5c71d | 633 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters); |
92adf4f6 | 634 | return kFALSE; |
635 | } | |
636 | ||
637 | if(chi2perTPCcluster > fMaxChi2PerTPCcluster) { | |
1ac5c71d | 638 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster); |
92adf4f6 | 639 | return kFALSE; |
640 | } | |
641 | ||
642 | if(extCov[0] > fMaxCov0) { | |
1ac5c71d | 643 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0); |
92adf4f6 | 644 | return kFALSE; |
645 | } | |
646 | ||
647 | if(extCov[2] > fMaxCov2) { | |
1ac5c71d | 648 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2); |
92adf4f6 | 649 | return kFALSE; |
650 | } | |
651 | ||
652 | if(extCov[5] > fMaxCov5) { | |
1ac5c71d | 653 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5); |
92adf4f6 | 654 | return kFALSE; |
655 | } | |
656 | ||
657 | if(extCov[9] > fMaxCov9) { | |
1ac5c71d | 658 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9); |
92adf4f6 | 659 | return kFALSE; |
660 | } | |
661 | ||
662 | if(extCov[14] > fMaxCov14) { | |
1ac5c71d | 663 | if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14); |
92adf4f6 | 664 | return kFALSE; |
665 | } | |
f9afc48d | 666 | |
92adf4f6 | 667 | return kTRUE; |
668 | ||
669 | } | |
670 | ||
671 | //________________________________________________________________________ | |
672 | Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom) | |
673 | { | |
674 | // Test some kinematical criteria for each kink | |
675 | ||
676 | AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1); | |
677 | const TVector3 motherMfromKink(kink->GetMotherP()); | |
678 | const TVector3 daughterMKink(kink->GetDaughterP()); | |
679 | Float_t qt=kink->GetQt(); | |
680 | ||
681 | Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.); | |
682 | Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.); | |
683 | ||
684 | Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2); | |
685 | ||
686 | Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658); | |
687 | Float_t p1XM= motherMfromKink.Px(); | |
688 | Float_t p1YM= motherMfromKink.Py(); | |
689 | Float_t p1ZM= motherMfromKink.Pz(); | |
690 | Float_t p2XM= daughterMKink.Px(); | |
691 | Float_t p2YM= daughterMKink.Py(); | |
692 | Float_t p2ZM= daughterMKink.Pz(); | |
693 | Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM))); | |
694 | Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag()); | |
695 | ||
696 | if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) { | |
697 | ||
698 | if(trackMom.Mag()<=1.1) { | |
699 | return kTRUE; | |
700 | } | |
701 | else | |
702 | if (kinkAngle<maxDecAngKmu) { | |
703 | return kTRUE; | |
704 | } | |
705 | } | |
706 | return kFALSE; | |
707 | } |