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