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