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