Add fast merging option (Diego)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronHelper.cxx
CommitLineData
ffbede40 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// Dielectron helper functions wrapped in a namespace
18//
19//
20// Authors:
21// Jens Wiechula <Jens.Wiechula@cern.ch>
ba15fdfb 22// Frederick Kramer <Frederick.Kramer@cern.ch>
ffbede40 23//
24
25
26
27
28#include <TError.h>
29#include <TMath.h>
30#include <TObjString.h>
31#include <TObjArray.h>
32#include <TVectorD.h>
33
1201a1a9 34#include <AliVEvent.h>
ba15fdfb 35#include <AliVParticle.h>
1201a1a9 36#include <AliKFParticle.h>
ba15fdfb 37#include <AliESDtrackCuts.h>
38#include <AliESDEvent.h>
39#include <AliMCEvent.h>
40#include <AliAODEvent.h>
41#include <AliAODTracklets.h>
42#include <AliMultiplicity.h>
43#include <AliStack.h>
1201a1a9 44
ffbede40 45#include "AliDielectronHelper.h"
46
47//_____________________________________________________________________________
48TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
49{
50 //
51 // Make logarithmic binning
52 // the user has to delete the array afterwards!!!
53 //
54
55 //check limits
56 if (xmin<1e-20 || xmax<1e-20){
57 Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
58 return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax);
59 }
60 if (xmax<xmin){
61 Double_t tmp=xmin;
62 xmin=xmax;
63 xmax=tmp;
64 }
65 TVectorD *binLim=new TVectorD(nbinsX+1);
66 Double_t first=xmin;
67 Double_t last=xmax;
68 Double_t expMax=TMath::Log(last/first);
69 for (Int_t i=0; i<nbinsX+1; ++i){
70 (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
71 }
72 return binLim;
73}
74
75//_____________________________________________________________________________
76TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
77{
78 //
79 // Make linear binning
80 // the user has to delete the array afterwards!!!
81 //
82 if (xmax<xmin){
83 Double_t tmp=xmin;
84 xmin=xmax;
85 xmax=tmp;
86 }
87 TVectorD *binLim=new TVectorD(nbinsX+1);
88 Double_t first=xmin;
89 Double_t last=xmax;
90 Double_t binWidth=(last-first)/nbinsX;
91 for (Int_t i=0; i<nbinsX+1; ++i){
92 (*binLim)[i]=first+binWidth*(Double_t)i;
93 }
94 return binLim;
95}
96
97//_____________________________________________________________________________
98TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
99{
100 //
101 // Make arbitrary binning, bins separated by a ','
102 //
103 TString limits(bins);
104 if (limits.IsNull()){
105 Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
106 return 0x0;
107 }
108
109 TObjArray *arr=limits.Tokenize(",");
110 Int_t nLimits=arr->GetEntries();
111 if (nLimits<2){
112 Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
113 delete arr;
114 return 0x0;
115 }
116
117 TVectorD *binLimits=new TVectorD(nLimits);
118 for (Int_t iLim=0; iLim<nLimits; ++iLim){
119 (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
120 }
121
122 delete arr;
123 return binLimits;
124}
125
ba15fdfb 126
127//_____________________________________________________________________________
128Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
129 // determination of Nch
130 if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
131
132 AliStack *stack = ((AliMCEvent*)ev)->Stack();
133
134 if (!stack) return -1;
135
136 Int_t nParticles = stack->GetNtrack();
137 Int_t nCh = 0;
138
139 // count..
140 for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
141 if (!stack->IsPhysicalPrimary(iMc)) continue;
142
143 TParticle* particle = stack->Particle(iMc);
144 if (!particle) continue;
145 if (particle->GetPDG()->Charge() == 0) continue;
146
147 Float_t eta = particle->Eta();
148 if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
149 }
150
151 return nCh;
152}
153
154
155//_____________________________________________________________________________
156Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev){
157 // Compute the collision multiplicity based on AOD or ESD tracklets
158 // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
159
160 if (!ev) return -1;
161
162 Int_t nTracklets = 0;
163 Int_t nAcc = 0;
164 Double_t etaRange = 1.6;
165
166 if (ev->IsA() == AliAODEvent::Class()) {
167 AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
168 nTracklets = tracklets->GetNumberOfTracklets();
169 for (Int_t nn = 0; nn < nTracklets; nn++) {
170 Double_t theta = tracklets->GetTheta(nn);
171 Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
172 if (TMath::Abs(eta) < etaRange) nAcc++;
173 }
174 } else if (ev->IsA() == AliESDEvent::Class()) {
175 nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
176 for (Int_t nn = 0; nn < nTracklets; nn++) {
177 Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
178 if (TMath::Abs(eta) < etaRange) nAcc++;
179 }
180 } else return -1;
181
182 return nAcc;
183}
184
185
186//_____________________________________________________________________________
5cc8c825 187Int_t AliDielectronHelper::GetNacc(const AliVEvent */*ev*/){
ba15fdfb 188 // put a robust Nacc definition here
189
190 return -1;
5cc8c825 191/*
ba15fdfb 192 if (!ev || ev->IsA()!=AliESDEvent::Class()) return -1;
193
194 // basic track cuts for the N_acc definition
195 AliESDtrackCuts esdTC;
196 esdTC.SetMaxDCAToVertexZ(3.0);
197 esdTC.SetMaxDCAToVertexXY(1.0);
198 esdTC.SetEtaRange( -0.9 , 0.9 );
199 esdTC.SetAcceptKinkDaughters(kFALSE);
200 esdTC.SetRequireITSRefit(kTRUE);
201 esdTC.SetRequireTPCRefit(kTRUE);
202 esdTC.SetMinNClustersTPC(70);
203 esdTC.SetMaxChi2PerClusterTPC(4);
204 esdTC.SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
205
206 Int_t nRecoTracks = ev->GetNumberOfTracks();
207 Int_t nAcc = 0;
208
209 for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
210 AliVParticle* candidate = ev->GetTrack(iTrack);
211 if (esdTC.IsSelected(candidate)) nAcc++;
212 }
213
214 return nAcc;
5cc8c825 215*/
ba15fdfb 216}
217
218
1201a1a9 219//_____________________________________________________________________________
220void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
221 // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
222 if (!kfParticle) return;
223 Double_t dx = 0.;
224 Double_t dy = 0.;
225 Double_t dz = 0.;
226
227 if (ev){
228 dx = ev->GetPrimaryVertex()->GetX()-0.;
229 dy = ev->GetPrimaryVertex()->GetY()-0.;
230 dz = ev->GetPrimaryVertex()->GetZ()-0.;
231 }
232
233 kfParticle->X() = kfParticle->GetX() - dx;
234 kfParticle->Y() = kfParticle->GetY() - dy;
235 kfParticle->Z() = kfParticle->GetZ() - dz;
236
237
238 // Rotate the kf particle
239 Double_t c = cos(angle);
240 Double_t s = sin(angle);
241
242 Double_t mA[8][ 8];
243 for( Int_t i=0; i<8; i++ ){
244 for( Int_t j=0; j<8; j++){
245 mA[i][j] = 0;
246 }
247 }
248 for( int i=0; i<8; i++ ){
249 mA[i][i] = 1;
250 }
251 mA[0][0] = c; mA[0][1] = s;
252 mA[1][0] = -s; mA[1][1] = c;
253 mA[3][3] = c; mA[3][4] = s;
254 mA[4][3] = -s; mA[4][4] = c;
255
256 Double_t mAC[8][8];
257 Double_t mAp[8];
258
259 for( Int_t i=0; i<8; i++ ){
260 mAp[i] = 0;
261 for( Int_t k=0; k<8; k++){
262 mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
263 }
264 }
265
266 for( Int_t i=0; i<8; i++){
267 kfParticle->Parameter(i) = mAp[i];
268 }
269
270 for( Int_t i=0; i<8; i++ ){
271 for( Int_t j=0; j<8; j++ ){
272 mAC[i][j] = 0;
273 for( Int_t k=0; k<8; k++ ){
274 mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
275 }
276 }
277 }
278
279 for( Int_t i=0; i<8; i++ ){
280 for( Int_t j=0; j<=i; j++ ){
281 Double_t xx = 0;
282 for( Int_t k=0; k<8; k++){
283 xx+= mAC[i][k]*mA[j][k];
284 }
285 kfParticle->Covariance(i,j) = xx;
286 }
287 }
288
289 kfParticle->X() = kfParticle->GetX() + dx;
290 kfParticle->Y() = kfParticle->GetY() + dy;
291 kfParticle->Z() = kfParticle->GetZ() + dz;
292
293}
294