1 /*************************************************************************
2 * Copyright(c) 1998-2009, 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 **************************************************************************/
17 // Dielectron helper functions wrapped in a namespace
21 // Jens Wiechula <Jens.Wiechula@cern.ch>
22 // Frederick Kramer <Frederick.Kramer@cern.ch>
23 // Julian Book <Julian.Book@cern.ch>
30 #include <TObjString.h>
31 #include <TObjArray.h>
37 #include <AliVEvent.h>
38 #include <AliVParticle.h>
39 #include <AliKFParticle.h>
40 #include <AliESDtrackCuts.h>
41 #include <AliESDEvent.h>
42 #include <AliMCEvent.h>
43 #include <AliAODEvent.h>
44 #include <AliAODTracklets.h>
45 #include <AliMultiplicity.h>
48 #include "AliDielectronVarCuts.h"
49 #include "AliDielectronTrackCuts.h"
50 #include "AliDielectronVarManager.h"
51 #include "AliDielectronHelper.h"
53 //_____________________________________________________________________________
54 TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
57 // Make logarithmic binning
58 // the user has to delete the array afterwards!!!
62 if (xmin<1e-20 || xmax<1e-20){
63 Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
64 return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax);
71 TVectorD *binLim=new TVectorD(nbinsX+1);
74 Double_t expMax=TMath::Log(last/first);
75 for (Int_t i=0; i<nbinsX+1; ++i){
76 (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
81 //_____________________________________________________________________________
82 TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
85 // Make linear binning
86 // the user has to delete the array afterwards!!!
93 TVectorD *binLim=new TVectorD(nbinsX+1);
96 Double_t binWidth=(last-first)/nbinsX;
97 for (Int_t i=0; i<nbinsX+1; ++i){
98 (*binLim)[i]=first+binWidth*(Double_t)i;
103 //_____________________________________________________________________________
104 TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
107 // Make arbitrary binning, bins separated by a ','
109 TString limits(bins);
110 if (limits.IsNull()){
111 Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
115 TObjArray *arr=limits.Tokenize(",");
116 Int_t nLimits=arr->GetEntries();
118 Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
123 TVectorD *binLimits=new TVectorD(nLimits);
124 for (Int_t iLim=0; iLim<nLimits; ++iLim){
125 (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
132 //_____________________________________________________________________________
133 void AliDielectronHelper::GetMaxPtAndPhi(const AliVEvent *ev, Double_t &ptMax, Double_t &phiOfptMax){
135 // find the highest pt and its phi in the event
137 for(Int_t itrk=0; itrk<ev->GetNumberOfTracks(); itrk++) {
138 AliVParticle *part= ev->GetTrack(itrk);
139 if(part && part->Pt() > ptMax) {
141 phiOfptMax = part->Phi();
148 //_____________________________________________________________________________
149 Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
150 // determination of Nch
151 if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
153 AliStack *stack = ((AliMCEvent*)ev)->Stack();
155 if (!stack) return -1;
157 Int_t nParticles = stack->GetNtrack();
161 for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
162 if (!stack->IsPhysicalPrimary(iMc)) continue;
164 TParticle* particle = stack->Particle(iMc);
165 if (!particle) continue;
166 if (particle->GetPDG()->Charge() == 0) continue;
168 Float_t eta = particle->Eta();
169 if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
176 //_____________________________________________________________________________
177 Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){
178 // Compute the collision multiplicity based on AOD or ESD tracklets
179 // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
183 Int_t nTracklets = 0;
186 if (ev->IsA() == AliAODEvent::Class()) {
187 AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
188 if (!tracklets) return -1;
189 nTracklets = tracklets->GetNumberOfTracklets();
190 for (Int_t nn = 0; nn < nTracklets; nn++) {
191 Double_t theta = tracklets->GetTheta(nn);
192 Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
193 if (TMath::Abs(eta) < etaRange) nAcc++;
195 } else if (ev->IsA() == AliESDEvent::Class()) {
196 nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
197 for (Int_t nn = 0; nn < nTracklets; nn++) {
198 Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
199 if (TMath::Abs(eta) < etaRange) nAcc++;
208 //________________________________________________________________
209 Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
211 // Correct the number of accepted tracklets based on the period and event vertex
214 Int_t runNo = event->GetRunNumber();
216 Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
217 Double_t refMult = 0.0; // reference multiplicity
219 if(runNo>114930 && runNo<117223) period = 0;
220 if(runNo>119158 && runNo<120830) period = 1;
221 if(runNo>122373 && runNo<126438) period = 2;
222 if(runNo>127711 && runNo<130841) period = 3;
223 if(period<0 || period>3) return uncorrectedNacc;
225 if(type<0 || type>8) return uncorrectedNacc;
226 if(type == 0) refMult = 5.0; // SPD tracklets in |eta|<0.5
227 if(type == 1) refMult = 9.5; // SPD tracklets in |eta|<1.0
228 if(type == 2) refMult = 13.0; // SPD tracklets in |eta|<1.6
229 if(type == 3) refMult = 6.0; // ITSTPC+ in |eta|<0.5
230 if(type == 4) refMult = 12.0; // ITSTPC+ in |eta|<1.0
231 if(type == 5) refMult = 16.0; // ITSTPC+ in |eta|<1.6
232 if(type == 6) refMult = 6.0; // ITSSA+ in |eta|<0.5
233 if(type == 7) refMult = 12.0; // ITSSA+ in |eta|<1.0
234 if(type == 8) refMult = 15.0; // ITSSA+ in |eta|<1.6
236 if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;
238 TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
239 if(!estimatorAvg) return uncorrectedNacc;
241 Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
243 Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
245 Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
247 return correctedNacc;
250 //_____________________________________________________________________________
251 Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
252 // put a robust Nacc definition here
256 AliDielectronVarCuts varCuts;
257 varCuts.AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0);
258 varCuts.AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0);
259 varCuts.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
260 varCuts.AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0);
261 varCuts.AddCut(AliDielectronVarManager::kNclsTPC, 70.0, 160.0);
262 varCuts.AddCut(AliDielectronVarManager::kKinkIndex0, -0.5, 0.5); //noKinks
264 AliDielectronTrackCuts trkCuts;
265 trkCuts.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
266 trkCuts.SetRequireITSRefit(kTRUE);
267 trkCuts.SetRequireTPCRefit(kTRUE);
269 Int_t nRecoTracks = ev->GetNumberOfTracks();
272 for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
273 AliVTrack *track = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
274 if (!track) continue;
275 if (!trkCuts.IsSelected(track)) continue;
276 if (!varCuts.IsSelected(track)) continue;
283 //_____________________________________________________________________________
284 Double_t AliDielectronHelper::GetITSTPCMatchEff(const AliVEvent *ev){
285 // recalulate the its-tpc matching efficiecy
289 AliDielectronVarCuts varCutsTPC;
290 varCutsTPC.AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0);
291 varCutsTPC.AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0);
292 varCutsTPC.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
293 varCutsTPC.AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0);
294 varCutsTPC.AddCut(AliDielectronVarManager::kNclsTPC, 50.0, 160.0);
295 AliDielectronTrackCuts trkCutsTPC;
296 trkCutsTPC.SetRequireTPCRefit(kTRUE);
298 AliDielectronVarCuts varCutsITS;
299 varCutsITS.AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
300 AliDielectronTrackCuts trkCutsITS;
301 trkCutsITS.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
302 trkCutsITS.SetRequireITSRefit(kTRUE);
305 Int_t nRecoTracks = ev->GetNumberOfTracks();
306 Double_t nTPC = 0, nITS = 0;
308 for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
309 AliVTrack *track = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
310 if (!track) continue;
312 if(!trkCutsITS.IsSelected(track)) continue;
313 if(!varCutsITS.IsSelected(track)) continue;
316 if(!trkCutsTPC.IsSelected(track)) continue;
317 if(!varCutsTPC.IsSelected(track)) continue;
322 // printf(" tracks TPC %.3e ITS %.3e = %.5f \n",nTPC,nITS,(nITS>0. ? nTPC/nITS : -1));
323 return (nITS>0. ? nTPC/nITS : -1);
327 //_____________________________________________________________________________
328 void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
329 // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
330 if (!kfParticle) return;
336 dx = ev->GetPrimaryVertex()->GetX()-0.;
337 dy = ev->GetPrimaryVertex()->GetY()-0.;
338 dz = ev->GetPrimaryVertex()->GetZ()-0.;
341 kfParticle->X() = kfParticle->GetX() - dx;
342 kfParticle->Y() = kfParticle->GetY() - dy;
343 kfParticle->Z() = kfParticle->GetZ() - dz;
346 // Rotate the kf particle
347 Double_t c = cos(angle);
348 Double_t s = sin(angle);
351 for( Int_t i=0; i<8; i++ ){
352 for( Int_t j=0; j<8; j++){
356 for( int i=0; i<8; i++ ){
359 mA[0][0] = c; mA[0][1] = s;
360 mA[1][0] = -s; mA[1][1] = c;
361 mA[3][3] = c; mA[3][4] = s;
362 mA[4][3] = -s; mA[4][4] = c;
367 for( Int_t i=0; i<8; i++ ){
369 for( Int_t k=0; k<8; k++){
370 mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
374 for( Int_t i=0; i<8; i++){
375 kfParticle->Parameter(i) = mAp[i];
378 for( Int_t i=0; i<8; i++ ){
379 for( Int_t j=0; j<8; j++ ){
381 for( Int_t k=0; k<8; k++ ){
382 mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
387 for( Int_t i=0; i<8; i++ ){
388 for( Int_t j=0; j<=i; j++ ){
390 for( Int_t k=0; k<8; k++){
391 xx+= mAC[i][k]*mA[j][k];
393 kfParticle->Covariance(i,j) = xx;
397 kfParticle->X() = kfParticle->GetX() + dx;
398 kfParticle->Y() = kfParticle->GetY() + dy;
399 kfParticle->Z() = kfParticle->GetZ() + dz;
403 //_____________________________________________________________________________
404 Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
406 // counting number of mother particles generated in given eta range and 2 particle decay
407 if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
409 AliStack *stack = ((AliMCEvent*)ev)->Stack();
411 if (!stack) return -1;
413 Int_t nParticles = stack->GetNtrack();
417 for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
419 TParticle* particle = stack->Particle(iMc);
420 if (!particle) continue;
421 if (particle->GetPdgCode() != pdgMother) continue;
422 if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
424 if (particle->GetNDaughters() != 2) continue;
426 if (particle->GetFirstDaughter()>=nParticles ||
427 particle->GetFirstDaughter()<0 ) continue;
429 TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
430 if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter) continue;
431 if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
434 if (particle->GetLastDaughter()>=nParticles ||
435 particle->GetLastDaughter()<0 ) continue;
437 TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
438 if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter) continue;
439 if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
443 if(particle->IsPrimary() != prim) continue;