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