7418b8002151059af44ebf8fa4437409472dd077
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronHelper.cxx
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> 
22 //   Frederick Kramer <Frederick.Kramer@cern.ch> 
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 #include <TF1.h>
34 #include <TRandom.h>
35 #include <TProfile.h>
36
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>
46 #include <AliStack.h>
47
48 #include "AliDielectronVarCuts.h"
49 #include "AliDielectronTrackCuts.h"
50 #include "AliDielectronVarManager.h"
51 #include "AliDielectronHelper.h"
52
53 //_____________________________________________________________________________
54 TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
55 {
56   //
57   // Make logarithmic binning
58   // the user has to delete the array afterwards!!!
59   //
60   
61   //check limits
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);
65   }
66   if (xmax<xmin){
67     Double_t tmp=xmin;
68     xmin=xmax;
69     xmax=tmp;
70   }
71   TVectorD *binLim=new TVectorD(nbinsX+1);
72   Double_t first=xmin;
73   Double_t last=xmax;
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);
77   }
78   return binLim;
79 }
80
81 //_____________________________________________________________________________
82 TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
83 {
84   //
85   // Make linear binning
86   // the user has to delete the array afterwards!!!
87   //
88   if (xmax<xmin){
89     Double_t tmp=xmin;
90     xmin=xmax;
91     xmax=tmp;
92   }
93   TVectorD *binLim=new TVectorD(nbinsX+1);
94   Double_t first=xmin;
95   Double_t last=xmax;
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;
99   }
100   return binLim;
101 }
102
103 //_____________________________________________________________________________
104 TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
105 {
106   //
107   // Make arbitrary binning, bins separated by a ','
108   //
109   TString limits(bins);
110   if (limits.IsNull()){
111     Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
112     return 0x0;
113   }
114   
115   TObjArray *arr=limits.Tokenize(",");
116   Int_t nLimits=arr->GetEntries();
117   if (nLimits<2){
118     Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
119     delete arr;
120     return 0x0;
121   }
122   
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();
126   }
127   
128   delete arr;
129   return binLimits;
130 }
131
132
133 //_____________________________________________________________________________
134 Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
135   // determination of Nch
136   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
137
138   AliStack *stack = ((AliMCEvent*)ev)->Stack();
139
140   if (!stack) return -1;
141
142   Int_t nParticles = stack->GetNtrack();
143   Int_t nCh = 0;
144
145   // count..
146   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
147     if (!stack->IsPhysicalPrimary(iMc)) continue;
148
149     TParticle* particle = stack->Particle(iMc);
150     if (!particle) continue;
151     if (particle->GetPDG()->Charge() == 0) continue;
152
153     Float_t eta = particle->Eta();
154     if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
155   }
156
157   return nCh;
158 }
159
160
161 //_____________________________________________________________________________
162 Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){
163   // Compute the collision multiplicity based on AOD or ESD tracklets
164   // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
165
166   if (!ev) return -1;
167
168   Int_t nTracklets = 0;
169   Int_t nAcc = 0;
170   
171   if (ev->IsA() == AliAODEvent::Class()) {
172     AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
173     nTracklets = tracklets->GetNumberOfTracklets();
174     for (Int_t nn = 0; nn < nTracklets; nn++) {
175       Double_t theta = tracklets->GetTheta(nn);
176       Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
177       if (TMath::Abs(eta) < etaRange) nAcc++;
178     }
179   } else if (ev->IsA() == AliESDEvent::Class()) {
180     nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
181     for (Int_t nn = 0; nn < nTracklets; nn++) {
182       Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
183       if (TMath::Abs(eta) < etaRange) nAcc++;
184     }
185   } else return -1;
186
187   return nAcc;
188 }
189
190
191
192 //________________________________________________________________
193 Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
194   //
195   // Correct the number of accepted tracklets based on the period and event vertex
196   //
197
198   Int_t runNo = event->GetRunNumber();
199
200   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
201   Double_t refMult = 0.0;   // reference multiplicity
202   
203   if(runNo>114930 && runNo<117223) period = 0;
204   if(runNo>119158 && runNo<120830) period = 1;
205   if(runNo>122373 && runNo<126438) period = 2;
206   if(runNo>127711 && runNo<130841) period = 3;
207   if(period<0 || period>3) return uncorrectedNacc;
208
209   if(type<0 || type>8) return uncorrectedNacc;
210   if(type == 0) refMult = 5.0;         // SPD tracklets in |eta|<0.5 
211   if(type == 1) refMult = 9.5;         // SPD tracklets in |eta|<1.0
212   if(type == 2) refMult = 13.0;        // SPD tracklets in |eta|<1.6
213   if(type == 3) refMult = 6.0;         // ITSTPC+ in |eta|<0.5
214   if(type == 4) refMult = 12.0;        // ITSTPC+ in |eta|<1.0
215   if(type == 5) refMult = 16.0;        // ITSTPC+ in |eta|<1.6
216   if(type == 6) refMult = 6.0;         // ITSSA+ in |eta|<0.5
217   if(type == 7) refMult = 12.0;        // ITSSA+ in |eta|<1.0
218   if(type == 8) refMult = 15.0;        // ITSSA+ in |eta|<1.6
219
220   if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;
221
222   TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
223   if(!estimatorAvg) return uncorrectedNacc;
224
225   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
226
227   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
228
229   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
230
231   return correctedNacc;
232 }
233
234 //_____________________________________________________________________________
235 Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
236   // put a robust Nacc definition here
237
238   if (!ev) return -1;
239   
240   AliDielectronVarCuts *varCuts = new AliDielectronVarCuts("VarCuts","VarCuts");
241   varCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
242   varCuts->AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
243   varCuts->AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
244   varCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
245   varCuts->AddCut(AliDielectronVarManager::kNclsTPC,     70.0, 160.0);
246   varCuts->AddCut(AliDielectronVarManager::kKinkIndex0,  -0.5,   0.5);   //noKinks
247     
248   AliDielectronTrackCuts *trkCuts = new AliDielectronTrackCuts("TrkCuts","TrkCuts");
249   trkCuts->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
250   trkCuts->SetRequireITSRefit(kTRUE);
251   trkCuts->SetRequireTPCRefit(kTRUE);
252
253   Int_t nRecoTracks = ev->GetNumberOfTracks();
254   Int_t nAcc = 0;
255   
256   for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
257     AliVParticle* candidate = ev->GetTrack(iTrack);
258     if (!candidate) continue;
259     AliVTrack *track        = static_cast<AliVTrack*>(candidate);
260     if (!track) continue;
261     if (varCuts->IsSelected(track) && trkCuts->IsSelected(track)) 
262       nAcc++;
263   }
264   
265   delete varCuts;
266   delete trkCuts;
267
268   return nAcc;
269 }
270
271 //_____________________________________________________________________________
272 void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
273   // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
274   if (!kfParticle) return;
275   Double_t dx = 0.;
276   Double_t dy = 0.;
277   Double_t dz = 0.;
278
279   if (ev){
280     dx = ev->GetPrimaryVertex()->GetX()-0.;
281     dy = ev->GetPrimaryVertex()->GetY()-0.;
282     dz = ev->GetPrimaryVertex()->GetZ()-0.;
283   }
284   
285   kfParticle->X() = kfParticle->GetX() - dx;
286   kfParticle->Y() = kfParticle->GetY() - dy;
287   kfParticle->Z() = kfParticle->GetZ() - dz;
288   
289   
290   // Rotate the kf particle
291   Double_t c = cos(angle);
292   Double_t s = sin(angle);
293   
294   Double_t mA[8][ 8];
295   for( Int_t i=0; i<8; i++ ){
296     for( Int_t j=0; j<8; j++){
297       mA[i][j] = 0;
298     }
299   }
300   for( int i=0; i<8; i++ ){
301     mA[i][i] = 1;
302   }
303   mA[0][0] =  c;  mA[0][1] = s;
304   mA[1][0] = -s;  mA[1][1] = c;
305   mA[3][3] =  c;  mA[3][4] = s;
306   mA[4][3] = -s;  mA[4][4] = c;
307   
308   Double_t mAC[8][8];
309   Double_t mAp[8];
310   
311   for( Int_t i=0; i<8; i++ ){
312     mAp[i] = 0;
313     for( Int_t k=0; k<8; k++){
314       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
315     }
316   }
317   
318   for( Int_t i=0; i<8; i++){
319     kfParticle->Parameter(i) = mAp[i];
320   }
321   
322   for( Int_t i=0; i<8; i++ ){
323     for( Int_t j=0; j<8; j++ ){
324       mAC[i][j] = 0;
325       for( Int_t k=0; k<8; k++ ){
326         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
327       }
328     }
329   }
330   
331   for( Int_t i=0; i<8; i++ ){
332     for( Int_t j=0; j<=i; j++ ){
333       Double_t xx = 0;
334       for( Int_t k=0; k<8; k++){
335         xx+= mAC[i][k]*mA[j][k];
336       }
337       kfParticle->Covariance(i,j) = xx;
338     }
339   }
340   
341   kfParticle->X() = kfParticle->GetX() + dx;
342   kfParticle->Y() = kfParticle->GetY() + dy;
343   kfParticle->Z() = kfParticle->GetZ() + dz;
344   
345 }
346
347 //_____________________________________________________________________________
348 Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
349   // counting number of mother particles generated in given eta range and 2 particle decay
350   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
351   
352   AliStack *stack = ((AliMCEvent*)ev)->Stack();
353   
354   if (!stack) return -1;
355   
356   Int_t nParticles = stack->GetNtrack();
357   Int_t nMothers   = 0;
358   
359   // count..
360   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
361     
362     TParticle* particle = stack->Particle(iMc);
363     if (!particle) continue;
364     if (particle->GetPdgCode() != pdgMother)               continue;
365     if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
366
367     if (particle->GetNDaughters() != 2)                 continue;
368     // 1st daugther
369     if (particle->GetFirstDaughter()>=nParticles ||
370         particle->GetFirstDaughter()<0             ) continue;
371     
372     TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
373     if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter)     continue;
374     if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
375     
376     // 2nd daughter
377     if (particle->GetLastDaughter()>=nParticles ||
378         particle->GetLastDaughter()<0             ) continue;
379
380     TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
381     if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter)     continue;
382     if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
383     
384     // primary
385     if (prim != -1) {
386       if(particle->IsPrimary() != prim) continue;
387     }
388     nMothers++;
389   }
390   return nMothers;
391 }
392
393