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