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
  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