]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliESDv0Analysis.C
Updates for embedding on AOD level
[u/mrichter/AliRoot.git] / STEER / AliESDv0Analysis.C
index bf2dd7d2c7e648619da495c95bcb3c6ebb0e0792..2351f87a1c0363f9b3c7a84027ce717f4d90e0fd 100644 (file)
@@ -8,41 +8,46 @@
 
 #if !defined( __CINT__) || defined(__MAKECINT__)
   #include <Riostream.h>
+  #include <TROOT.h>
   #include <TTree.h>
-  #include "TFile.h"
-  #include "TH1F.h"
-  #include "TH2F.h"
-  #include "TCanvas.h"
-  #include "TStopwatch.h"
-  #include "TParticle.h"
+  #include <TFile.h>
+  #include <TH1F.h>
+  #include <TCanvas.h>
 
-  #include "AliRun.h"
+  #include "AliESDEvent.h"
+  #include "AliESDv0.h"
+#endif
 
-  #include "AliESD.h"
+extern TROOT *gROOT;
 
-#endif
+Int_t AliESDv0Analysis(const Char_t *dir=".") { 
+   TH1F *hm=(TH1F*)gROOT->FindObject("hm");
+   if (!hm) {
+      hm=new TH1F("hm","Lambda+LambdaBar Effective Mass",60,1.065,1.165);
+      hm->SetXTitle("Mass (GeV/c**2)");
+   }
+   Char_t fname[100];
+   sprintf(fname,"%s/AliESDs.root",dir);
+   TFile *ef=TFile::Open(fname);
+   if (!ef||!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+   cerr<<"\n****** "<<fname<<" ******\n";
 
-extern AliRun *gAlice;
+   AliESDEvent* event = new AliESDEvent();
 
-Int_t AliESDv0Analysis(Int_t nev=1) { 
-   TH1F *hm=new TH1F("hm","Effective Mass",40,1.065,1.165);
-   hm->SetXTitle("Mass (GeV/c**2)");
 
-   TFile *ef=TFile::Open("AliESDs.root");
-   if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
-   AliESD* event = new AliESD;
    TTree* tree = (TTree*) ef->Get("esdTree");
    if (!tree) {cerr<<"no ESD tree found\n"; return 1;};
-   tree->SetBranchAddress("ESD", &event);
+   event->ReadFromTree(tree);
 
-   TStopwatch timer;
    Int_t rc=0,n=0;
 
    //****** Tentative particle type "concentrations"
-   Double_t c[5]={0.0, 0.0, 0.1, 0.1, 0.1};
+   Double_t c[5]={0.0, 0.0, 1, 0, 1};
+   AliPID pid;
+   pid.SetPriors(c);
 
    //******* The loop over events
-   while (tree->GetEvent(n)) {
+    while (tree->GetEvent(n)) {
 
      cerr<<"Processing event number : "<<n++<<endl;
 
@@ -51,40 +56,64 @@ Int_t AliESDv0Analysis(Int_t nev=1) {
 
      while (nv0--) {
        AliESDv0 *v0=event->GetV0(nv0);
-       Int_t pi=v0->GetPindex();
-       AliESDtrack *t=event->GetTrack(pi);
-       Int_t isProton=1;
-       if ((t->GetStatus()&AliESDtrack::kESDpid)!=0) {
-        Double_t r[10]; t->GetESDpid(r);
-         Double_t rcc=0.;
-         Int_t i;
-         for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
-         if (rcc==0.) continue;
-         //Here we apply Bayes' formula
-         Double_t w[10];
-         for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
-
-         if (w[4]<w[3]) isProton=0;
-         if (w[4]<w[2]) isProton=0;
-         if (w[4]<w[1]) isProton=0;
-        if (w[4]<w[0]) isProton=0;
-       }
-       if (!isProton) continue;
+
+       Int_t protonIdx=v0->GetPindex();
+       Int_t pionIdx  =v0->GetNindex();
+      
        v0->ChangeMassHypothesis(3122);
        Double_t mass=v0->GetEffMass();
+       if (mass>1.17) {  //check also the LambdaBar hypothesis
+          v0->ChangeMassHypothesis(-3122);
+          mass=v0->GetEffMass();
+          if (mass>1.17) continue;
+          Int_t tmp=protonIdx; protonIdx=pionIdx; pionIdx=tmp;
+       } 
+
+       AliESDtrack *protonTrk=event->GetTrack(protonIdx);
+       AliESDtrack *pionTrk  =event->GetTrack(pionIdx);
+
+       if (protonTrk->GetP()<0.5) continue;
+
+       // Check if the "proton track" is a proton
+       if ((protonTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
+        Double_t r[10]; protonTrk->GetESDpid(r);
+         pid.SetProbabilities(r);
+         Double_t pp=pid.GetProbability(AliPID::kProton);
+         if (pp < pid.GetProbability(AliPID::kElectron)) continue;
+         if (pp < pid.GetProbability(AliPID::kMuon)) continue;
+         if (pp < pid.GetProbability(AliPID::kPion)) continue;
+         if (pp < pid.GetProbability(AliPID::kKaon)) continue;
+       }
+       //Check if the "pion track" is a pion
+       if ((pionTrk->GetStatus()&AliESDtrack::kESDpid)!=0) {
+        Double_t r[10]; pionTrk->GetESDpid(r);
+         pid.SetProbabilities(r);
+         Double_t ppi=pid.GetProbability(AliPID::kPion);
+         if (ppi < pid.GetProbability(AliPID::kElectron)) continue;
+         if (ppi < pid.GetProbability(AliPID::kMuon)) continue;
+         if (ppi < pid.GetProbability(AliPID::kKaon)) continue;
+         if (ppi < pid.GetProbability(AliPID::kProton)) continue;
+        }
+
        hm->Fill(mass);
      } 
+
    }
 
    delete event;
+   delete tree;
    ef->Close();
 
-   timer.Stop(); timer.Print();
-
-   hm->Draw();
+   TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
+   if (!c1) {
+      c1=new TCanvas();
+   }
 
-   ef->Close();
+   if (hm->GetEntries()>100) hm->Fit("gaus","","",1.11,1.12);
+   else hm->Draw();
 
+   c1->Update();
 
    return rc;
 }