Add fast merging option (Diego)
[u/mrichter/AliRoot.git] / PWG3 / 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
34 #include <AliVEvent.h>
35 #include <AliVParticle.h>
36 #include <AliKFParticle.h>
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>
44
45 #include "AliDielectronHelper.h"
46
47 //_____________________________________________________________________________
48 TVectorD* 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 //_____________________________________________________________________________
76 TVectorD* 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 //_____________________________________________________________________________
98 TVectorD* 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
126
127 //_____________________________________________________________________________
128 Int_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 //_____________________________________________________________________________
156 Int_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 //_____________________________________________________________________________
187 Int_t AliDielectronHelper::GetNacc(const AliVEvent */*ev*/){
188   // put a robust Nacc definition here
189
190   return -1;
191 /*  
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;
215 */
216 }
217
218
219 //_____________________________________________________________________________
220 void 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