PWG0/AliPWG0Helper.cxx

   1 /* $Id$ */
   2
   3 #include <AliPWG0Helper.h>
   4
   5 #include <TParticle.h>
   6 #include <TParticlePDG.h>
   7 #include <TH1.h>
   8 #include <TH2.h>
   9 #include <TH3.h>
  10 #include <TList.h>
  11 #include <TTree.h>
  12 #include <TBranch.h>
  13 #include <TLeaf.h>
  14
  15 #include <AliHeader.h>
  16 #include <AliStack.h>
  17 #include <AliLog.h>
  18
  19 #include <AliLog.h>
  20 #include <AliESD.h>
  21 #include <AliESDEvent.h>
  22 #include <AliESDVertex.h>
  23
  24 #include <AliGenEventHeader.h>
  25 #include <AliGenPythiaEventHeader.h>
  26 #include <AliGenCocktailEventHeader.h>
  27 #include <AliGenDPMjetEventHeader.h>
  28
  29 //____________________________________________________________________
  30 ClassImp(AliPWG0Helper)
  31
  32 Int_t AliPWG0Helper::fgLastProcessType = -1;
  33
  34 //____________________________________________________________________
  35 Bool_t AliPWG0Helper::IsEventTriggered(const AliESD* aEsd, Trigger trigger)
  36 {
  37   // see function with ULong64_t argument
  38
  39   ULong64_t triggerMask = aEsd->GetTriggerMask();
  40   return IsEventTriggered(triggerMask, trigger);
  41 }
  42
  43 //____________________________________________________________________
  44 Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
  45 {
  46   // check if the event was triggered
  47   //
  48   // this function needs the branch fTriggerMask
  49
  50   // definitions from p-p.cfg
  51   ULong64_t spdFO = (1 << 14);
  52   ULong64_t v0left = (1 << 11);
  53   ULong64_t v0right = (1 << 12);
  54
  55   switch (trigger)
  56   {
  57     case kMB1:
  58     {
  59       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
  60         return kTRUE;
  61       break;
  62     }
  63     case kMB2:
  64     {
  65       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
  66         return kTRUE;
  67       break;
  68     }
  69     case kSPDFASTOR:
  70     {
  71       if (triggerMask & spdFO)
  72         return kTRUE;
  73       break;
  74     }
  75   }
  76
  77   return kFALSE;
  78 }
  79
  80 //____________________________________________________________________
  81 const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
  82 {
  83   // Get the vertex from the ESD and returns it if the vertex is valid
  84   //
  85   // Second argument decides which vertex is used (this selects
  86   // also the quality criteria that are applied)
  87
  88   const AliESDVertex* vertex = 0;
  89   Float_t requiredZResolution = -1;
  90   if (analysisMode == kSPD || analysisMode == kTPCITS)
  91   {
  92     vertex = aEsd->GetPrimaryVertexSPD();
  93     requiredZResolution = 0.1;
  94     if (debug)
  95       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
  96   }
  97   else if (analysisMode == kTPC)
  98   {
  99     vertex = aEsd->GetPrimaryVertexTPC();
 100     requiredZResolution = 0.6;
 101     if (debug)
 102       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
 103   }
 104   else
 105     Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode);
 106
 107   if (!vertex) {
 108     if (debug)
 109       Printf("AliPWG0Helper::GetVertex: No vertex found in ESD");
 110     return 0;
 111   }
 112
 113   // check Ncontributors
 114   if (vertex->GetNContributors() <= 0) {
 115     if (debug)
 116       Printf("AliPWG0Helper::GetVertex: NContributors() <= 0");
 117     return 0;
 118   }
 119
 120   // check resolution
 121   Double_t zRes = vertex->GetZRes();
 122
 123   if (zRes == 0 || zRes > requiredZResolution) {
 124     if (debug)
 125       Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
 126     return 0;
 127   }
 128
 129   if (debug)
 130     Printf("AliPWG0Helper::GetVertex: Returning valid vertex");
 131
 132   return vertex;
 133 }
 134
 135 //____________________________________________________________________
 136 Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
 137 {
 138   //
 139   // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
 140   // shall be counted as a primary particle
 141   //
 142   // This function or a equivalent should be available in some common place of AliRoot
 143   //
 144   // WARNING: Call this function only for particles that are among the particles from the event generator!
 145   // --> stack->Particle(id) with id < stack->GetNprimary()
 146
 147   // if the particle has a daughter primary, we do not want to count it
 148   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
 149   {
 150     if (adebug)
 151       printf("Dropping particle because it has a daughter among the primaries.\n");
 152     return kFALSE;
 153   }
 154
 155   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
 156
 157   // skip quarks and gluon
 158   if (pdgCode <= 10 || pdgCode == 21)
 159   {
 160     if (adebug)
 161       printf("Dropping particle because it is a quark or gluon.\n");
 162     return kFALSE;
 163   }
 164
 165   if (strcmp(aParticle->GetName(),"XXX") == 0)
 166   {
 167     Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
 168     return kFALSE;
 169   }
 170
 171   TParticlePDG* pdgPart = aParticle->GetPDG();
 172
 173   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
 174   {
 175     Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
 176     return kFALSE;
 177   }
 178
 179   if (pdgPart->Charge() == 0)
 180   {
 181     if (adebug)
 182       printf("Dropping particle because it is not charged.\n");
 183     return kFALSE;
 184   }
 185
 186   return kTRUE;
 187 }
 188
 189 //____________________________________________________________________
 190 void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save)
 191 {
 192   // create projections of 3d hists to all 2d combinations
 193   // the histograms are not returned, just use them from memory or use this to create them in a file
 194
 195   TH1* proj = hist->Project3D("yx");
 196   proj->SetXTitle(hist->GetXaxis()->GetTitle());
 197   proj->SetYTitle(hist->GetYaxis()->GetTitle());
 198   if (save)
 199     proj->Write();
 200
 201   proj = hist->Project3D("zx");
 202   proj->SetXTitle(hist->GetXaxis()->GetTitle());
 203   proj->SetYTitle(hist->GetZaxis()->GetTitle());
 204   if (save)
 205     proj->Write();
 206
 207   proj = hist->Project3D("zy");
 208   proj->SetXTitle(hist->GetYaxis()->GetTitle());
 209   proj->SetYTitle(hist->GetZaxis()->GetTitle());
 210   if (save)
 211     proj->Write();
 212 }
 213
 214 //____________________________________________________________________
 215 void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save)
 216 {
 217   // create projections of the 3d hists divides them
 218   // axis decides to which plane, if axis is 0 to all planes
 219   // the histograms are not returned, just use them from memory or use this to create them in a file
 220
 221   if (axis == 0)
 222   {
 223     CreateDividedProjections(hist, hist2, "yx", putErrors, save);
 224     CreateDividedProjections(hist, hist2, "zx", putErrors, save);
 225     CreateDividedProjections(hist, hist2, "zy", putErrors, save);
 226
 227     return;
 228   }
 229
 230   TH1* proj = hist->Project3D(axis);
 231
 232   if (strlen(axis) == 2)
 233   {
 234     proj->SetYTitle(GetAxisTitle(hist, axis[0]));
 235     proj->SetXTitle(GetAxisTitle(hist, axis[1]));
 236   }
 237   else if (strlen(axis) == 1)
 238     proj->SetXTitle(GetAxisTitle(hist, axis[0]));
 239
 240   TH1* proj2 = hist2->Project3D(axis);
 241   if (strlen(axis) == 2)
 242   {
 243     proj2->SetYTitle(GetAxisTitle(hist2, axis[0]));
 244     proj2->SetXTitle(GetAxisTitle(hist2, axis[1]));
 245   }
 246   else if (strlen(axis) == 1)
 247     proj2->SetXTitle(GetAxisTitle(hist2, axis[0]));
 248
 249   TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName())));
 250   //printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX()));
 251   //printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY()));
 252   division->Divide(proj, proj2, 1, 1, "B");
 253   division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle()));
 254
 255   if (putErrors)
 256   {
 257     division->Sumw2();
 258     if (division->GetDimension() == 1)
 259     {
 260       Int_t nBins = division->GetNbinsX();
 261       for (Int_t i = 1; i <= nBins; ++i)
 262         if (proj2->GetBinContent(i) != 0)
 263           division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i));
 264     }
 265     else if (division->GetDimension() == 2)
 266     {
 267       Int_t nBinsX = division->GetNbinsX();
 268       Int_t nBinsY = division->GetNbinsY();
 269       for (Int_t i = 1; i <= nBinsX; ++i)
 270         for (Int_t j = 1; j <= nBinsY; ++j)
 271           if (proj2->GetBinContent(i, j) != 0)
 272             division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j));
 273     }
 274   }
 275
 276   if (save)
 277   {
 278     proj->Write();
 279     proj2->Write();
 280     division->Write();
 281   }
 282 }
 283
 284 //____________________________________________________________________
 285 const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
 286 {
 287   // returns the title of the axis given in axis (x, y, z)
 288
 289   if (axis == 'x')
 290     return hist->GetXaxis()->GetTitle();
 291   else if (axis == 'y')
 292     return hist->GetYaxis()->GetTitle();
 293   else if (axis == 'z')
 294     return hist->GetZaxis()->GetTitle();
 295
 296   return 0;
 297 }
 298
 299
 300 AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 301
 302   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
 303
 304   if (!pythiaGenHeader) {
 305     printf("AliPWG0Helper::GetProcessType : Unknown gen Header type). \n");
 306     return kInvalidProcess;
 307   }
 308
 309
 310   Int_t pythiaType = pythiaGenHeader->ProcessType();
 311   fgLastProcessType = pythiaType;
 312   MCProcessType globalType = kInvalidProcess;
 313
 314
 315   if (adebug) {
 316     printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaType);
 317   }
 318
 319
 320   if(pythiaType==92||pythiaType==93){
 321     globalType = kSD;
 322   }
 323   else if(pythiaType==94){
 324     globalType = kDD;
 325   }
 326   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
 327   else {
 328     globalType = kND;
 329   }
 330   return globalType;
 331 }
 332
 333
 334 AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 335   //
 336   // get the process type of the event.
 337   //
 338
 339   // can only read pythia headers, either directly or from cocktalil header
 340   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
 341
 342   if (!dpmJetGenHeader) {
 343     printf("AliPWG0Helper::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
 344     return kInvalidProcess;
 345   }
 346
 347   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
 348   fgLastProcessType = dpmJetType;
 349   MCProcessType globalType = kInvalidProcess;
 350
 351
 352   if (adebug) {
 353     printf("AliPWG0Helper::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
 354   }
 355
 356
 357   if(dpmJetType == 1){ // this is explicitly inelastic
 358     globalType = kND;
 359   }
 360   else if(dpmJetType==5||dpmJetType==6){
 361     globalType = kSD;
 362   }
 363   else if(dpmJetType==7||dpmJetType==4){// DD and double pomeron
 364     globalType = kDD;
 365   }
 366   return globalType;
 367 }
 368
 369
 370 AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
 371   //
 372   // get the process type of the event.
 373   //
 374
 375
 376   // Check for simple headers first
 377
 378   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 379   if (pythiaGenHeader) {
 380     return GetPythiaEventProcessType(pythiaGenHeader,adebug);
 381   }
 382
 383   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader());
 384   if (dpmJetGenHeader) {
 385     return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
 386   }
 387
 388
 389   // check for cocktail
 390
 391   AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 392   if (!genCocktailHeader) {
 393     printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 394     return kInvalidProcess;
 395   }
 396
 397   TList* headerList = genCocktailHeader->GetHeaders();
 398   if (!headerList) {
 399     return kInvalidProcess;
 400   }
 401
 402   for (Int_t i=0; i<headerList->GetEntries(); i++) {
 403
 404     pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 405     if (pythiaGenHeader) {
 406       return GetPythiaEventProcessType(pythiaGenHeader,adebug);
 407     }
 408
 409     dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i));
 410     if (dpmJetGenHeader) {
 411       return GetDPMjetEventProcessType(dpmJetGenHeader,adebug);
 412     }
 413   }
 414   return kInvalidProcess;
 415 }
 416
 417
 418
 419 //____________________________________________________________________
 420 TParticle* AliPWG0Helper::FindPrimaryMother(AliStack* stack, Int_t label)
 421 {
 422   //
 423   // Finds the first mother among the primary particles of the particle identified by <label>,
 424   // i.e. the primary that "caused" this particle
 425   //
 426
 427   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
 428   if (motherLabel < 0)
 429     return 0;
 430
 431   return stack->Particle(motherLabel);
 432 }
 433
 434 //____________________________________________________________________
 435 Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
 436 {
 437   //
 438   // Finds the first mother among the primary particles of the particle identified by <label>,
 439   // i.e. the primary that "caused" this particle
 440   //
 441   // returns its label
 442   //
 443
 444   Int_t nPrim  = stack->GetNprimary();
 445
 446   while (label >= nPrim)
 447   {
 448     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
 449
 450     TParticle* particle = stack->Particle(label);
 451     if (!particle)
 452     {
 453       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
 454       return -1;
 455     }
 456
 457     // find mother
 458     if (particle->GetMother(0) < 0)
 459     {
 460       AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
 461       return -1;
 462     }
 463
 464     label = particle->GetMother(0);
 465   }
 466
 467   return label;
 468 }
 469
 470 //____________________________________________________________________
 471 void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
 472 {
 473   //
 474   // normalizes a 1-d histogram to its bin width
 475   //
 476
 477   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 478   {
 479     hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i));
 480     hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i));
 481   }
 482 }
 483
 484 //____________________________________________________________________
 485 void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
 486 {
 487   //
 488   // normalizes a 2-d histogram to its bin width (x width * y width)
 489   //
 490
 491   for (Int_t i=1; i<=hist->GetNbinsX(); ++i)
 492     for (Int_t j=1; j<=hist->GetNbinsY(); ++j)
 493     {
 494       Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j);
 495       hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor);
 496       hist->SetBinError(i, j, hist->GetBinError(i, j) / factor);
 497     }
 498 }
 499
 500 //____________________________________________________________________
 501 void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
 502 {
 503   //
 504   // Prints the given configuration
 505   //
 506
 507   TString str(">>>> Running with ");
 508
 509   switch (analysisMode)
 510   {
 511     case kInvalid: str += "invalid setting"; break;
 512     case kSPD : str += "SPD-only"; break;
 513     case kTPC : str += "TPC-only"; break;
 514     case kTPCITS : str += "Global tracking"; break;
 515   }
 516
 517   str += " and trigger ";
 518
 519   switch (trigger)
 520   {
 521     case kMB1 : str += "MB1"; break;
 522     case kMB2 : str += "MB2"; break;
 523     case kSPDFASTOR : str += "SPD FASTOR"; break;
 524   }
 525
 526   str += " <<<<";
 527
 528   Printf("%s", str.Data());
 529 }
 530