making explicit the file format when saving a Canvas
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKinkLikeSign.cxx
CommitLineData
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//----------------------------------------------------------------------------------------------------------------
15// class AliResonanceKinkLikeSign
16// Example of an analysis task for producing a like-sign background for resonances having at least one
17// kaon-kink in their decay products.
7ccf0419 18// Background is calculated from a positive kaon kink and a negative track.
10eaad41 19//-----------------------------------------------------------------------------------------------------------------
20
21#include "TChain.h"
22#include "TTree.h"
23#include "TVector3.h"
24#include "TF1.h"
25#include "TH1D.h"
7ccf0419 26#include "TH2D.h"
f27d6e67 27#include <TDatabasePDG.h>
28#include <TParticlePDG.h>
10eaad41 29#include "TLorentzVector.h"
30#include "AliAnalysisTaskSE.h"
31#include "AliAnalysisManager.h"
32
33#include "AliESDInputHandler.h"
34
35#include "AliResonanceKinkLikeSign.h"
36#include "AliESDkink.h"
37
38ClassImp(AliResonanceKinkLikeSign)
39
40//________________________________________________________________________
41AliResonanceKinkLikeSign::AliResonanceKinkLikeSign()
7ccf0419 42 : AliAnalysisTaskSE(), fESD(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0)
10eaad41 43{
44 // Constructor
45}
46
47//________________________________________________________________________
48AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name)
7ccf0419 49 : AliAnalysisTaskSE(name), fESD(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0)
10eaad41 50
51{
52 // Constructor
53
54 // Define input and output slots here
55 // Input slot #0 works with a TChain
56 DefineInput(0, TChain::Class());
57 DefineOutput(1, TList::Class());
58}
59
60//________________________________________________________________________
6dd40f14 61void AliResonanceKinkLikeSign::ConnectInputData(Option_t *option)
10eaad41 62{
63 // Connect ESD or AOD here
64 // Called once
6dd40f14 65 AliAnalysisTaskSE::ConnectInputData(option);
10eaad41 66 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
67 if (!tree) {
68 Printf("ERROR: Could not read chain from input slot 0");
69 } else {
f27d6e67 70
10eaad41 71 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
72
73 if (!esdH) {
74 Printf("ERROR: Could not get ESDInputHandler");
75 } else
76 fESD = esdH->GetEvent();
77 }
78}
79
80//________________________________________________________________________
81void AliResonanceKinkLikeSign::CreateOutputObjects()
82{
83 // Create histograms
84 // Called once
85
86 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);
87 f1->SetParameter(0,0.493677);
88 f1->SetParameter(1,0.9127037);
89 f1->SetParameter(2,TMath::Pi());
90
91 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);
92 f2->SetParameter(0,0.13957018);
93 f2->SetParameter(1,0.2731374);
94 f2->SetParameter(2,TMath::Pi());
95
96 //OpenFile(1); //uncomment for proof analysis
97
98 // for K*0(892)
7ccf0419 99 fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", 60,0.6,1.2);
100 fLikeSignInvmassPt=new TH2D("fLikeSignInvmassPt"," ", 60,0.6,1.2, 100,0.0,10.0);
10eaad41 101
102 // for phi(1020)
7ccf0419 103 // fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", 70,0.99,1.088);
10eaad41 104
105
106 fListOfHistos=new TList();
7ccf0419 107 fListOfHistos->Add(fPosKaonLikeSign);
108 fListOfHistos->Add(fLikeSignInvmassPt);
10eaad41 109
110}
111
112//________________________________________________________________________
6dd40f14 113void AliResonanceKinkLikeSign::Exec(Option_t *option)
10eaad41 114{
115 // Main loop
116 // Called for each event
117
f27d6e67 118 Double_t daughter1Mass=TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass();
119 Double_t daughter2Mass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
10eaad41 120
f2a16059 121 if (option) Printf("WARNING: an option was initially required");
122
10eaad41 123 if (!fESD) {
124 Printf("ERROR: fESD not available");
125 return;
126 }
127
128 const AliESDVertex* vertex = GetEventVertex(fESD);
129 if(!vertex) return;
130
131 Double_t ptrackpos[3], ptrackneg[3];
132
133 TLorentzVector p4pos;
134 TLorentzVector p4neg;
135 TLorentzVector p4comb;
136
137 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
138 AliESDtrack* trackpos = fESD->GetTrack(iTracks);
139 if (!trackpos) {
1ac5c71d 140 if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
10eaad41 141 continue;
142 }
143
7ccf0419 144 if (trackpos->GetSign() < 0) continue;
10eaad41 145
146 trackpos->GetPxPyPz(ptrackpos);
147
148 Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
149
150 Float_t bpos[2];
151 Float_t bCovpos[3];
152 trackpos->GetImpactParameters(bpos,bCovpos);
153
154 if (bCovpos[0]<=0 || bCovpos[2]<=0) {
1ac5c71d 155 if (fDebug > 0) Printf("Estimated b resolution lower or equal zero!");
10eaad41 156 bCovpos[0]=0; bCovpos[2]=0;
157 }
158
159 Float_t dcaToVertexXYpos = bpos[0];
160 Float_t dcaToVertexZpos = bpos[1];
161
162 if(nSigmaToVertex>=4) continue;
163 if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
164
165 TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
166
167 if(posTrackMom.Perp()<=0.25) continue;
168
7ccf0419 169 //Uncomment the following block if the Like Sign is made of K+ kink + positive track
10eaad41 170
04c3c355 171 Int_t indexKink=trackpos->GetKinkIndex(0);
7ccf0419 172 Int_t kaonKinkFlag=0;
04c3c355 173 if(indexKink<0){
7ccf0419 174
04c3c355 175 AliESDkink *kink=fESD->GetKink(TMath::Abs(indexKink)-1);
7ccf0419 176 const TVector3 motherMfromKink(kink->GetMotherP());
177 const TVector3 daughterMKink(kink->GetDaughterP());
178 Float_t qt=kink->GetQt();
179
180 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
181 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
182
183 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
184
185 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
186 Float_t p1XM= motherMfromKink.Px();
187 Float_t p1YM= motherMfromKink.Py();
188 Float_t p1ZM= motherMfromKink.Pz();
189 Float_t p2XM= daughterMKink.Px();
190 Float_t p2YM= daughterMKink.Py();
191 Float_t p2ZM= daughterMKink.Pz();
192 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
193 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
194
195 if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
196
197 if(posTrackMom.Mag()<=1.1) {
198 kaonKinkFlag=1;
199 }
200 else
201 if (kinkAngle<maxDecAngKmu) {
202 kaonKinkFlag=1;
203 }
204 }
205
206 } //End Kink Information
207
208 if(kaonKinkFlag==0) continue;
209 if(kaonKinkFlag==1) p4pos.SetVectM(posTrackMom,daughter1Mass);
210
211 // Comment the following statements till the "for" if the Like Sign of K+ kink + positive track is needed
212
213// UInt_t status=trackpos->GetStatus();
214// if((status&AliESDtrack::kTPCrefit)==0) continue;
215// if(trackpos->GetTPCclusters(0)<50) continue;
216// if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
217// Double_t extCovPos[15];
218// trackpos->GetExternalCovariance(extCovPos);
219// if(extCovPos[0]>2) continue;
220// if(extCovPos[2]>2) continue;
221// if(extCovPos[5]>0.5) continue;
222// if(extCovPos[9]>0.5) continue;
223// if(extCovPos[14]>2) continue;
10eaad41 224//
7ccf0419 225// p4pos.SetVectM(posTrackMom,daughter1Mass);
10eaad41 226
227 for (Int_t j=0; j<fESD->GetNumberOfTracks(); j++) {
228 if(iTracks==j) continue;
229 AliESDtrack* trackneg=fESD->GetTrack(j);
7ccf0419 230 if (trackneg->GetSign() < 0) continue;
10eaad41 231
232 trackneg->GetPxPyPz(ptrackneg);
233 Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
234
235 Float_t bneg[2];
236 Float_t bCovneg[3];
237 trackneg->GetImpactParameters(bneg,bCovneg);
238 if (bCovneg[0]<=0 || bCovneg[2]<=0) {
1ac5c71d 239 if (fDebug > 0) Printf("Estimated b resolution lower or equal zero!");
10eaad41 240 bCovneg[0]=0; bCovneg[2]=0;
241 }
242
243 Float_t dcaToVertexXYneg = bneg[0];
244 Float_t dcaToVertexZneg = bneg[1];
245
246 if(negSigmaToVertex>=4) continue;
247 if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
248
249 TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
250
251 if(negTrackMom.Perp()<=0.25) continue;
252
253 UInt_t statusneg=trackneg->GetStatus();
254
255 if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
256
257 if(trackneg->GetTPCclusters(0)<50) continue;
258 if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
259 Double_t extCovneg[15];
260 trackneg->GetExternalCovariance(extCovneg);
261 if(extCovneg[0]>2) continue;
262 if(extCovneg[2]>2) continue;
263 if(extCovneg[5]>0.5) continue;
264 if(extCovneg[9]>0.5) continue;
265 if(extCovneg[14]>2) continue;
266
267 p4neg.SetVectM(negTrackMom, daughter2Mass);
268
269 p4comb=p4pos;
270 p4comb+=p4neg;
7ccf0419 271 fPosKaonLikeSign->Fill(p4comb.M());
272 fLikeSignInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
10eaad41 273
274 } //inner track loop
275
276 } //outer track loop
277
278 PostData(1, fListOfHistos);
279}
280
10eaad41 281//____________________________________________________________________//
282
283Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
284 // Calculates the number of sigma to the vertex.
285
286 Float_t b[2];
287 Float_t bRes[2];
288 Float_t bCov[3];
289
290 esdTrack->GetImpactParameters(b,bCov);
291
292 if (bCov[0]<=0 || bCov[2]<=0) {
293 //AliDebug(1, "Estimated b resolution lower or equal zero!");
294 bCov[0]=0; bCov[2]=0;
295 }
296 bRes[0] = TMath::Sqrt(bCov[0]);
297 bRes[1] = TMath::Sqrt(bCov[2]);
298
299 if (bRes[0] == 0 || bRes[1] ==0) return -1;
300
301 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
302
303 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
304
305 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
306
307 return d;
308}
309
310//____________________________________________________________________//
311
312const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const
313
314{
315 // Get the vertex
316
317 const AliESDVertex* vertex = esd->GetPrimaryVertex();
318
319 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
320 else
321 {
322 vertex = esd->GetPrimaryVertexSPD();
323 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
324 else
325 return 0;
326 }
327}
328