Refining the kRatio PID mode
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Feb 2010 21:21:40 +0000 (21:21 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 13 Feb 2010 21:21:40 +0000 (21:21 +0000)
PWG2/SPECTRA/AliAnalysisTaskProtons.cxx
PWG2/SPECTRA/AliAnalysisTaskProtons.h
PWG2/SPECTRA/AliProtonAnalysisBase.cxx
PWG2/SPECTRA/AliProtonAnalysisBase.h
PWG2/SPECTRA/macros/configProtonAnalysisBaseObject.C
PWG2/SPECTRA/macros/runProtonAnalysis.C

index 771ebc4..e662f12 100644 (file)
@@ -37,7 +37,7 @@ ClassImp(AliAnalysisTaskProtons)
 AliAnalysisTaskProtons::AliAnalysisTaskProtons()
   : AliAnalysisTask(), fESD(0), fAOD(0), fMC(0),
     fListAnalysis(0), fListQA(0), fHistEventStats(0), 
-    fProtonAnalysis(0) {
+  fProtonAnalysis(0), fCutCanvas(0) {
   //Dummy constructor
   
 }
@@ -46,7 +46,7 @@ AliAnalysisTaskProtons::AliAnalysisTaskProtons()
 AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name) 
   : AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0),
     fListAnalysis(0), fListQA(0), fHistEventStats(0), 
-    fProtonAnalysis(0) {
+    fProtonAnalysis(0), fCutCanvas(0) {
   // Constructor
   
   // Define input and output slots here
@@ -55,6 +55,7 @@ AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name)
   // Output slot #0 writes into a TList container
   DefineOutput(0, TList::Class());
   DefineOutput(1, TList::Class());
+  DefineOutput(2, TCanvas::Class());
 }
 
 //________________________________________________________________________
@@ -120,6 +121,8 @@ void AliAnalysisTaskProtons::CreateOutputObjects() {
   fListQA->SetName("fListQA");
   fListQA->Add(fProtonAnalysis->GetQAList());
   fListQA->Add(dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVertexQAList());
+
+  fCutCanvas = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetListOfCuts();
 }
 
 //________________________________________________________________________
@@ -235,6 +238,7 @@ void AliAnalysisTaskProtons::Exec(Option_t *) {
   // Post output data.
   PostData(0, fListAnalysis);
   PostData(1, fListQA);
+  PostData(2, fCutCanvas);
 }      
 
 //________________________________________________________________________
@@ -260,9 +264,9 @@ void AliAnalysisTaskProtons::Terminate(Option_t *) {
   c1->cd(2)->SetLeftMargin(0.15); c1->cd(2)->SetBottomMargin(0.15);  
   if (fHistYPtAntiProtons) fHistYPtAntiProtons->DrawCopy("colz");
 
-  TCanvas *c2 = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetListOfCuts();
+  /*TCanvas *c2 = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetListOfCuts();
   TFile *flocal = TFile::Open("ListOfCuts.root","recreate");
   c2->Write();
-  flocal->Close();
+  flocal->Close();*/
 }
 
index 089dca1..93df763 100644 (file)
@@ -9,6 +9,8 @@
 //-------------------------------------------------------------------------
 
 class TList;
+class TCanvas;
+
 class AliESDEvent;
 class AliAODEvent;
 class AliMCEvent;
@@ -41,6 +43,7 @@ class AliAnalysisTaskProtons : public AliAnalysisTask {
   TH1F   *fHistEventStats; //event statistics
 
   AliProtonAnalysis *fProtonAnalysis; //analysis object 
+  TCanvas *fCutCanvas; //Tcanvas with the analysis parameters (book-keeping)
   
   AliAnalysisTaskProtons(const AliAnalysisTaskProtons&); // not implemented
   AliAnalysisTaskProtons& operator=(const AliAnalysisTaskProtons&); // not implemented
index 9c99829..649f767 100644 (file)
@@ -72,16 +72,16 @@ AliProtonAnalysisBase::AliProtonAnalysisBase() :
   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
   fMinTPCdEdxPointsFlag(kFALSE),
   fFunctionProbabilityFlag(kFALSE), 
-  fNSigma(0),
+  fNSigma(0), fNRatio(0),
   fElectronFunction(0), fMuonFunction(0),
   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
   fDebugMode(kFALSE), fListVertexQA(new TList()) {
   //Default constructor
   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
-  for(Int_t i = 0; i < 24; i++) {
+  /*for(Int_t i = 0; i < 24; i++) {
     fdEdxMean[i] = 0.0;
     fdEdxSigma[i] = 0.0;
-  }
+    }*/
   fListVertexQA->SetName("fListVertexQA");
   TH1F *gHistVx = new TH1F("gHistVx",
                           "Vx distribution;V_{x} [cm];Entries",
@@ -685,24 +685,29 @@ TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
   l.DrawLatex(0.1,0.66,listOfCuts.Data());
   listOfCuts = "PID mode: "; 
   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
-  if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})"; 
-  if(fProtonPIDMode == kSigma1) {
-    listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
-    listOfCuts += " #sigma";
+  if(fProtonPIDMode == kRatio) {
+    listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.}) > ";
+    listOfCuts += fNRatio;
   }
-  if(fProtonPIDMode == kSigma2) {
-    listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
+  if(fProtonPIDMode == kSigma) {
+    listOfCuts += "N_{#sigma} area: "; listOfCuts += fNSigma;
     listOfCuts += " #sigma";
   }
+  //if(fProtonPIDMode == kSigma2) {
+  //listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
+  //listOfCuts += " #sigma";
+  //}
   l.DrawLatex(0.1,0.58,listOfCuts.Data());
   listOfCuts = "Accepted vertex diamond: "; 
-  l.DrawLatex(0.1,0.5,listOfCuts.Data());
+  l.DrawLatex(0.1,0.52,listOfCuts.Data());
   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
-  l.DrawLatex(0.6,0.5,listOfCuts.Data());
+  l.DrawLatex(0.6,0.52,listOfCuts.Data());
   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
-  l.DrawLatex(0.6,0.4,listOfCuts.Data());
+  l.DrawLatex(0.6,0.45,listOfCuts.Data());
   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
-  l.DrawLatex(0.6,0.3,listOfCuts.Data());
+  l.DrawLatex(0.6,0.38,listOfCuts.Data());
+  listOfCuts = "N_{contributors} > "; listOfCuts += fMinNumOfContributors; 
+  l.DrawLatex(0.6,0.31,listOfCuts.Data());
   listOfCuts = "Phase space: "; 
   l.DrawLatex(0.1,0.2,listOfCuts.Data());
   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";  
@@ -878,11 +883,38 @@ Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
   }
   //Ratio of the measured over the theoretical dE/dx a la STAR
   else if(fProtonPIDMode == kRatio) {
-    Printf("The kRatio mode is not implemented yet!!!");
-    return kFALSE;
-  }
+    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+    if(tpcTrack) {
+      gPt = tpcTrack->Pt();
+      gP = tpcTrack->P();
+      gEta = tpcTrack->Eta();
+    }
+    Double_t fAlephParameters[5];
+    if(fAnalysisMC) {
+      fAlephParameters[0] = 2.15898e+00/50.;
+      fAlephParameters[1] = 1.75295e+01;
+      fAlephParameters[2] = 3.40030e-09;
+      fAlephParameters[3] = 1.96178e+00;
+      fAlephParameters[4] = 3.91720e+00;
+    }
+    else {
+      fAlephParameters[0] = 0.0283086;
+      fAlephParameters[1] = 2.63394e+01;
+      fAlephParameters[2] = 5.04114e-11;
+      fAlephParameters[3] = 2.12543e+00;
+      fAlephParameters[4] = 4.88663e+00;
+    }
+    
+    AliESDpid *fESDpid = new AliESDpid(); 
+    AliTPCPIDResponse tpcResponse = fESDpid->GetTPCResponse(); 
+    tpcResponse.SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
+    Double_t normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse.GetExpectedSignal(gP,AliPID::kProton));
+    
+    if(normalizeddEdx >= fNRatio)
+      return kTRUE;
+  }//kRatio PID mode
   //Definition of an N-sigma area around the dE/dx vs P band
-  else if(fProtonPIDMode == kSigma1) {
+  else if(fProtonPIDMode == kSigma) {
     Double_t fAlephParameters[5];
     if(fAnalysisMC) {
       fAlephParameters[0] = 2.15898e+00/50.;
@@ -909,9 +941,9 @@ Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
   
     if(nsigma <= fNSigma) 
       return kTRUE;
-  }//kSigma1 PID method
+  }//kSigma PID method
   //Another definition of an N-sigma area around the dE/dx vs P band
-  else if(fProtonPIDMode == kSigma2) {
+  /*else if(fProtonPIDMode == kSigma2) {
     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
     if(tpcTrack) {
       gPt = tpcTrack->Pt();
@@ -941,48 +973,10 @@ Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
 
     if(normalizeddEdx >= -0.15)
       return kTRUE;
-  }
+  }*/
 
   return kFALSE;
 }
 
-//________________________________________________________________________
-void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
-  // This function is used in case the kSigma1 or kSigma2 PID mode is selected
-  // It takes as an argument the name of the ascii file (for the time being) 
-  // that is generated as a prior process.
-  // This ascii file has three columns: The min. P value (bins of 50MeV/c) 
-  // the mean and the sigma of the dE/dx distributions for protons coming 
-  // from a gaussian fit.
-  ifstream in;
-  in.open(filename);
-
-  Double_t gPtMin = 0.0;
-  Int_t iCounter = 0;
-  while(in.good()) {
-    in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
-    if(fDebugMode)
-      Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
-    iCounter += 1;
-  }
-}
-
-//________________________________________________________________________
-Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
-  // This is the Bethe-Bloch function normalised to 1 at the minimum
-  // We renormalize it based on the MC information
-  // WARNING: To be revised soon!!!
-  // This is just a temporary fix
-  Double_t normalization = 49.2;
-  Double_t bg2=bg*bg;
-  Double_t beta2 = bg2/(1.+ bg2);
-  Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
-  //
-  const Float_t kmeanCorrection =0.1;
-  Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
-  bb *= meanCorrection;
-
-  return normalization*bb; 
-}
 
 
index a54a3a3..93d3ada 100644 (file)
@@ -30,7 +30,7 @@ class AliProtonAnalysisBase : public TObject {
  public:
   enum TriggerMode { kMB1 = 0, kMB2, kSPDFASTOR };
   enum AnalysisMode { kInvalid = -1, kTPC = 0, kHybrid, kGlobal };
-  enum PIDMode { kBayesian = 0, kRatio, kSigma1, kSigma2 };
+  enum PIDMode { kBayesian = 0, kRatio, kSigma};
 
   AliProtonAnalysisBase();
   virtual ~AliProtonAnalysisBase();
@@ -247,7 +247,8 @@ class AliProtonAnalysisBase : public TObject {
   Bool_t IsProton(AliESDtrack *track);
   void SetNSigma(Int_t nsigma) {fNSigma = nsigma;}  
   Int_t GetNSigma() const {return fNSigma;}
-  void SetdEdxBandInfo(const char* filename);
+  void SetRatio(Double_t ratio) {fNRatio = ratio;}
+  Double_t GetRatio() {return fNRatio;}
   void SetPriorProbabilities(Double_t * const partFrac) {
     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) fPartFrac[i] = partFrac[i];} 
   void SetPriorProbabilityFunctions(TF1 *const felectron, 
@@ -261,7 +262,7 @@ class AliProtonAnalysisBase : public TObject {
   }
   Bool_t IsPriorProbabilityFunctionUsed() const {return fFunctionProbabilityFlag;}
   Double_t GetParticleFraction(Int_t i, Double_t p);
-  Double_t Bethe(Double_t bg) const;
+  //Double_t Bethe(Double_t bg) const;
 
   void SetDebugMode() {fDebugMode = kTRUE;}
   Bool_t GetDebugMode() const {return fDebugMode;}
@@ -320,8 +321,7 @@ class AliProtonAnalysisBase : public TObject {
   //pid
   Bool_t fFunctionProbabilityFlag; //flag: kTRUE if functions used
   Int_t fNSigma; //N-sigma cut in the dE/dx band
-  Double_t fdEdxMean[24]; //mean values of the dE/dx distributions for the proton band - P slices
-  Double_t fdEdxSigma[24]; //sigma values of the dE/dx distributions for the proton band - P slices
+  Double_t fNRatio; //min value of the ratio of the measured dE/dx vs the expected
   Double_t fPartFrac[10]; //prior probabilities
   TF1  *fElectronFunction; //momentum dependence of the prior probs
   TF1  *fMuonFunction; //momentum dependence of the prior probs
index 87f6b50..ea13eec 100644 (file)
@@ -96,16 +96,11 @@ AliProtonAnalysisBase *GetProtonAnalysisBaseObject(const char* analysisLevel = "
       break;\r
     case "Ratio":\r
       baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kRatio);\r
+      baseAnalysis->SetRatio(-0.2);\r
       break;\r
-    case "Sigma1":\r
+    case "Sigma":\r
       baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kSigma1);\r
       baseAnalysis->SetNSigma(4);\r
-      baseAnalysis->SetdEdxBandInfo("$ALICE_ROOT/PWG2/data/protonsdEdxInfo.dat");\r
-      break;\r
-    case "Sigma2":\r
-      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kSigma2);\r
-      baseAnalysis->SetNSigma(3);\r
-      baseAnalysis->SetdEdxBandInfo("$ALICE_ROOT/PWG2/data/protonsdEdxInfo.dat");\r
       break;\r
     default:\r
       break;\r
index 0eae932..50566bc 100644 (file)
@@ -1,6 +1,6 @@
 void runProtonAnalysis(Bool_t kAnalyzeMC = kTRUE,
                       const char* esdAnalysisType = "Hybrid",
-                      const char* pidMode = "Sigma1",
+                      const char* pidMode = "Ratio",
                       Bool_t kUseOnlineTrigger = kTRUE,
                       Bool_t kUseOfflineTrigger = kTRUE) {
   //Macro to run the proton analysis tested for local, proof & GRID.
@@ -28,8 +28,7 @@ void runProtonAnalysis(Bool_t kAnalyzeMC = kTRUE,
   //ESD analysis type can be one of the three: "TPC", "Hybrid", "Global"
   //PID mode can be one of the four: "Bayesian" (standard Bayesian approach) 
   //   "Ratio" (ratio of measured over expected/theoretical dE/dx a la STAR) 
-  //   "Sigma1" (N-sigma area around the fitted dE/dx vs P band)
-  //   "Sigma2" (same as previous but taking into account the No of TPC points)
+  //   "Sigma" (N-sigma area around the fitted dE/dx vs P band)
   TStopwatch timer;
   timer.Start();
   
@@ -57,11 +56,18 @@ void runLocal(const char* mode = "ESD",
              Bool_t kUseOfflineTrigger = kTRUE,
              const char* path = "/home/pchrist/ALICE/Alien/Tutorial/November2007/Tags") {
   TString smode = mode;
+  TString cutFilename = "ListOfCuts."; cutFilename += mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
+    cutFilename += "."; cutFilename += analysisType;
     outputFilename += "."; outputFilename += analysisType;
   }
-  outputFilename += ".root";
+  if(pidMode) {
+    cutFilename += "."; cutFilename += pidMode;
+    outputFilename += "."; outputFilename += pidMode;
+  }
+ cutFilename += ".root";
+ outputFilename += ".root";
 
   //____________________________________________________//
   //_____________Setting up the par files_______________//
@@ -129,11 +135,16 @@ void runLocal(const char* mode = "ESD",
                                                             TList::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cutCanvas",
+                                                            TCanvas::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                            cutFilename.Data());
 
   //____________________________________________//
   mgr->ConnectInput(taskProtons,0,cinput1);
   mgr->ConnectOutput(taskProtons,0,coutput1);
   mgr->ConnectOutput(taskProtons,1,coutput2);
+  mgr->ConnectOutput(taskProtons,2,coutput3);
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("local",chain);
@@ -148,12 +159,19 @@ void runInteractive(const char* mode = "ESD",
                    Bool_t kUseOfflineTrigger = kTRUE,
                    const char* collectionName = "tag.xml") {
   gSystem->Load("libProofPlayer.so");
-
+  
   TString smode = mode;
+  TString cutFilename = "ListOfCuts."; cutFilename += mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
+    cutFilename += "."; cutFilename += analysisType;
     outputFilename += "."; outputFilename += analysisType;
   }
+  if(pidMode) {
+    cutFilename += "."; cutFilename += pidMode;
+    outputFilename += "."; outputFilename += pidMode;
+  }
+  cutFilename += ".root";
   outputFilename += ".root";
 
   printf("*** Connect to AliEn ***\n");
@@ -229,11 +247,16 @@ void runInteractive(const char* mode = "ESD",
                                                             TList::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cutCanvas",
+                                                            TCanvas::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                            cutFilename.Data());
 
   //____________________________________________//
   mgr->ConnectInput(taskProtons,0,cinput1);
   mgr->ConnectOutput(taskProtons,0,coutput1);
   mgr->ConnectOutput(taskProtons,1,coutput2);
+  mgr->ConnectOutput(taskProtons,2,coutput3);
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("local",chain);
@@ -248,10 +271,17 @@ void runBatch(const char* mode = "ESD",
              Bool_t kUseOfflineTrigger = kTRUE,
              const char *collectionfile = "wn.xml") {
   TString smode = mode;
+  TString cutFilename = "ListOfCuts."; cutFilename += mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
+    cutFilename += "."; cutFilename += analysisType;
     outputFilename += "."; outputFilename += analysisType;
   }
+  if(pidMode) {
+    cutFilename += "."; cutFilename += pidMode;
+    outputFilename += "."; outputFilename += pidMode;
+  }
+  cutFilename += ".root";
   outputFilename += ".root";
 
   printf("*** Connect to AliEn ***\n");
@@ -319,11 +349,16 @@ void runBatch(const char* mode = "ESD",
                                                             TList::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cutCanvas",
+                                                            TCanvas::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                            cutFilename.Data());
   
   //____________________________________________//
   mgr->ConnectInput(taskProtons,0,cinput1);
   mgr->ConnectOutput(taskProtons,0,coutput1);
   mgr->ConnectOutput(taskProtons,1,coutput2);
+  mgr->ConnectOutput(taskProtons,2,coutput3);
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("grid",chain);
@@ -339,10 +374,17 @@ void runProof(const char* mode = "ESD",
              Int_t stats = 0, Int_t startingPoint = 0,
              const char* dataset = 0x0) {  
   TString smode = mode;
+  TString cutFilename = "ListOfCuts."; cutFilename += mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
+    cutFilename += "."; cutFilename += analysisType;
     outputFilename += "."; outputFilename += analysisType;
   }
+  if(pidMode) {
+    cutFilename += "."; cutFilename += pidMode;
+    outputFilename += "."; outputFilename += pidMode;
+  }
+  cutFilename += ".root";
   outputFilename += ".root";
 
   gEnv->SetValue("XSec.GSI.DelegProxy","2");
@@ -400,11 +442,16 @@ void runProof(const char* mode = "ESD",
                                                             TList::Class(),
                                                            AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("cutCanvas",
+                                                            TCanvas::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                            cutFilename.Data());
 
   //____________________________________________//
   mgr->ConnectInput(taskProtons,0,cinput1);
   mgr->ConnectOutput(taskProtons,0,coutput1);
   mgr->ConnectOutput(taskProtons,1,coutput2);
+  mgr->ConnectOutput(taskProtons,3,coutput3);
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();