1 /*************************************************************************
2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2008
18 //=============================================================================
19 // AliUnicorEvent-AliESD interface //
20 //=============================================================================
26 #include "AliESDVZERO.h"
27 #include "AliESDEvent.h"
28 #include "AliMultiplicity.h"
29 #include "AliUnicorEventAliceESD.h"
31 ClassImp(AliUnicorEventAliceESD)
33 //=============================================================================
34 AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fViper(0), fESD(esd)//, fPhysicsSelection(0)
38 // printf("%s object created\n",ClassName());
39 if (!fESD) fESD = new AliESDEvent();
40 // fPhysicsSelection = new AliPhysicsSelection();
42 // V0 percentile graph for centrality determination
44 // TFile::Open("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_137161_v4.root","read");
45 TFile::Open("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_137366_v2.root","read");
46 const TH1F *hi = (const TH1F*) gFile->Get("hmultV0_percentile");
47 //const TH1D *hi = (const TH1D*) gFile->Get("hNtracks_percentile");
48 //const TH1F *hi = (const TH1F*) gFile->Get("hNclusters1_percentile");
49 fViper = new TGraph((const TH1*) hi);
52 //=============================================================================
53 AliUnicorEventAliceESD::~AliUnicorEventAliceESD() {
58 //=============================================================================
59 Bool_t AliUnicorEventAliceESD::Good() const
63 // if (!fPhysicsSelection->IsCollisionCandidate(fESD)) return kFALSE;
64 const AliESDVertex *vtx = fESD->GetPrimaryVertex();
65 if (!vtx->GetStatus()) return kFALSE;
66 if (fabs(Zver())>1) return kFALSE;
67 if (NGoodParticles()<9) return kFALSE;
70 //=============================================================================
71 Double_t AliUnicorEventAliceESD::Centrality() const {
73 // centrality between 0 (central) and 1 (very peripheral)
75 // V0 multiplicity, not linearized
77 AliESDVZERO *v0 = fESD->GetVZEROData();
78 float multa = v0->GetMTotV0A();
79 float multc = v0->GetMTotV0C();
80 double cent = fViper->Eval(multa+multc)/100.0;
84 // double cent = fViper->Eval(((AliVEvent*) fESD)->GetNumberOfTracks())/100.0;
86 // number of clusters in second layer of SPD
88 // double nspd2 = fESD->GetMultiplicity()->GetNumberOfITSClusters(1);
89 // double zv = fESD->GetPrimaryVertexSPD()->GetZ();
90 // const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
92 // float corr = 1 + zv*(pars[1] + zv*pars[2]);
93 // nspd2 = corr>0? nspd2/corr : -1;
94 // double cent = fViper->Eval(nspd2)/100.0;
98 //=============================================================================
99 Bool_t AliUnicorEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const
101 // track cuts and particle id cuts; pidi=0 means take all species
102 // consider using the standard ESDcut
104 // track quality cuts
106 AliESDtrack *track = fESD->GetTrack(i);
107 if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0; // TPC refit
108 if (!track->IsOn(AliESDtrack::kITSrefit)) return 0; // ITS refit
109 if (track->GetTPCNcls() < 90) return 0; // number of TPC clusters
110 if (track->GetKinkIndex(0) > 0) return 0; // no kink daughters
111 const AliExternalTrackParam *tp = GetTrackParam(i);
112 if (!tp) return 0; // track param
113 if (fabs(tp->Eta())>0.8) return 0; // fiducial pseudorapidity
115 // double pi9 = TMath::Pi()/9.0;
116 // double eta = tp->Eta();
117 // double phi = ParticlePhi(i);
118 // if (eta>0 && phi>-8*pi9 && phi<-7*pi9) return 0; // A10
119 // if (eta>0 && phi> 1*pi9 && phi< 2*pi9) return 0; // A01
120 // if (tp->Pt()<0.2) return 0; // lower pt cutoff
123 track->GetImpactParametersTPC(r,z);
124 if (fabs(z)>1.0) return 0; // impact parameter in z
125 if (fabs(r)>1.0) return 0; // impact parameter in xy
128 //TBits shared = track->GetTPCSharedMap();
129 //if (shared.CountBits()) return 0; // no shared clusters; pragmatic but dangerous
131 if (pidi==0) return 1;
135 if (!track->IsOn(AliESDtrack::kTPCpid)) return 0;
136 Double_t p[AliPID::kSPECIES];
138 Int_t q = tp->Charge();
140 if (pidi == -211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==-1;
141 else if (pidi == 211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==1;
142 else if (pidi == -321) return p[AliPID::kKaon]>0.5 && q==-1;
143 else if (pidi == 321) return p[AliPID::kKaon]>0.5 && q==1;
144 else if (pidi == -2212) return p[AliPID::kProton]>0.5 && q==-1;
145 else if (pidi == 2212) return p[AliPID::kProton]>0.5 && q==1;
146 else if (pidi == -11) return p[AliPID::kElectron]>0.5 && q==1;
147 else if (pidi == 11) return p[AliPID::kElectron]>0.5 && q==-1;
150 //=============================================================================
151 Bool_t AliUnicorEventAliceESD::PairGood(double p0, double the0, double phi0, double z0,
152 double p1, double the1, double phi1, double z1) const {
154 // two-track separation cut
156 double dthe = the1-the0;
157 if (fabs(dthe)>0.010) return kTRUE;
158 double B = -0.5; // magnetic field in T, needed for helix; should be gotten in the proper way
159 double pt0 = p0*sin(the0);
160 double pt1 = p1*sin(the1);
161 double r = 1.2; // we will calculate the distance between the two tracks at r=1.2 m
162 // for (double r=0.8; r<2.5; r+=0.2) {
163 double si0 = -0.3*B*z0*r/2/pt0;
164 double si1 = -0.3*B*z1*r/2/pt1;
165 if (fabs(si0)>=1.0) return kTRUE; // could be done better
166 if (fabs(si1)>=1.0) return kTRUE;
167 double dphi = phi1 - phi0 + asin(si1) - asin(si0);
168 dphi = TVector2::Phi_mpi_pi(dphi);
169 if (fabs(dphi)<0.020) return kFALSE;
173 //=============================================================================