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     if (!tracklets) return -1;
 174     nTracklets = tracklets->GetNumberOfTracklets();
 175     for (Int_t nn = 0; nn < nTracklets; nn++) {
 176       Double_t theta = tracklets->GetTheta(nn);
 177       Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
 178       if (TMath::Abs(eta) < etaRange) nAcc++;
 179     }
 180   } else if (ev->IsA() == AliESDEvent::Class()) {
 181     nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
 182     for (Int_t nn = 0; nn < nTracklets; nn++) {
 183       Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
 184       if (TMath::Abs(eta) < etaRange) nAcc++;
 185     }
 186   } else return -1;
 187
 188   return nAcc;
 189 }
 190
 191
 192
 193 //________________________________________________________________
 194 Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
 195   //
 196   // Correct the number of accepted tracklets based on the period and event vertex
 197   //
 198
 199   Int_t runNo = event->GetRunNumber();
 200
 201   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 202   Double_t refMult = 0.0;   // reference multiplicity
 203
 204   if(runNo>114930 && runNo<117223) period = 0;
 205   if(runNo>119158 && runNo<120830) period = 1;
 206   if(runNo>122373 && runNo<126438) period = 2;
 207   if(runNo>127711 && runNo<130841) period = 3;
 208   if(period<0 || period>3) return uncorrectedNacc;
 209
 210   if(type<0 || type>8) return uncorrectedNacc;
 211   if(type == 0) refMult = 5.0;         // SPD tracklets in |eta|<0.5
 212   if(type == 1) refMult = 9.5;         // SPD tracklets in |eta|<1.0
 213   if(type == 2) refMult = 13.0;        // SPD tracklets in |eta|<1.6
 214   if(type == 3) refMult = 6.0;         // ITSTPC+ in |eta|<0.5
 215   if(type == 4) refMult = 12.0;        // ITSTPC+ in |eta|<1.0
 216   if(type == 5) refMult = 16.0;        // ITSTPC+ in |eta|<1.6
 217   if(type == 6) refMult = 6.0;         // ITSSA+ in |eta|<0.5
 218   if(type == 7) refMult = 12.0;        // ITSSA+ in |eta|<1.0
 219   if(type == 8) refMult = 15.0;        // ITSSA+ in |eta|<1.6
 220
 221   if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;
 222
 223   TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
 224   if(!estimatorAvg) return uncorrectedNacc;
 225
 226   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 227
 228   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 229
 230   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 231
 232   return correctedNacc;
 233 }
 234
 235 //_____________________________________________________________________________
 236 Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
 237   // put a robust Nacc definition here
 238
 239   if (!ev) return -1;
 240
 241   AliDielectronVarCuts *varCuts = new AliDielectronVarCuts("VarCuts","VarCuts");
 242   varCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
 243   varCuts->AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
 244   varCuts->AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
 245   varCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
 246   varCuts->AddCut(AliDielectronVarManager::kNclsTPC,     70.0, 160.0);
 247   varCuts->AddCut(AliDielectronVarManager::kKinkIndex0,  -0.5,   0.5);   //noKinks
 248
 249   AliDielectronTrackCuts *trkCuts = new AliDielectronTrackCuts("TrkCuts","TrkCuts");
 250   trkCuts->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
 251   trkCuts->SetRequireITSRefit(kTRUE);
 252   trkCuts->SetRequireTPCRefit(kTRUE);
 253
 254   Int_t nRecoTracks = ev->GetNumberOfTracks();
 255   Int_t nAcc = 0;
 256
 257   for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
 258     AliVTrack *track        = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
 259     if (!track) continue;
 260     if (varCuts->IsSelected(track) && trkCuts->IsSelected(track))
 261       nAcc++;
 262   }
 263
 264   delete varCuts;
 265   delete trkCuts;
 266
 267   return nAcc;
 268 }
 269
 270 //_____________________________________________________________________________
 271 void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
 272   // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
 273   if (!kfParticle) return;
 274   Double_t dx = 0.;
 275   Double_t dy = 0.;
 276   Double_t dz = 0.;
 277
 278   if (ev){
 279     dx = ev->GetPrimaryVertex()->GetX()-0.;
 280     dy = ev->GetPrimaryVertex()->GetY()-0.;
 281     dz = ev->GetPrimaryVertex()->GetZ()-0.;
 282   }
 283
 284   kfParticle->X() = kfParticle->GetX() - dx;
 285   kfParticle->Y() = kfParticle->GetY() - dy;
 286   kfParticle->Z() = kfParticle->GetZ() - dz;
 287
 288
 289   // Rotate the kf particle
 290   Double_t c = cos(angle);
 291   Double_t s = sin(angle);
 292
 293   Double_t mA[8][ 8];
 294   for( Int_t i=0; i<8; i++ ){
 295     for( Int_t j=0; j<8; j++){
 296       mA[i][j] = 0;
 297     }
 298   }
 299   for( int i=0; i<8; i++ ){
 300     mA[i][i] = 1;
 301   }
 302   mA[0][0] =  c;  mA[0][1] = s;
 303   mA[1][0] = -s;  mA[1][1] = c;
 304   mA[3][3] =  c;  mA[3][4] = s;
 305   mA[4][3] = -s;  mA[4][4] = c;
 306
 307   Double_t mAC[8][8];
 308   Double_t mAp[8];
 309
 310   for( Int_t i=0; i<8; i++ ){
 311     mAp[i] = 0;
 312     for( Int_t k=0; k<8; k++){
 313       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
 314     }
 315   }
 316
 317   for( Int_t i=0; i<8; i++){
 318     kfParticle->Parameter(i) = mAp[i];
 319   }
 320
 321   for( Int_t i=0; i<8; i++ ){
 322     for( Int_t j=0; j<8; j++ ){
 323       mAC[i][j] = 0;
 324       for( Int_t k=0; k<8; k++ ){
 325         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
 326       }
 327     }
 328   }
 329
 330   for( Int_t i=0; i<8; i++ ){
 331     for( Int_t j=0; j<=i; j++ ){
 332       Double_t xx = 0;
 333       for( Int_t k=0; k<8; k++){
 334         xx+= mAC[i][k]*mA[j][k];
 335       }
 336       kfParticle->Covariance(i,j) = xx;
 337     }
 338   }
 339
 340   kfParticle->X() = kfParticle->GetX() + dx;
 341   kfParticle->Y() = kfParticle->GetY() + dy;
 342   kfParticle->Z() = kfParticle->GetZ() + dz;
 343
 344 }
 345
 346 //_____________________________________________________________________________
 347 Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
 348   // TODO: add AODs
 349   // counting number of mother particles generated in given eta range and 2 particle decay
 350   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
 351
 352   AliStack *stack = ((AliMCEvent*)ev)->Stack();
 353
 354   if (!stack) return -1;
 355
 356   Int_t nParticles = stack->GetNtrack();
 357   Int_t nMothers   = 0;
 358
 359   // count..
 360   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
 361
 362     TParticle* particle = stack->Particle(iMc);
 363     if (!particle) continue;
 364     if (particle->GetPdgCode() != pdgMother)               continue;
 365     if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
 366
 367     if (particle->GetNDaughters() != 2)                 continue;
 368     // 1st daugther
 369     if (particle->GetFirstDaughter()>=nParticles ||
 370         particle->GetFirstDaughter()<0             ) continue;
 371
 372     TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
 373     if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter)     continue;
 374     if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
 375
 376     // 2nd daughter
 377     if (particle->GetLastDaughter()>=nParticles ||
 378         particle->GetLastDaughter()<0             ) continue;
 379
 380     TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
 381     if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter)     continue;
 382     if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
 383
 384     // primary
 385     if (prim != -1) {
 386       if(particle->IsPrimary() != prim) continue;
 387     }
 388     nMothers++;
 389   }
 390   return nMothers;
 391 }
 392
 393