]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/ThermalFits/AliParticleYield.cxx
Added monitoring
[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / AliParticleYield.cxx
index fa50beb63c580d4f5c8614d12dfef0f72333b233..70311bfd82d675e97f6aa84198c0e25ece18a0f2 100644 (file)
@@ -129,15 +129,15 @@ fCentr(part.fCentr),
 fIsSum(part.fIsSum),
 fTag(part.fTag){
   // Copy constructor
+  
 }
 
 AliParticleYield::~AliParticleYield() {
 
 }
 
-TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const char * separators){
-  // Read the table from an ASCII File and returns a tree of particles. See ReadFromASCIIFile for detailed info on the format
-  TClonesArray * arr = ReadFromASCIIFile(fileName, separators);
+TTree * AliParticleYield::GetTreeFromArray(TClonesArray * arr) {
+  // Returns a tree from an array of Tparticles
   AliParticleYield * part = 0;
   TTree * tree = new TTree ("treePart", "Particle Yields and Ratios");
   tree->Branch("particles", &part);
@@ -145,13 +145,23 @@ TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const c
   while ((part = (AliParticleYield*) iterPart.Next())){
     tree->Fill();
   }
-  
+  if(part) delete part;
+  return tree;
+
+
+
+}
+
+
+TTree * AliParticleYield::ReadFromASCIIFileAsTree(const char * fileName, const char * separators){
+  // Read the table from an ASCII File and returns a tree of particles. See ReadFromASCIIFile for detailed info on the format
+  TClonesArray * arr = ReadFromASCIIFile(fileName, separators);
+  TTree * tree = GetTreeFromArray(arr);
   delete arr;
-  delete part;
   return tree;
 }
 
-TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TString selection) {
+TClonesArray * AliParticleYield::GetEntriesMatchingSelection(TTree * tree, TCut selection) {
   // Returns an array of particles from a tree created with ReadFromASCIIFileAsTree matching the selection. You can use the normal tree sintax for the selection, e.g. "fCentr == \"V0M0010\" && fStatus == 0".
 
   TClonesArray * arr = new TClonesArray("AliParticleYield");
@@ -197,7 +207,7 @@ TClonesArray * AliParticleYield::ReadFromASCIIFile(const char * fileName, const
   Int_t ipart = 0;
   std::cout << "Reading " << fileName << std::endl;
   
-  while (line.ReadLine(filestream) ) {
+  while (line.ReadLine(filestream) ) {    
     // Strip trailing and leading whitespaces
     line = line.Strip(TString::kLeading,  ' ');
     line = line.Strip(TString::kTrailing, ' ');
@@ -324,58 +334,115 @@ const char * AliParticleYield::GetLatexName(Int_t pdg) const {
 
   switch (pdg) {
   case 211:
-    if (fIsSum) return "(#pi^{+} + #pi^{-})";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#pi^{+} + #pi^{-})/2";
+      return "(#pi^{+} + #pi^{-})";
+    }
     return "#pi^{+}";
+    break;
   case -211:
     return "#pi^{-}";
+    break;
   case 321:
-    if (fIsSum) return "(K^{+} + K^{-})";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K^{+} + K^{-})/2";
+      return "(K^{+} + K^{-})";
+    }
     return "K^{+}";
+    break;
   case -321:
     return "K^{-}";
+    break;
   case 2212:
-    if (fIsSum) return "(p + #bar{p})";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(p + #bar{p})/2";
+      return "(p + #bar{p})";
+    }
     return "p";
+    break;
   case -2212:
     return "#bar^{p}";
+    break;
   case 3122:
-    if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
-    return "#Omega^{-}";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Lambda + #bar{#Lambda})/2";
+      return "(#Lambda + #bar{#Lambda})";
+    }
+    return "#Lambda";
+    break;
   case -3122:
-    return "#bar{#Omega}^{+}";
+    return "#bar{#Lamnba}";
+    break;
   case 3312:
-    if (fIsSum) return "(#Xi^{-} + #bar{#Xi}^{+})";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Xi^{-} + #bar{#Xi}^{+})/2";
+      return "(#Xi^{-} + #bar{#Xi}^{+})";
+    }
     return "#Xi^{-}";
+    break;
   case -3312:
     return "#bar{#Xi}^{+}";
+    break;
   case 3334:
-    if (fIsSum) return "(#Omega^{-} + #bar{#Omega}^{+})";
+    if (fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(#Omega^{-} + #bar{#Omega}^{+})/2";
+      return "(#Omega^{-} + #bar{#Omega}^{+})";
+    }
     return "#Omega^{-}";
+    break;
   case -3334:
     return "#bar{#Omega}^{+}";
+    break;
   case 310:
     return "K^{0}_{S}";
+    break;
   case 333:
     return "#phi";
+    break;
   case 313:
-    return "K^{*}";
+    if(fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(K* + #bar{K*})/2";
+      return "(K* + #bar{K*})";
+    }
+    return "K*";
+    break;
+  case -313:
+    return "#bar{K*}";
+    break;
   case 1000010020:
-    if(fIsSum) return "(d + #bar{d})";
+    if(fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(d + #bar{d})/2";
+      return "(d + #bar{d})";
+    }
     return "d";// Deuteron
+    break;
   case -1000010020:
     return "#bar{d}";// Deuteron
+    break;
   case 1000020030:
-    if(fIsSum) return "(^{3}He + #bar{^{3}He})";
+    if(fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "(^{3}He + #bar{^{3}He})/2";
+      return "(^{3}He + #bar{^{3}He})";
+    }
     return "^{3}He";
+    break;
   case -1000020030:
     return "#bar{^{3}He}";
+    break;
   case 1010010030:
-    if(fIsSum) return "^{3}_{#Lambda}H + #bar{^{3}_{#Lambda}H}";
-    return "^{3}_{#Lambda}H";
+    if(fIsSum) {
+      if (fMeasurementType & kTypeAveragePartAntiPart ) return "({}^{3}_{#Lambda}H + {}^{3}_{#Lambda}#bar{H})/2";
+      return "{}^{3}_{#Lambda}H + {}^{3}_{#Lambda}#bar{H}";
+    }
+    return "{}^{3}_{#Lambda}H";
+    break;
   case -1010010030:
-    return "#bar{^{3}_{#Lambda}H}";    
+    return "{}^{3}_{#Lambda}#bar{H}";    
+    break;
+  default:
+    AliWarning("Latex Name not know for this particle");
   }
-  AliWarning("Latex Name not know for this particle");
+
   return fPartName.Data();
 
 }
@@ -578,29 +645,68 @@ Bool_t AliParticleYield::CheckTypeConsistency() const {
     isOK= kFALSE;
   } else if (!(fMeasurementType & kTypeOnlyTotError) && (!GetStatError() || !GetSystError())) {
     AliError("Stat or syst errors are null");
+    isOK = kFALSE;
   } 
+  if((fMeasurementType & kTypeLinearInterpolation) && (fMeasurementType & kTypeAverageAndRefit) && (fMeasurementType & kTypeExtrPionRatio)) {
+    AliError("Only one out of the \"Liner interpolation\",  \"Average and refit\", \"Extrapolated with constant ratio to pions\" bits can be set"); 
+    isOK = kFALSE;
+  }
+  if((fMeasurementType & kTypeAveragePartAntiPart) && !fIsSum) {
+    AliError("Average part antipart set, but fIsSum is 0! This type bit should only be set for sums.");
+    isOK = kFALSE;
+  }
 
   return isOK;
 }
 
-void AliParticleYield::Print (Option_t *) const {
+void AliParticleYield::Print (Option_t *opt) const {
   // Print
-  Printf("-------------------------------");
-  CheckTypeConsistency();
-  Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fIsSum ? "particle + antiparticle" : "", fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
-  Printf("Status: %s, %s", kStatusString[fStatus], fMeasurementType & kTypeLinearInterpolation ? "Interpolated" : "Measured");
-  Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data()); 
-  if(fMeasurementType & kTypeOnlyTotError) {
-    Printf("%f +- %f (total error)", fYield, fSystError);
-  }
-  else {
-    Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
-  }
-  if(fNormErrorNeg) {
-    Printf("Normalization uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);    
+  // Available options: 
+  //  - short
+  //  - justvalue (does not print normalization error)
+  TString sopt(opt);
+  if(sopt.Contains("short")) {
+    printf("[%s]: %f +- %f +- %f ", fPartName.Data(), fYield, fStatError, fSystError);
+    if(fNormErrorNeg) {
+      printf("(+%f-%f)", fNormErrorPos, fNormErrorNeg);    
+    }else if(fNormErrorPos) {
+      printf("(+-%f)", fNormErrorPos);    
+    }
+    printf("[0x%8.8x,%d]\n", fMeasurementType, fStatus);
   }
-  else {
-    Printf("Normalization uncertainty: %f", fNormErrorPos);
+  else if (sopt.Contains("justvalue")) {
+    Printf("%f +- %f +- %f ", fYield, fStatError, fSystError);
+    
+  } else {
+    Printf("-------------------------------");
+    CheckTypeConsistency();
+    TString sumType = "";
+    if      (fIsSum && (fMeasurementType & kTypeAveragePartAntiPart)) sumType = "(particle + antiparticle)/2";
+    else if (fIsSum && !(fMeasurementType & kTypeAveragePartAntiPart)) sumType = "particle + antiparticle";
+    if(fMeasurementType & kTypeParticleRatio) {
+      Printf("%s [%s] (%d/%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, fPdgCode2, sumType.Data(), fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
+    }
+    else{ 
+      Printf("%s [%s] (%d) %s %s", fPartName.Data(), GetLatexName(), fPdgCode, sumType.Data(), fTag.Length() ? Form("[%s]", fTag.Data()) : "" );
+    }
+    TString measurementType = IsTypeMeasured() ? "Measured" : ""; 
+    if(fMeasurementType & kTypeLinearInterpolation) measurementType += "Interpolated";
+    if(fMeasurementType & kTypeAverageAndRefit) measurementType     += "Averaged+Refitted";
+    if(fMeasurementType & kTypeExtrPionRatio)   measurementType     += "Extrapolated assuming constant ratio to pions";
+    Printf("Status: %s, %s", kStatusString[fStatus],  measurementType.Data());
+    Printf("%s , sqrt(s) = %2.2f GeV, %2.2f < y < %2.2f %s", kSystemString[fCollisionSystem], fSqrtS, fYMin, fYMax, fCentr.Data()); 
+    if(fMeasurementType & kTypeOnlyTotError) {
+      Printf("%f +- %f (total error)", fYield, fSystError);
+    }
+    else {
+      Printf("%f +- %f (stat) +- %f (syst)", fYield, fStatError, fSystError);
+    }
+    if(fNormErrorNeg) {
+      Printf("Normalization uncertainty: +%f-%f", fNormErrorPos, fNormErrorNeg);    
+    }
+    else {
+      Printf("Normalization uncertainty: %f", fNormErrorPos);
+    }
   }
 }
 
@@ -732,12 +838,18 @@ AliParticleYield * AliParticleYield::FindParticle(TClonesArray * arr, Int_t pdg,
     }  
 
   }
+  if(!foundPart) {
+    Printf("ERROR<AliParticleYield::FindParticle>: Cannot find %d (System %d, sqrts = %2.2f GeV, %s, %s, Status:%d, pdg2:%d)", 
+           pdg, system, sqrts, centrality.Data(), isSum ? "part+antipart" : "", status, pdg2);
+  }
   return foundPart;
 
 }
 
-void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield*part2) {
+void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield*part2, const char * pdgSep) {
   // Combines metadata from part1 and part2
+  // pdgSep is a separator to be added in the name and pdg (e.g. + for a sum, / for a ratio)
+
   Int_t pdg1 = part1->GetPdgCode();
   Int_t pdg2 = pdg1 == part2->GetPdgCode() ? part1->GetPdgCode2() : part2->GetPdgCode();
   Int_t   system = part1->GetCollisionSystem() == part2->GetCollisionSystem() ? part2->GetCollisionSystem() : -1; 
@@ -745,11 +857,14 @@ void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield
   Int_t ymin = part1->GetYMin() == part2->GetYMin() ? part2->GetYMin() : -1000; 
   Int_t ymax = part1->GetYMax() == part2->GetYMax() ? part2->GetYMax() : -1000; 
   Int_t status = part1->GetStatus() == part2->GetStatus() ? part2->GetStatus() : -1; 
-  Int_t type = part1->GetMeasurementType() == part2->GetMeasurementType() ? part2->GetMeasurementType() : -1; 
+  Int_t type = part1->GetMeasurementType() | part2->GetMeasurementType();
   
-  TString centr = part1->GetCentr() == part2->GetCentr() ? part2->GetCentr() : part1->GetCentr()+"/"+part2->GetCentr(); 
-  TString tag = part1->GetTag() == part2->GetTag() ? part2->GetTag() : part1->GetTag()+"/"+part2->GetTag(); 
-  Int_t issum = part1->GetIsSum() == part2->GetIsSum() ? part2->GetIsSum() : -1000; 
+  TString centr = part1->GetCentr() == part2->GetCentr() ? part2->GetCentr() : part1->GetCentr()+pdgSep+part2->GetCentr(); 
+  TString tag = part1->GetTag() == part2->GetTag() ? part2->GetTag() : part1->GetTag()+pdgSep+part2->GetTag(); 
+  TString name = part1->GetPartName()+pdgSep+part2->GetPartName();
+
+  Int_t issum = part1->GetIsSum() || part2->GetIsSum() ? 1 : 0; 
+
 
   SetPdgCode(pdg1);
   SetPdgCode2(pdg2);
@@ -762,12 +877,13 @@ void AliParticleYield::CombineMetadata(AliParticleYield *part1, AliParticleYield
   SetCentr(centr);
   SetTag(tag);
   SetIsSum(issum);
+  SetPartName(name);
 
 }
 
 AliParticleYield * AliParticleYield::Add   (AliParticleYield * part1, AliParticleYield * part2, Double_t correlatedError , Option_t * opt){
 
-  // Computes the ratio of 2 particles.
+  // Computes the sum of 2 particles.
   // Valid options:
   //  - NQ: Propagates normalization errors quadratically (by default they are propagated linearly)
   //  - SL: propagates STATISTICAL errors linearly
@@ -798,8 +914,8 @@ AliParticleYield * AliParticleYield::Add   (AliParticleYield * part1, AliParticl
   part->SetStatError(stat);
   part->SetSystError(syst);
   part->SetNormError(norm);
-  part->CombineMetadata(part1, part2);
-
+  part->CombineMetadata(part1, part2, "+");
+  part->SetIsSum(1); // CombineMetadata inherits this form part1 and part2
   return part;
 
 
@@ -829,6 +945,7 @@ AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliPartic
     return 0;    
   }
 
+
   TString sopt(opt);
   sopt.ToUpper();
   if(part1->IsTypeRatio() || part2->IsTypeRatio()){
@@ -845,7 +962,10 @@ AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliPartic
   }
   
   if(correlatedError) {
+    std::cout << "Subtracting correlated error " << correlatedError << std::endl;
+    std::cout << "Before : " << syst << "[" << syst/value*100 <<"%]"<< std::endl;    
     syst = TMath::Sqrt(syst/value*syst/value -correlatedError*correlatedError)*value; // FIXME: this line was never tested
+    std::cout << "After  : " << syst << "[" << syst/value*100 <<"%]"<< std::endl;    
   }
 
   AliParticleYield * part = new AliParticleYield();
@@ -853,7 +973,9 @@ AliParticleYield * AliParticleYield::Divide (AliParticleYield * part1, AliPartic
   part->SetStatError(stat);
   part->SetSystError(syst);
   part->SetNormError(norm);
-  part->CombineMetadata(part1, part2);
+  part->CombineMetadata(part1, part2, "/");
+
+  part->SetMeasurementType(part->GetMeasurementType() | kTypeParticleRatio);// Set ratio bit
 
   return part;