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 //
  23
  24
  25
  26
  27 #include <TError.h>
  28 #include <TMath.h>
  29 #include <TObjString.h>
  30 #include <TObjArray.h>
  31 #include <TVectorD.h>
  32
  33 #include <AliVEvent.h>
  34 #include <AliKFParticle.h>
  35
  36 #include "AliDielectronHelper.h"
  37
  38 //_____________________________________________________________________________
  39 TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
  40 {
  41   //
  42   // Make logarithmic binning
  43   // the user has to delete the array afterwards!!!
  44   //
  45
  46   //check limits
  47   if (xmin<1e-20 || xmax<1e-20){
  48     Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
  49     return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax);
  50   }
  51   if (xmax<xmin){
  52     Double_t tmp=xmin;
  53     xmin=xmax;
  54     xmax=tmp;
  55   }
  56   TVectorD *binLim=new TVectorD(nbinsX+1);
  57   Double_t first=xmin;
  58   Double_t last=xmax;
  59   Double_t expMax=TMath::Log(last/first);
  60   for (Int_t i=0; i<nbinsX+1; ++i){
  61     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
  62   }
  63   return binLim;
  64 }
  65
  66 //_____________________________________________________________________________
  67 TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
  68 {
  69   //
  70   // Make linear binning
  71   // the user has to delete the array afterwards!!!
  72   //
  73   if (xmax<xmin){
  74     Double_t tmp=xmin;
  75     xmin=xmax;
  76     xmax=tmp;
  77   }
  78   TVectorD *binLim=new TVectorD(nbinsX+1);
  79   Double_t first=xmin;
  80   Double_t last=xmax;
  81   Double_t binWidth=(last-first)/nbinsX;
  82   for (Int_t i=0; i<nbinsX+1; ++i){
  83     (*binLim)[i]=first+binWidth*(Double_t)i;
  84   }
  85   return binLim;
  86 }
  87
  88 //_____________________________________________________________________________
  89 TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
  90 {
  91   //
  92   // Make arbitrary binning, bins separated by a ','
  93   //
  94   TString limits(bins);
  95   if (limits.IsNull()){
  96     Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
  97     return 0x0;
  98   }
  99
 100   TObjArray *arr=limits.Tokenize(",");
 101   Int_t nLimits=arr->GetEntries();
 102   if (nLimits<2){
 103     Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
 104     delete arr;
 105     return 0x0;
 106   }
 107
 108   TVectorD *binLimits=new TVectorD(nLimits);
 109   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 110     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 111   }
 112
 113   delete arr;
 114   return binLimits;
 115 }
 116
 117 //_____________________________________________________________________________
 118 void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
 119   // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
 120   if (!kfParticle) return;
 121   Double_t dx = 0.;
 122   Double_t dy = 0.;
 123   Double_t dz = 0.;
 124
 125   if (ev){
 126     dx = ev->GetPrimaryVertex()->GetX()-0.;
 127     dy = ev->GetPrimaryVertex()->GetY()-0.;
 128     dz = ev->GetPrimaryVertex()->GetZ()-0.;
 129   }
 130
 131   kfParticle->X() = kfParticle->GetX() - dx;
 132   kfParticle->Y() = kfParticle->GetY() - dy;
 133   kfParticle->Z() = kfParticle->GetZ() - dz;
 134
 135
 136   // Rotate the kf particle
 137   Double_t c = cos(angle);
 138   Double_t s = sin(angle);
 139
 140   Double_t mA[8][ 8];
 141   for( Int_t i=0; i<8; i++ ){
 142     for( Int_t j=0; j<8; j++){
 143       mA[i][j] = 0;
 144     }
 145   }
 146   for( int i=0; i<8; i++ ){
 147     mA[i][i] = 1;
 148   }
 149   mA[0][0] =  c;  mA[0][1] = s;
 150   mA[1][0] = -s;  mA[1][1] = c;
 151   mA[3][3] =  c;  mA[3][4] = s;
 152   mA[4][3] = -s;  mA[4][4] = c;
 153
 154   Double_t mAC[8][8];
 155   Double_t mAp[8];
 156
 157   for( Int_t i=0; i<8; i++ ){
 158     mAp[i] = 0;
 159     for( Int_t k=0; k<8; k++){
 160       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
 161     }
 162   }
 163
 164   for( Int_t i=0; i<8; i++){
 165     kfParticle->Parameter(i) = mAp[i];
 166   }
 167
 168   for( Int_t i=0; i<8; i++ ){
 169     for( Int_t j=0; j<8; j++ ){
 170       mAC[i][j] = 0;
 171       for( Int_t k=0; k<8; k++ ){
 172         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
 173       }
 174     }
 175   }
 176
 177   for( Int_t i=0; i<8; i++ ){
 178     for( Int_t j=0; j<=i; j++ ){
 179       Double_t xx = 0;
 180       for( Int_t k=0; k<8; k++){
 181         xx+= mAC[i][k]*mA[j][k];
 182       }
 183       kfParticle->Covariance(i,j) = xx;
 184     }
 185   }
 186
 187   kfParticle->X() = kfParticle->GetX() + dx;
 188   kfParticle->Y() = kfParticle->GetY() + dy;
 189   kfParticle->Z() = kfParticle->GetZ() + dz;
 190
 191 }
 192