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   // TODO: add AODs
 350   // counting number of mother particles generated in given eta range and 2 particle decay
 351   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
 352
 353   AliStack *stack = ((AliMCEvent*)ev)->Stack();
 354
 355   if (!stack) return -1;
 356
 357   Int_t nParticles = stack->GetNtrack();
 358   Int_t nMothers   = 0;
 359
 360   // count..
 361   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
 362
 363     TParticle* particle = stack->Particle(iMc);
 364     if (!particle) continue;
 365     if (particle->GetPdgCode() != pdgMother)               continue;
 366     if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
 367
 368     if (particle->GetNDaughters() != 2)                 continue;
 369     // 1st daugther
 370     if (particle->GetFirstDaughter()>=nParticles ||
 371         particle->GetFirstDaughter()<0             ) continue;
 372
 373     TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
 374     if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter)     continue;
 375     if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
 376
 377     // 2nd daughter
 378     if (particle->GetLastDaughter()>=nParticles ||
 379         particle->GetLastDaughter()<0             ) continue;
 380
 381     TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
 382     if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter)     continue;
 383     if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
 384
 385     // primary
 386     if (prim != -1) {
 387       if(particle->IsPrimary() != prim) continue;
 388     }
 389     nMothers++;
 390   }
 391   return nMothers;
 392 }
 393
 394