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 //_____________________________________________________________________________
  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 void AliDielectronHelper::GetMaxPtAndPhi(const AliVEvent *ev, Double_t &ptMax, Double_t &phiOfptMax){
 134   //
 135   // find the highest pt and its phi in the event
 136   //
 137   for(Int_t itrk=0; itrk<ev->GetNumberOfTracks(); itrk++) {
 138     AliVParticle *part= ev->GetTrack(itrk);
 139     if(part && part->Pt() > ptMax) {
 140       ptMax      = part->Pt();
 141       phiOfptMax = part->Phi();
 142     }
 143   }
 144
 145
 146 }
 147
 148 //_____________________________________________________________________________
 149 Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
 150   // determination of Nch
 151   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
 152
 153   AliStack *stack = ((AliMCEvent*)ev)->Stack();
 154
 155   if (!stack) return -1;
 156
 157   Int_t nParticles = stack->GetNtrack();
 158   Int_t nCh = 0;
 159
 160   // count..
 161   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
 162     if (!stack->IsPhysicalPrimary(iMc)) continue;
 163
 164     TParticle* particle = stack->Particle(iMc);
 165     if (!particle) continue;
 166     if (particle->GetPDG()->Charge() == 0) continue;
 167
 168     Float_t eta = particle->Eta();
 169     if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
 170   }
 171
 172   return nCh;
 173 }
 174
 175
 176 //_____________________________________________________________________________
 177 Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){
 178   // Compute the collision multiplicity based on AOD or ESD tracklets
 179   // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
 180
 181   if (!ev) return -1;
 182
 183   Int_t nTracklets = 0;
 184   Int_t nAcc = 0;
 185
 186   if (ev->IsA() == AliAODEvent::Class()) {
 187     AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
 188     if (!tracklets) return -1;
 189     nTracklets = tracklets->GetNumberOfTracklets();
 190     for (Int_t nn = 0; nn < nTracklets; nn++) {
 191       Double_t theta = tracklets->GetTheta(nn);
 192       Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
 193       if (TMath::Abs(eta) < etaRange) nAcc++;
 194     }
 195   } else if (ev->IsA() == AliESDEvent::Class()) {
 196     nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
 197     for (Int_t nn = 0; nn < nTracklets; nn++) {
 198       Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
 199       if (TMath::Abs(eta) < etaRange) nAcc++;
 200     }
 201   } else return -1;
 202
 203   return nAcc;
 204 }
 205
 206
 207
 208 //________________________________________________________________
 209 Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
 210   //
 211   // Correct the number of accepted tracklets based on the period and event vertex
 212   //
 213
 214   Int_t runNo = event->GetRunNumber();
 215
 216   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 217   Double_t refMult = 0.0;   // reference multiplicity
 218
 219   if(runNo>114930 && runNo<117223) period = 0;
 220   if(runNo>119158 && runNo<120830) period = 1;
 221   if(runNo>122373 && runNo<126438) period = 2;
 222   if(runNo>127711 && runNo<130841) period = 3;
 223   if(period<0 || period>3) return uncorrectedNacc;
 224
 225   if(type<0 || type>8) return uncorrectedNacc;
 226   if(type == 0) refMult = 5.0;         // SPD tracklets in |eta|<0.5
 227   if(type == 1) refMult = 9.5;         // SPD tracklets in |eta|<1.0
 228   if(type == 2) refMult = 13.0;        // SPD tracklets in |eta|<1.6
 229   if(type == 3) refMult = 6.0;         // ITSTPC+ in |eta|<0.5
 230   if(type == 4) refMult = 12.0;        // ITSTPC+ in |eta|<1.0
 231   if(type == 5) refMult = 16.0;        // ITSTPC+ in |eta|<1.6
 232   if(type == 6) refMult = 6.0;         // ITSSA+ in |eta|<0.5
 233   if(type == 7) refMult = 12.0;        // ITSSA+ in |eta|<1.0
 234   if(type == 8) refMult = 15.0;        // ITSSA+ in |eta|<1.6
 235
 236   if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;
 237
 238   TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
 239   if(!estimatorAvg) return uncorrectedNacc;
 240
 241   Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
 242
 243   Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
 244
 245   Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
 246
 247   return correctedNacc;
 248 }
 249
 250 //_____________________________________________________________________________
 251 Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
 252   // put a robust Nacc definition here
 253
 254   if (!ev) return -1;
 255
 256   AliDielectronVarCuts varCuts;
 257   varCuts.AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
 258   varCuts.AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
 259   varCuts.AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
 260   varCuts.AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
 261   varCuts.AddCut(AliDielectronVarManager::kNclsTPC,     70.0, 160.0);
 262   varCuts.AddCut(AliDielectronVarManager::kKinkIndex0,  -0.5,   0.5);   //noKinks
 263
 264   AliDielectronTrackCuts trkCuts;
 265   trkCuts.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
 266   trkCuts.SetRequireITSRefit(kTRUE);
 267   trkCuts.SetRequireTPCRefit(kTRUE);
 268
 269   Int_t nRecoTracks = ev->GetNumberOfTracks();
 270   Int_t nAcc = 0;
 271
 272   for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
 273     AliVTrack *track        = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
 274     if (!track) continue;
 275     if (!trkCuts.IsSelected(track)) continue;
 276     if (!varCuts.IsSelected(track)) continue;
 277     nAcc++;
 278   }
 279
 280   return nAcc;
 281 }
 282
 283 //_____________________________________________________________________________
 284 Double_t AliDielectronHelper::GetITSTPCMatchEff(const AliVEvent *ev){
 285   // recalulate the its-tpc matching efficiecy
 286
 287   if (!ev) return -1;
 288
 289   AliDielectronVarCuts varCutsTPC;
 290   varCutsTPC.AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
 291   varCutsTPC.AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
 292   varCutsTPC.AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
 293   varCutsTPC.AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
 294   varCutsTPC.AddCut(AliDielectronVarManager::kNclsTPC,     50.0, 160.0);
 295   AliDielectronTrackCuts trkCutsTPC;
 296   trkCutsTPC.SetRequireTPCRefit(kTRUE);
 297
 298   AliDielectronVarCuts varCutsITS;
 299   varCutsITS.AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
 300   AliDielectronTrackCuts trkCutsITS;
 301   trkCutsITS.SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
 302   trkCutsITS.SetRequireITSRefit(kTRUE);
 303
 304
 305   Int_t nRecoTracks = ev->GetNumberOfTracks();
 306   Double_t nTPC = 0, nITS = 0;
 307
 308   for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
 309     AliVTrack *track        = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
 310     if (!track) continue;
 311
 312     if(!trkCutsITS.IsSelected(track)) continue;
 313     if(!varCutsITS.IsSelected(track)) continue;
 314     nITS+=1.;
 315
 316     if(!trkCutsTPC.IsSelected(track)) continue;
 317     if(!varCutsTPC.IsSelected(track)) continue;
 318     nTPC+=1.;
 319
 320   }
 321
 322   //  printf(" tracks TPC %.3e ITS %.3e = %.5f \n",nTPC,nITS,(nITS>0. ? nTPC/nITS : -1));
 323   return (nITS>0. ? nTPC/nITS : -1);
 324 }
 325
 326
 327 //_____________________________________________________________________________
 328 void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
 329   // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
 330   if (!kfParticle) return;
 331   Double_t dx = 0.;
 332   Double_t dy = 0.;
 333   Double_t dz = 0.;
 334
 335   if (ev){
 336     dx = ev->GetPrimaryVertex()->GetX()-0.;
 337     dy = ev->GetPrimaryVertex()->GetY()-0.;
 338     dz = ev->GetPrimaryVertex()->GetZ()-0.;
 339   }
 340
 341   kfParticle->X() = kfParticle->GetX() - dx;
 342   kfParticle->Y() = kfParticle->GetY() - dy;
 343   kfParticle->Z() = kfParticle->GetZ() - dz;
 344
 345
 346   // Rotate the kf particle
 347   Double_t c = cos(angle);
 348   Double_t s = sin(angle);
 349
 350   Double_t mA[8][ 8];
 351   for( Int_t i=0; i<8; i++ ){
 352     for( Int_t j=0; j<8; j++){
 353       mA[i][j] = 0;
 354     }
 355   }
 356   for( int i=0; i<8; i++ ){
 357     mA[i][i] = 1;
 358   }
 359   mA[0][0] =  c;  mA[0][1] = s;
 360   mA[1][0] = -s;  mA[1][1] = c;
 361   mA[3][3] =  c;  mA[3][4] = s;
 362   mA[4][3] = -s;  mA[4][4] = c;
 363
 364   Double_t mAC[8][8];
 365   Double_t mAp[8];
 366
 367   for( Int_t i=0; i<8; i++ ){
 368     mAp[i] = 0;
 369     for( Int_t k=0; k<8; k++){
 370       mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
 371     }
 372   }
 373
 374   for( Int_t i=0; i<8; i++){
 375     kfParticle->Parameter(i) = mAp[i];
 376   }
 377
 378   for( Int_t i=0; i<8; i++ ){
 379     for( Int_t j=0; j<8; j++ ){
 380       mAC[i][j] = 0;
 381       for( Int_t k=0; k<8; k++ ){
 382         mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
 383       }
 384     }
 385   }
 386
 387   for( Int_t i=0; i<8; i++ ){
 388     for( Int_t j=0; j<=i; j++ ){
 389       Double_t xx = 0;
 390       for( Int_t k=0; k<8; k++){
 391         xx+= mAC[i][k]*mA[j][k];
 392       }
 393       kfParticle->Covariance(i,j) = xx;
 394     }
 395   }
 396
 397   kfParticle->X() = kfParticle->GetX() + dx;
 398   kfParticle->Y() = kfParticle->GetY() + dy;
 399   kfParticle->Z() = kfParticle->GetZ() + dz;
 400
 401 }
 402
 403 //_____________________________________________________________________________
 404 Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
 405   // TODO: add AODs
 406   // counting number of mother particles generated in given eta range and 2 particle decay
 407   if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
 408
 409   AliStack *stack = ((AliMCEvent*)ev)->Stack();
 410
 411   if (!stack) return -1;
 412
 413   Int_t nParticles = stack->GetNtrack();
 414   Int_t nMothers   = 0;
 415
 416   // count..
 417   for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
 418
 419     TParticle* particle = stack->Particle(iMc);
 420     if (!particle) continue;
 421     if (particle->GetPdgCode() != pdgMother)               continue;
 422     if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
 423
 424     if (particle->GetNDaughters() != 2)                 continue;
 425     // 1st daugther
 426     if (particle->GetFirstDaughter()>=nParticles ||
 427         particle->GetFirstDaughter()<0             ) continue;
 428
 429     TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
 430     if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter)     continue;
 431     if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
 432
 433     // 2nd daughter
 434     if (particle->GetLastDaughter()>=nParticles ||
 435         particle->GetLastDaughter()<0             ) continue;
 436
 437     TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
 438     if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter)     continue;
 439     if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
 440
 441     // primary
 442     if (prim != -1) {
 443       if(particle->IsPrimary() != prim) continue;
 444     }
 445     nMothers++;
 446   }
 447   return nMothers;
 448 }
 449
 450