function to reset stat error on corrections
[u/mrichter/AliRoot.git] / PWG0 / AliPWG0Helper.cxx
index 90c4a92..66ff632 100644 (file)
 
 #include <AliHeader.h>
 #include <AliStack.h>
-#include <AliLog.h>
 
 #include <AliLog.h>
 #include <AliESD.h>
 #include <AliESDEvent.h>
 #include <AliESDVertex.h>
+#include <AliVertexerTracks.h>
+#include <AliMultiplicity.h>
 
 #include <AliGenEventHeader.h>
 #include <AliGenPythiaEventHeader.h>
 #include <AliGenCocktailEventHeader.h>
 #include <AliGenDPMjetEventHeader.h>
+#include <AliESDVZERO.h>
 
 //____________________________________________________________________
 ClassImp(AliPWG0Helper)
 
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsEventTriggered(const AliESD* aEsd, Trigger trigger)
-{
-  // see function with ULong64_t argument
-
-  ULong64_t triggerMask = aEsd->GetTriggerMask();
-  return IsEventTriggered(triggerMask, trigger);
-}
+Int_t AliPWG0Helper::fgLastProcessType = -1;
 
 //____________________________________________________________________
-Bool_t AliPWG0Helper::IsEventTriggered(ULong64_t triggerMask, Trigger trigger)
+Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug)
 {
-  // check if the event was triggered
-  //
-  // this function needs the branch fTriggerMask
-  
+    // Checks if a vertex meets the needed quality criteria
 
-  // definitions from p-p.cfg
-  ULong64_t spdFO = (1 << 14);
-  ULong64_t v0left = (1 << 11);
-  ULong64_t v0right = (1 << 12);
-
-  switch (trigger)
+  Float_t requiredZResolution = -1;
+  if (analysisMode & kSPD || analysisMode & kTPCITS)
   {
-    case kMB1:
-    {
-      if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
-        return kTRUE;
-      break;
-    }
-    case kMB2:
-    {
-      if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
-        return kTRUE;
-      break;
-    }
+    requiredZResolution = 0.1;
   }
+  else if (analysisMode & kTPC)
+    requiredZResolution = 10.;
 
-  return kFALSE;
-}
-
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESD* aEsd)
-{
-  // see function with AliESDVertex argument
-
-  const AliESDVertex* vtxESD = aEsd->GetVertex();
-  return IsVertexReconstructed(vtxESD);
-}
-
-//____________________________________________________________________
-Bool_t AliPWG0Helper::IsVertexReconstructed(const AliESDVertex* vtxESD)
-{
-  // checks if the vertex is reasonable
-  //
-  // this function needs the branches fSPDVertex*
-
-  if (!vtxESD)
-    return kFALSE;
-
-  // the vertex should be reconstructed
-  if (strcmp(vtxESD->GetName(), "default")==0)
-    return kFALSE;
-
-  // check Ncontributors
-  if (vtxESD->GetNContributors() <= 0)
-    return kFALSE;
-
-  Double_t vtx_res[3];
-  vtx_res[0] = vtxESD->GetXRes();
-  vtx_res[1] = vtxESD->GetYRes();
-  vtx_res[2] = vtxESD->GetZRes();
+  // check resolution
+  Double_t zRes = vertex->GetZRes();
 
-  if (vtx_res[2]==0 || vtx_res[2]>0.1)
+  if (zRes > requiredZResolution) {
+    if (debug)
+      Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
     return kFALSE;
+  }
 
   return kTRUE;
 }
 
 //____________________________________________________________________
-const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug)
+const AliESDVertex* AliPWG0Helper::GetVertex(AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug, Bool_t bRedoTPC)
 {
   // Get the vertex from the ESD and returns it if the vertex is valid
   //
@@ -117,18 +67,30 @@ const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMo
   // also the quality criteria that are applied)
 
   const AliESDVertex* vertex = 0;
-  Float_t requiredZResolution = -1;
-  if (analysisMode == kSPD || analysisMode == kTPCITS)
+  if (analysisMode & kSPD || analysisMode & kTPCITS)
   {
     vertex = aEsd->GetPrimaryVertexSPD();
-    requiredZResolution = 0.1;
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning SPD vertex");
   }
-  else if (analysisMode == kTPC) 
+  else if (analysisMode & kTPC)
   {
+    if(bRedoTPC){
+      if (debug)
+        Printf("AliPWG0Helper::GetVertex: Redoing vertex");
+      Double_t kBz = aEsd->GetMagneticField();
+      AliVertexerTracks vertexer(kBz);
+      vertexer.SetTPCMode();
+      AliESDVertex *vTPC = vertexer.FindPrimaryVertex(aEsd);
+      aEsd->SetPrimaryVertexTPC(vTPC);
+      for (Int_t i=0; i<aEsd->GetNumberOfTracks(); i++) {
+       AliESDtrack *t = aEsd->GetTrack(i);
+       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
+      }
+      delete vTPC;
+    }
+
     vertex = aEsd->GetPrimaryVertexTPC();
-    requiredZResolution = 0.6;
     if (debug)
       Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks");
   }
@@ -143,22 +105,26 @@ const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMo
 
   // check Ncontributors
   if (vertex->GetNContributors() <= 0) {
-    if (debug)
-      Printf("AliPWG0Helper::GetVertex: NContributors() <= 0");
+    if (debug){
+      Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors());
+      Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices());
+      vertex->Print();
+    }
     return 0;
   }
 
   // check resolution
   Double_t zRes = vertex->GetZRes();
-
-  if (zRes == 0 || zRes > requiredZResolution) {
-    if (debug)
-      Printf("AliPWG0Helper::GetVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution);
+  if (zRes == 0) {
+    Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0.");
     return 0;
   }
 
   if (debug)
-    Printf("AliPWG0Helper::GetVertex: Returning valid vertex");
+  {
+    Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle());
+    vertex->Print();
+  }
 
   return vertex;
 }
@@ -184,6 +150,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
   }
 
   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
+  
 
   // skip quarks and gluon
   if (pdgCode <= 10 || pdgCode == 21)
@@ -193,10 +160,17 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
     return kFALSE;
   }
 
+  Int_t status = aParticle->GetStatusCode();
+  // skip non final state particles..
+  if(status!=1){
+    if (adebug)
+      printf("Dropping particle because it is not a final state particle.\n");
+    return kFALSE;
+  }
+
   if (strcmp(aParticle->GetName(),"XXX") == 0)
   {
-    if (adebug)
-      printf("WARNING: There is a particle named XXX.\n");
+    Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
     return kFALSE;
   }
 
@@ -204,8 +178,7 @@ Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimari
 
   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
   {
-    if (adebug)
-      printf("WARNING: There is a particle with an unknown particle class (pdg code %d).\n", pdgCode);
+    Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
     return kFALSE;
   }
 
@@ -330,7 +303,7 @@ const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis)
 }
 
 
-Int_t AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
 
   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
 
@@ -341,6 +314,7 @@ Int_t AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_
 
 
   Int_t pythiaType = pythiaGenHeader->ProcessType();
+  fgLastProcessType = pythiaType;
   MCProcessType globalType = kInvalidProcess;  
 
 
@@ -363,7 +337,7 @@ Int_t AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_
 }
 
 
-Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
   //
   // get the process type of the event.
   //
@@ -377,6 +351,7 @@ Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_
   }
 
   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
+  fgLastProcessType = dpmJetType;
   MCProcessType globalType = kInvalidProcess;  
 
 
@@ -398,7 +373,7 @@ Int_t AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_
 }
 
 
-Int_t AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
+AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) {
   //
   // get the process type of the event.
   //
@@ -498,45 +473,6 @@ Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
   return label;
 }
 
-void AliPWG0Helper::SetBranchStatusRecursive(TTree* tree, char *bname, Bool_t status, Bool_t debug)
-{
-  // Function  to switch on/off all data members of a top level branch
-  // this is needed for branches without a trailing dot ".", for those
-  // the root functionality with regular expressions does not work.
-  // Usage e.g.
-  // chain->SetBranchStatus("*", 0); 
-  // SetBranchStatusRecursive(chain,"SPDVertex",1);
-  // You need to give the full name of the top level branch zou want to access
-  //==========================================================================
-  // Author Christian.Klein-Boesing@cern.ch
-
-  if (!tree)
-    return;
-
-  TBranch *br = tree->GetBranch(bname);
-  if(!br) {
-    Printf("AliPWG0Helper::SetBranchStatusRecursive: Branch %s not found", bname);
-  }
-
-  TObjArray *leaves = tree->GetListOfLeaves();
-  Int_t nleaves = leaves->GetEntries();
-  TLeaf *leaf = 0;
-  TBranch *branch = 0;
-  TBranch *mother = 0;
-  for (Int_t i=0;i<nleaves;i++)  {
-    // the matched entry e.g. SPDVertex is its own Mother
-    leaf = (TLeaf*)leaves->UncheckedAt(i);
-    branch = (TBranch*)leaf->GetBranch();
-    mother = branch->GetMother();
-    if (mother==br){
-      if (debug)
-        Printf(">>>> AliPWG0Helper::SetBranchStatusRecursive: Setting status %s %s to %d", mother->GetName(), leaf->GetName(), status);
-      if (status) branch->ResetBit(kDoNotProcess);
-      else        branch->SetBit(kDoNotProcess);
-    }
-  }
-}
-
 //____________________________________________________________________
 void AliPWG0Helper::NormalizeToBinWidth(TH1* hist)
 {
@@ -568,31 +504,33 @@ void AliPWG0Helper::NormalizeToBinWidth(TH2* hist)
 }
 
 //____________________________________________________________________
-void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, Trigger trigger)
+void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger)
 {
   //
   // Prints the given configuration
   //
 
-  TString str(">>>> Running with ");
-
-  switch (analysisMode)
-  {
-    case kInvalid: str += "invalid setting"; break;
-    case kSPD : str += "SPD-only"; break;
-    case kTPC : str += "TPC-only"; break;
-    case kTPCITS : str += "Global tracking"; break;
-  }
+  TString str(">>>> Running with >");
 
-  str += " and trigger ";
+  if (analysisMode & kSPD)
+    str += "SPD-only";
+    
+  if (analysisMode & kTPC)
+     str += "TPC-only";
+    
+  if (analysisMode & kTPCITS)
+     str += "Global tracking";
 
-  switch (trigger)
+  if (analysisMode & kFieldOn)
   {
-    case kMB1 : str += "MB1"; break;
-    case kMB2 : str += "MB2"; break;
+     str += " (with field)";
   }
-
-  str += " <<<<";
+  else
+     str += " (WITHOUT field)";
+  
+  str += "< and trigger >";
+  str += AliTriggerAnalysis::GetTriggerName(trigger);
+  str += "< <<<<";
 
   Printf("%s", str.Data());
 }