new analysis to part corr for electron finding and b tagging
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Jul 2010 16:08:14 +0000 (16:08 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Jul 2010 16:08:14 +0000 (16:08 +0000)
PWG4/PWG4PartCorrDepLinkDef.h
PWG4/PartCorrDep/AliAnaBtag.cxx [new file with mode: 0644]
PWG4/PartCorrDep/AliAnaBtag.h [new file with mode: 0644]
PWG4/libPWG4PartCorrDep.pkg
PWG4/macros/AddTaskPartCorr.C
PWG4/macros/ConfigAnalysisBTag.C [new file with mode: 0644]

index 66ad6195a00e4c759b9c0163ea5d9356a5221038..75eb68debe00475bc36f20b50a343d7e19e52a71 100755 (executable)
@@ -18,6 +18,7 @@
 #pragma link C++ class AliAnaParticleJetLeadingConeCorrelation+;
 #pragma link C++ class AliAnaCalorimeterQA+;
 #pragma link C++ class AliAnaElectron+;
+#pragma link C++ class AliAnaBtag+;
 #pragma link C++ class AliAnalysisTaskTaggedPhotons+;
 #pragma link C++ class AliAnaOmegaToPi0Gamma+;
 
diff --git a/PWG4/PartCorrDep/AliAnaBtag.cxx b/PWG4/PartCorrDep/AliAnaBtag.cxx
new file mode 100644 (file)
index 0000000..cab8564
--- /dev/null
@@ -0,0 +1,1093 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+* Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes hereby granted      *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+/* $Id: $ */\r
+\r
+//_________________________________________________________________________\r
+//\r
+// Class for the electron identification and b-tagging.\r
+// Clusters from EMCAL matched to tracks\r
+// and kept in the AOD. Few histograms produced.\r
+//\r
+// -- Author: T.R.P-R.Aronsson (Yale) J.L. Klay (Cal Poly), M. Heinz (Yale)\r
+//////////////////////////////////////////////////////////////////////////////\r
+  \r
+// --- ROOT system --- \r
+#include <TH2F.h>\r
+#include <TH3F.h>\r
+#include <TParticle.h>\r
+#include <TNtuple.h>\r
+#include <TClonesArray.h>\r
+//#include <TObjString.h>\r
+//#include <Riostream.h>\r
+\r
+// --- Analysis system --- \r
+#include "AliAnaBtag.h" \r
+#include "AliCaloTrackReader.h"\r
+#include "AliMCAnalysisUtils.h"\r
+#include "AliAODCaloCluster.h"\r
+#include "AliFiducialCut.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODPid.h"\r
+#include "AliCaloPID.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliStack.h"\r
+#include "AliExternalTrackParam.h"\r
+#include "AliESDv0.h"\r
+#include "AliESDtrack.h"\r
+#include "AliAODJet.h"\r
+#include "AliAODEvent.h"\r
+#include "AliGenPythiaEventHeader.h"\r
+//#include "iostream.h"\r
+ClassImp(AliAnaBtag)\r
+  \r
+//____________________________________________________________________________\r
+AliAnaBtag::AliAnaBtag() \r
+: AliAnaPartCorrBaseClass(),fCalorimeter(""),\r
+  fpOverEmin(0.),fpOverEmax(0.),fResidualCut(0.),fMinClusEne(0.),\r
+  fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),\r
+  fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),\r
+  fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),fWriteNtuple(0),\r
+\r
+  fEleNtuple(0),\r
+  fhEmcalElectrons(0),fhTRDElectrons(0),fhTPCElectrons(0),fhDVM1(0),fhDVM2(0),fhJets(0),fhJetsAllEtaPhi(0),fhJetsLeadingBElectronEtaPhi(0),fhDVM1EtaPhi(0),fhBJetElectronDetaDphi(0),fhClusterEnergy(0),fhTestalle(0),fhResidual(0),fhElectrons(0),fhTracks(0)\r
+{\r
+  //default ctor\r
+  \r
+  //Initialize parameters\r
+  InitParameters();\r
+\r
+}\r
+/*\r
+//____________________________________________________________________________\r
+AliAnaBtag::AliAnaBtag(const AliAnaBtag & g) \r
+  : AliAnaPartCorrBaseClass(g),fCalorimeter(g.fCalorimeter),\r
+    fpOverEmin(g.fpOverEmin),fpOverEmax(g.fpOverEmax),\r
+    fResidualCut(g.fResidualCut),fMinClusEne(g.fMinClusEne),\r
+    fDrCut(g.fDrCut),fPairDcaCut(g.fPairDcaCut),fDecayLenCut(g.fDecayLenCut),fImpactCut(g.fImpactCut),\r
+    fAssocPtCut(g.fAssocPtCut),fMassCut(g.fMassCut),fSdcaCut(g.fSdcaCut),fITSCut(g.fITSCut),\r
+    fNTagTrkCut(g.fNTagTrkCut),fIPSigCut(g.fIPSigCut),\r
+    fJetEtaCut(g.fJetEtaCut),fJetPhiMin(g.fJetPhiMin),fJetPhiMax(g.fJetPhiMax),\r
+    fWriteNtuple(g.fWriteNtuple),\r
+\r
+    fEleNtuple(g.fEleNtuple),\r
+    fhEmcalElectrons(g.fhEmcalElectrons),fhTRDElectrons(g.fhTRDElectrons),fhTPCElectrons(g.fhTPCElectrons),fhDVM1(g.fhDVM1),fhDVM2(g.fhDVM2),fhJets(g.fhJets),fhJetsAllEtaPhi(g.fhJetsAllEtaPhi),fhJetsLeadingBElectronEtaPhi(g.fhJetsLeadingBElectronEtaPhi),fhDVM1EtaPhi(g.fhDVM1EtaPhi),fhBJetElectronDetaDphi(g.fhBJetElectronDetaDphi),fhClusterEnergy(g.fhClusterEnergy),fhTestalle(g.fhTestalle),fhResidual(g.fhResidual),fhElectrons(g.fhElectrons),fhTracks(g.fhTracks)\r
+{\r
+  // cpy ctor\r
+  \r
+}\r
+\r
+//_________________________________________________________________________\r
+AliAnaBtag & AliAnaBtag::operator = (const AliAnaBtag & g)\r
+{\r
+  // assignment operator\r
+  \r
+  if(&g == this) return *this;\r
+  fCalorimeter = g.fCalorimeter;\r
+  fpOverEmin = g.fpOverEmin;\r
+  fpOverEmax = g.fpOverEmax;\r
+  fResidualCut = g.fResidualCut;\r
+  fMinClusEne = g.fMinClusEne;\r
+  fDrCut = g.fDrCut;\r
+  fPairDcaCut = g.fPairDcaCut;\r
+  fDecayLenCut = g.fDecayLenCut;\r
+  fImpactCut = g.fImpactCut;\r
+  fAssocPtCut = g.fAssocPtCut;\r
+  fMassCut = g.fMassCut;\r
+  fSdcaCut = g.fSdcaCut;\r
+  fITSCut = g.fITSCut;\r
+  fNTagTrkCut = g.fNTagTrkCut;\r
+  fIPSigCut = g.fIPSigCut;\r
+  fJetEtaCut = g.fJetEtaCut;\r
+  fJetPhiMin = g.fJetPhiMin;\r
+  fJetPhiMax = g.fJetPhiMax;\r
+  fWriteNtuple = g.fWriteNtuple;\r
+\r
+  fEleNtuple = g.fEleNtuple; \r
+\r
+  fhEmcalElectrons = g.fhEmcalElectrons;\r
+  fhTRDElectrons = g.fhTRDElectrons;\r
+  fhTPCElectrons = g.fhTPCElectrons;\r
+  fhDVM1 = g.fhDVM1;\r
+  fhDVM2 = g.fhDVM2;\r
+  fhJets = g.fhJets;\r
+  fhJetsAllEtaPhi = g.fhJetsAllEtaPhi;\r
+  fhJetsLeadingBElectronEtaPhi = g.fhJetsLeadingBElectronEtaPhi;\r
+  fhDVM1EtaPhi = g.fhDVM1EtaPhi;\r
+  fhBJetElectronDetaDphi = g.fhBJetElectronDetaDphi;\r
+  fhClusterEnergy = g.fhClusterEnergy;\r
+  fhTestalle = g.fhTestalle;\r
+  fhResidual = g.fhResidual;\r
+  fhElectrons = g.fhElectrons;\r
+  fhTracks = g.fhTracks;\r
+  return *this;\r
+  \r
+}\r
+*/\r
+//____________________________________________________________________________\r
+AliAnaBtag::~AliAnaBtag() \r
+{\r
+  //dtor\r
+\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+TList *  AliAnaBtag::GetCreateOutputObjects()\r
+{  \r
+  // Create histograms to be saved in output file and \r
+  // store them in outputContainer\r
+  TList * outputContainer = new TList() ; \r
+  outputContainer->SetName("ElectronHistos") ; \r
+\r
+\r
+  if(IsDataMC()){\r
+\r
+    //electron ntuple for further analysis\r
+    if(fWriteNtuple) {\r
+      fEleNtuple = new TNtuple("EleNtuple","Electron Ntuple","fluff");\r
+      outputContainer->Add(fEleNtuple) ;\r
+    }\r
+\r
+\r
+\r
+  }//Histos with MC\r
+\r
+\r
+    fhEmcalElectrons = new TH1F("fhEmcalElectrons","",400,0,400);\r
+    outputContainer->Add(fhEmcalElectrons);\r
+    \r
+    fhTRDElectrons = new TH1F("fhTRDElectrons","",400,0,400);\r
+    outputContainer->Add(fhTRDElectrons);\r
+\r
+    fhTPCElectrons = new TH1F("fhTPCElectrons","",400,0,400);\r
+    outputContainer->Add(fhTPCElectrons);\r
+\r
+    fhDVM1 = new TH1F("fhDVM1","",400,0,400);\r
+    outputContainer->Add(fhDVM1);\r
+\r
+    fhDVM2 = new TH1F("fhDVM2","",400,0,400);\r
+    outputContainer->Add(fhDVM2);\r
+\r
+\r
+    fhJets = new TH2F("fhJets","",400,0,400,20,0,20);\r
+    outputContainer->Add(fhJets);\r
+\r
+\r
+    //Big one:\r
+    //1 is All Jets\r
+    //2 empty\r
+    //3 Jets within pt 10 cut\r
+    //4 Jets With geometric cut\r
+    //5 Leading jet\r
+    //6 DVM jet\r
+    //9 All identified tracks as electrons\r
+    //10 is all track Pt()\r
+    //11 Tracks "in"\r
+    //12 is Cluster energy plot\r
+\r
+\r
+\r
+\r
+    fhJetsAllEtaPhi = new TH2F("fhJetsAllEtaPhi","",100,-2,2,100,-2,8);\r
+    outputContainer->Add(fhJetsAllEtaPhi);\r
+\r
+    fhJetsLeadingBElectronEtaPhi = new TH2F("fhJetsLeadingBElectronEtaPhi","",100,-5,5,200,-10,10);\r
+    outputContainer->Add(fhJetsLeadingBElectronEtaPhi);\r
+\r
+    fhDVM1EtaPhi = new TH2F("fhDVM1EtaPhi","",100,-2,2,100,-2,8);\r
+    outputContainer->Add(fhDVM1EtaPhi);\r
+\r
+    fhBJetElectronDetaDphi = new TH2F("fhBJetElectronDetaDphi","",100,-5,5,200,-10,10);\r
+    outputContainer->Add(fhBJetElectronDetaDphi);\r
+\r
+    fhClusterEnergy = new TH1F("fhClusterEnergy","",100,0,10);\r
+    outputContainer->Add(fhClusterEnergy);\r
+\r
+    fhTestalle = new TH1F("fhTestalle","",400,0,400);\r
+    outputContainer->Add(fhTestalle);\r
+\r
+    fhResidual = new TH1F("fhResidual","",500,0,5);\r
+    outputContainer->Add(fhResidual);\r
+\r
+    fhElectrons = new TH2F("fhElectrons","",200,0,100,20,0,20);\r
+    outputContainer->Add(fhElectrons);\r
+\r
+    fhTracks = new TH2F("fhTracks","",200,0,100,20,0,20);\r
+    outputContainer->Add(fhTracks);\r
+\r
+\r
+  \r
+  return outputContainer ;\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaBtag::Init()\r
+{\r
+\r
+  //do some initialization\r
+  if(fCalorimeter == "PHOS") {\r
+    printf("AliAnaBtag::Init() - !!STOP: You want to use PHOS in analysis but this is not (yet) supported!!\n!!Check the configuration file!!\n");\r
+    fCalorimeter = "EMCAL";\r
+  }\r
+  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){\r
+    printf("AliAnaBtag::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!!\n!!Check the configuration file!!\n");\r
+    abort();\r
+  }\r
+\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+void AliAnaBtag::InitParameters()\r
+{\r
+  \r
+  //Initialize the parameters of the analysis.\r
+  SetOutputAODClassName("AliAODPWG4Particle");\r
+  SetOutputAODName("PWG4Particle");\r
+\r
+  AddToHistogramsName("Btag_");\r
+\r
+  fCalorimeter = "EMCAL" ;\r
+  fpOverEmin = 0.5;\r
+  fpOverEmax = 1.2;\r
+  fResidualCut = 0.02;\r
+  fMinClusEne = 4.0;\r
+  //DVM B-tagging\r
+  fDrCut       = 1.0; \r
+  fPairDcaCut  = 0.02;\r
+  fDecayLenCut = 1.0;\r
+  fImpactCut   = 0.5;\r
+  fAssocPtCut  = 1.0;\r
+  fMassCut     = 1.5;\r
+  fSdcaCut     = 0.1;\r
+  fITSCut      = 4;\r
+  //IPSig B-tagging\r
+  fNTagTrkCut  = 2;\r
+  fIPSigCut    = 3.0;\r
+  //Jet fiducial cuts\r
+  fJetEtaCut = 0.3;\r
+  fJetPhiMin = 1.8;\r
+  fJetPhiMax = 2.9;\r
+}\r
+\r
+//__________________________________________________________________\r
+void  AliAnaBtag::MakeAnalysisFillAOD() \r
+{\r
+  //This reads in tracks, extrapolates to EMCAL, does p/E selectrons, identifies electron candidates\r
+  //After candidates are obtained, btagging and saving into AOD.\r
+\r
+  Double_t bfield = 0.;\r
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
+\r
+  //Select the calorimeter of the electron\r
+  if(fCalorimeter != "EMCAL") {\r
+    printf("This class not yet implemented for PHOS\n");\r
+    abort();\r
+  }\r
+  \r
+  TObjArray *cl = GetAODEMCAL();\r
+  \r
+  if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;\r
+  Int_t ntracks = GetAODCTS()->GetEntriesFast();\r
+  if(GetDebug() > 0)\r
+    printf("AliAnaBtag::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);\r
+\r
+  Int_t iCluster = -999;\r
+  Int_t ntot = cl->GetEntriesFast();\r
+\r
+  //CLUSTER STUFF\r
+  if(1){\r
+    for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
+      AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));\r
+      if(!clus) continue;\r
+      //Double_t x[3];\r
+      //clus->GetPosition(x);\r
+      //if(clus->E()>3.) fhJets->Fill(clus->E(),12);\r
+      fhClusterEnergy->Fill(clus->E());\r
+    }\r
+  }\r
+\r
+\r
+  for (Int_t itrk =  0; itrk <  ntracks; itrk++) {////////////// track loop\r
+    iCluster = -999; //start with no match\r
+    AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(itrk)) ;\r
+\r
+    Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
+    Bool_t dcaOkay = GetDCA(track,imp,cov);  //homegrown dca calculation until AOD is fixed\r
+    if(!dcaOkay) printf("AliAnaBtag::Problem computing DCA to primary vertex for track %d.  Skipping it...\n",itrk);\r
+\r
+    fhTracks->Fill(track->Pt(),1);\r
+\r
+    AliAODPid* pid = (AliAODPid*) track->GetDetPid();\r
+    if(pid == 0) {\r
+      if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);\r
+      continue;\r
+    } \r
+    else {\r
+      Double_t emcpos[3];\r
+      pid->GetEMCALPosition(emcpos);\r
+      Double_t emcmom[3];\r
+      pid->GetEMCALMomentum(emcmom);\r
+      \r
+      TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);\r
+      TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);\r
+      Double_t tphi = pos.Phi();\r
+      Double_t teta = pos.Eta();\r
+      Double_t tmom = mom.Mag();\r
+\r
+      Bool_t in = kFALSE;\r
+      if(mom.Phi()*180./TMath::Pi() > 80. && mom.Phi()*180./TMath::Pi() < 190. &&\r
+        mom.Eta() > -0.7 && mom.Eta() < 0.7) in = kTRUE; //Checks momentum?\r
+      //Also check the track\r
+      if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. &&\r
+        track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE;\r
+\r
+      if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD() - Track(Extrap) pt %2.2f(%2.2f), phi %2.2f(%2.2f), eta %2.2f(%2.2f) in fiducial cut %d\n",track->Pt(), mom.Pt(), track->Phi(), mom.Phi(), track->Eta(),mom.Eta(), in);\r
+\r
+      if(mom.Pt() > GetMinPt() && in) {\r
+       fhTracks->Fill(track->Pt(),3);\r
+       Double_t dEdx = pid->GetTPCsignal();\r
+               \r
+       //Int_t ntot = cl->GetEntriesFast();\r
+       Double_t res = 999.;\r
+       Double_t pOverE = -999.;\r
+       \r
+       Int_t pidProb = track->GetMostProbablePID();\r
+       Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;\r
+       Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;\r
+       Bool_t trkChgHad = kFALSE; if(pidProb == AliAODTrack::kPion || pidProb == AliAODTrack::kKaon || pidProb == AliAODTrack::kProton) trkChgHad = kTRUE;\r
+\r
+\r
+       //Track Matching!\r
+       Bool_t emcEle = kFALSE;      \r
+       double minRes=100.;\r
+       Double_t minR  = 99;\r
+        Double_t minPe =-1;\r
+        Double_t minEp =-1;\r
+        Double_t minMult = -1;\r
+        Double_t minPt   = -1;\r
+\r
+       for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
+         AliAODCaloCluster * clus = (AliAODCaloCluster*) (cl->At(iclus));\r
+         if(!clus) continue;\r
+\r
+         //As of 11-Oct-2009\r
+         //only select "good" clusters   \r
+//       if (clus->GetNCells()       < 2    ) continue;\r
+//           if (clus->GetNCells()       > 30   ) continue;\r
+//           if (clus->E()               < fMinClusEne ) continue;\r
+//           if (clus->GetDispersion()   > 1    ) continue;\r
+//           if (clus->GetM20()          > 0.4  ) continue;\r
+//           if (clus->GetM02()          > 0.4  ) continue;\r
+//           if (clus->GetM20()          < 0.03 ) continue;\r
+//           if (clus->GetM02()          < 0.03 ) continue;\r
+// new optimized from ben. 2010May\r
+         if (clus->GetNCells()       < 2    ) continue;\r
+          if (clus->GetNCells()       > 35   ) continue;\r
+          if (clus->E()               < 0 ) continue;\r
+          if (clus->GetDispersion()   > 1.08    ) continue;\r
+          if (clus->GetM20()          > 0.42  ) continue;\r
+          if (clus->GetM02()          > 0.4  ) continue;\r
+          if (clus->GetM20()          < 0 ) continue;\r
+          if (clus->GetM02()          < 0.06 ) continue;\r
+         \r
+         \r
+         Double_t x[3];\r
+         clus->GetPosition(x);\r
+         TVector3 cluspos(x[0],x[1],x[2]);\r
+         Double_t deta = teta - cluspos.Eta();\r
+         Double_t dphi = tphi - cluspos.Phi();\r
+         if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
+         if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
+\r
+         res = sqrt(dphi*dphi + deta*deta);\r
+         \r
+         if(res<minRes){\r
+           minRes=res;\r
+           \r
+         }\r
+         \r
+         if(res < 0.0275) { //   if(res < fResidualCut) { //Optimized from Ben\r
+           iCluster = iclus;\r
+           \r
+           Double_t energy = clus->E(); \r
+           if(energy > 0) pOverE = tmom/energy;\r
+           \r
+           if (res< minR) {\r
+              minR  = res;\r
+              minPe = pOverE;\r
+              minEp = energy/tmom;\r
+              minMult = clus->GetNCells() ;\r
+              minPt = track->Pt();\r
+            }\r
+\r
+           \r
+           \r
+           if(fWriteNtuple) {\r
+             fEleNtuple->Fill(1);\r
+           }\r
+\r
+         } else {\r
+             //unmatched\r
+         }//res cut\r
+       }//calo cluster loop\r
+\r
+       fhResidual->Fill(minRes);\r
+\r
+       if(minPe > 0.9 && minPe < 1.08) emcEle = kTRUE;//       if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;\r
+       \r
+\r
+       \r
+       if(emcEle)\r
+         fhEmcalElectrons->Fill(track->Pt());\r
+       if(trkEle)\r
+         fhTRDElectrons->Fill(track->Pt());\r
+       if(tpcEle)\r
+         fhTPCElectrons->Fill(track->Pt());\r
+       \r
+       \r
+       //Take all emcal electrons, but the others only if pT < 10 GeV\r
+       if(emcEle && (tpcEle || trkEle) ) {\r
+         fhTestalle->Fill(track->Pt());\r
+         //B-tagging\r
+         if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");\r
+         Int_t dvmbtag = GetDVMBtag(track);\r
+         \r
+         \r
+         if(dvmbtag>0){\r
+           fhDVM1->Fill(track->Pt());\r
+            fhDVM1EtaPhi->Fill(track->Eta(),track->Phi());       \r
+           \r
+          }\r
+         if(dvmbtag>1)\r
+           fhDVM2->Fill(track->Pt());\r
+         \r
+\r
+         //Create particle to save in AOD\r
+         Double_t eMass = 0.511/1000; //mass in GeV\r
+         Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);\r
+         AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);\r
+         tr.SetLabel(track->GetLabel());\r
+         tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters\r
+         tr.SetTrackLabel(track->GetID(),-1); //sets the indices of the original tracks tr.SetTrackLabel(track->GetID(),-1) instead of itrk;\r
+\r
+         if(emcEle) {//PID determined by EMCAL\r
+           tr.SetDetector(fCalorimeter);\r
+         } else {\r
+           tr.SetDetector("CTS"); //PID determined by CTS\r
+         }\r
+\r
+         if(GetReader()->GetAODCTSNormalInputEntries() <= itrk) tr.SetInputFileIndex(1);\r
+         //Make this preserve sign of particle\r
+         if(track->Charge() < 0) tr.SetPdg(11); //electron is 11\r
+         else  tr.SetPdg(-11); //positron is -11\r
+\r
+         tr.SetBtag(dvmbtag);\r
+         \r
+\r
+         //Check origin of the candiates\r
+         if(IsDataMC()){\r
+           tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex()));\r
+           if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD() - Origin of candidate %d\n",tr.GetTag());\r
+         }//Work with stack also   \r
+         \r
+         AddAODParticle(tr);\r
+         \r
+         if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD() - Electron selection cuts passed: pT %3.2f, pdg %d\n",tr.Pt(),tr.GetPdg());      \r
+       }//electron\r
+      }//pt, fiducial selection\r
+    }//pid check\r
+  }//track loop                         \r
+  \r
+\r
+\r
+\r
+  if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD()  End fill AODs \n");  \r
+  \r
+}\r
+\r
+//__________________________________________________________________\r
+void  AliAnaBtag::MakeAnalysisFillHistograms() \r
+{\r
+  //Do analysis and fill histograms\r
+\r
+  AliStack * stack = 0x0;\r
+//   TParticle * primary = 0x0;\r
+\r
+\r
+\r
+\r
+  if(IsDataMC()) {\r
+    if(GetReader()->ReadStack()){\r
+      stack =  GetMCStack() ;      \r
+      if(!stack)\r
+       printf("AliAnaBtag::MakeAnalysisFillHistograms() *** no stack ***: \n");\r
+      \r
+    }\r
+  }// is data and MC\r
+\r
+  ////////////////////////////////////\r
+  //Loop over jets and check for b-tag\r
+  ////////////////////////////////////\r
+  double maxjetEta=-4.;\r
+  double maxjetPhi=-4.;\r
+\r
+  Int_t njets = (GetReader()->GetOutputEvent())->GetNJets();\r
+  if(njets > 0) {\r
+    if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets.  Performing b-jet tag analysis\n",njets);\r
+\r
+    for(Int_t ijet = 0; ijet < njets ; ijet++) {\r
+      AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;\r
+\r
+      if(ijet==0){\r
+        maxjetEta=jet->Eta();\r
+        maxjetPhi=jet->Phi();\r
+      }\r
+\r
+      fhJets->Fill(jet->Pt(),1);\r
+      fhJetsAllEtaPhi->Fill(jet->Eta(),jet->Phi());\r
+\r
+      if(jet->Pt() < 0.) continue; //This has to be adjusted depending on pp or AA!\r
+      fhJets->Fill(jet->Pt(),3); //All jets after pt cut\r
+\r
+      //Geometric EMCAL cut\r
+      if(TMath::Abs(jet->Eta()) > fJetEtaCut) continue;\r
+      if(jet->Phi() < fJetPhiMin || jet->Phi() > fJetPhiMax) continue;\r
+      fhJets->Fill(jet->Pt(),4); //All jets after geometric cut\r
+\r
+      Bool_t leadJet  = kFALSE;\r
+      if (ijet==0){ \r
+       fhJets->Fill(jet->Pt(),5); //Leading jets\r
+       leadJet= kTRUE;\r
+      }\r
+\r
+\r
+      Bool_t dvmJet = kFALSE;  \r
+      TRefArray* rt = jet->GetRefTracks();\r
+      Int_t ntrk = rt->GetEntries();\r
+\r
+      for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
+       AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
+       Bool_t isDVM = CheckIfBjet(jetTrack);\r
+       if(isDVM) dvmJet = kTRUE;\r
+      }\r
+\r
+      if(dvmJet)\r
+       fhJets->Fill(jet->Pt(),6);\r
+\r
+\r
+\r
+      if(IsDataMC()) {\r
+       //determine tagging efficiency & mis-tagging rate\r
+       //using b-quarks from stack\r
+       Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());\r
+       Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());\r
+       if (isTrueBjet && GetDebug() > 0) printf("== True Bjet==\n");\r
+       if (isTrueDjet && GetDebug() > 0) printf("== True Charm-jet==\n");\r
+       if (dvmJet && GetDebug() > 0)     printf("== found DVM jet==\n");\r
+      }\r
+\r
+    } //jet loop\r
+  } //jets exist\r
+  \r
+\r
+\r
+\r
+  //Electron loop, read back electrons, fill histos\r
+  Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
+  if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);\r
+  \r
+  for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+    AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
+    Int_t pdg = ele->GetPdg();\r
+\r
+\r
+    if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue; //not necessary..\r
+\r
+    //MC tag of this electron\r
+    //    Int_t mctag = ele->GetTag();\r
+\r
+\r
+    fhElectrons->Fill(ele->Pt(),1); //All electrons\r
+    Bool_t photonic = kFALSE;\r
+    photonic = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
+    if(!photonic) fhElectrons->Fill(ele->Pt(),3); //nonphotonic electrons\r
+    if(photonic) fhElectrons->Fill(ele->Pt(),4);  //photonic electrons\r
+\r
+    //Fill electron histograms \r
+    Float_t phiele = ele->Phi();\r
+    Float_t etaele = ele->Eta();\r
+\r
+\r
+    if(ele->GetBtag()>0){ // removed bit tag shit\r
+      fhElectrons->Fill(ele->Pt(),5);\r
+      if(!photonic) fhElectrons->Fill(ele->Pt(),6);\r
+      if(photonic) fhElectrons->Fill(ele->Pt(),7);\r
+      fhJetsLeadingBElectronEtaPhi->Fill(maxjetEta,maxjetPhi); \r
+      double deta=etaele-maxjetEta;\r
+      double dphi=phiele-maxjetPhi;\r
+      \r
+      //double r = sqrt((dphi*dphi)+(deta*deta));\r
+      fhBJetElectronDetaDphi->Fill(deta,dphi);\r
+      \r
+    }\r
+    \r
+  }//electron aod loop\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Int_t AliAnaBtag::GetDVMBtag(AliAODTrack * tr )\r
+{\r
+  //This method uses the Displaced Vertex between electron-hadron\r
+  //pairs and the primary vertex to determine whether an electron is\r
+  //likely from a B hadron.\r
+\r
+  Int_t ncls1 = 0;\r
+  for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;\r
+\r
+\r
+  if (ncls1 < fITSCut) return 0;\r
+\r
+  Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
+  Bool_t dcaOkay = GetDCA(tr,imp,cov);  //homegrown dca calculation until AOD is fixed                  \r
+  if(!dcaOkay) {\r
+    printf("AliAnaBtag::Problem computing DCA to primary vertex for track %d",tr->GetID());\r
+    return 0;\r
+  }\r
+\r
+  if (TMath::Abs(imp[0])   > fImpactCut ) return 0;\r
+  if (TMath::Abs(imp[1])   > fImpactCut ) return 0;\r
+\r
+  Int_t nvtx1 = 0;\r
+  Int_t nvtx2 = 0;\r
+  Int_t nvtx3 = 0;\r
+\r
+  for (Int_t k2 =0; k2 < GetAODCTS()->GetEntriesFast() ; k2++) {\r
+    //loop over assoc\r
+    AliAODTrack* track2 = (AliAODTrack*)GetAODCTS()->At(k2);\r
+    Int_t id1 = tr->GetID();\r
+    Int_t id2 = track2->GetID();\r
+    if(id1 == id2) continue;\r
+\r
+    Int_t ncls2 = 0;\r
+    for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;\r
+    if (ncls2 < fITSCut) continue;\r
+\r
+    if(tr->Pt()<6.&&track2->Pt() < 0.4) continue;\r
+    if(tr->Pt()>6.&&tr->Pt()<10.&&track2->Pt() < 0.2) continue;\r
+    if(tr->Pt()>10.&&track2->Pt() < 0.6) continue;\r
+\r
+\r
+    Double_t dphi = tr->Phi() - track2->Phi();\r
+    if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
+    if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
+    Double_t deta = tr->Eta() - track2->Eta();\r
+    Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+\r
+    if(dr > fDrCut) continue;\r
+    \r
+    if(tr->Pt()<6.){\r
+      Double_t sDca = ComputeSignDca(tr, track2, 1.4,0.025);\r
+       if(sDca > 0.06) nvtx2++;\r
+    } \r
+    if(tr->Pt()>6.&&tr->Pt()<10.){\r
+      Double_t sDca = ComputeSignDca(tr, track2, 1.7,0.012);\r
+        if(sDca > 0.03) nvtx2++;\r
+    } \r
+    if(tr->Pt()>10.){\r
+      Double_t sDca = ComputeSignDca(tr, track2, 1.5,0.14);\r
+        if(sDca > 0.04) nvtx2++;\r
+    } \r
+\r
+\r
+  } //loop over hadrons\r
+\r
+  if(GetDebug() > 0) {\r
+    if (nvtx1>0) printf("result1 of btagging: %d \n",nvtx1);\r
+    if (nvtx2>0) printf("result2 of btagging: %d \n",nvtx2);\r
+    if (nvtx3>0) printf("result3 of btagging: %d \n",nvtx3);\r
+  }\r
+\r
+\r
+\r
+\r
+  return nvtx2;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Double_t AliAnaBtag::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , float masscut, double pdcacut)\r
+{\r
+  //Compute the signed dca between two tracks\r
+  //and return the result\r
+\r
+  Double_t signDca=-999.;\r
+  if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);\r
+\r
+  //=====Now calculate DCA between both tracks=======  \r
+  Double_t massE = 0.000511;\r
+  Double_t massK = 0.493677;\r
+\r
+  Double_t vertex[3] = {-999.,-999.,-999}; //vertex\r
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
+    GetReader()->GetVertex(vertex); //If only one file, get the vertex from there\r
+    //FIXME:  Add a check for whether file 2 is PYTHIA or HIJING\r
+    //If PYTHIA, then set the vertex from file 2, if not, use the\r
+    //vertex from file 1\r
+    if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex);\r
+  }\r
+  \r
+  TVector3 primV(vertex[0],vertex[1],vertex[2]) ;\r
+\r
+  if(GetDebug() > 5) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;\r
+\r
+  AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);\r
+  AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);\r
+\r
+  Double_t bfield[3];\r
+  param1->GetBxByBz(bfield);\r
+  bfield[0]=0;\r
+  bfield[1]=0;\r
+  bfield[2]=5.;\r
+  Double_t bz = 5.; // = param1->GetBz();\r
+  //cout<<"In ComputeSignDCA, bfield[3] and bz: "<<bfield[0]<<" "<<bfield[1]<<" "<<bfield[2]<<" "<<bz<<endl;\r
+  Double_t xplane1 = 0.; Double_t xplane2 = 0.;\r
+  Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2);\r
+\r
+  param1->PropagateToBxByBz(xplane1,bfield);\r
+  param2->PropagateToBxByBz(xplane2,bfield);\r
+\r
+  Int_t id1 = 0, id2 = 0;\r
+  AliESDv0 bvertex(*param1,id1,*param2,id2);\r
+  Double_t vx,vy,vz;\r
+  bvertex.GetXYZ(vx,vy,vz);\r
+\r
+  Double_t emom[3];\r
+  Double_t hmom[3];\r
+  param1->PxPyPz(emom);\r
+  param2->PxPyPz(hmom);\r
+  TVector3 emomAtB(emom[0],emom[1],emom[2]);\r
+  TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);\r
+  TVector3 secvtxpt(vx,vy,vz);\r
+  TVector3 decayvector(0,0,0);\r
+  decayvector = secvtxpt - primV; //decay vector from PrimVtx\r
+  Double_t decaylength = decayvector.Mag();\r
+\r
+  //printf("\t JLK pairDCA = %2.2f\n",pairdca);\r
+\r
+  if(GetDebug() > 0) {\r
+    printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );\r
+    printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );\r
+  }\r
+\r
+\r
+\r
+  if (emomAtB.Mag()>0 && pairdca < pdcacut && decaylength < fDecayLenCut ) {\r
+    TVector3 sumMom = emomAtB+hmomAtB;\r
+    Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);\r
+    Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);\r
+    Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);\r
+    Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
+    Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));\r
+    Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();\r
+\r
+//     if (masscut<1.1) fhDVMBtagQA2->Fill(sDca, mass);\r
+\r
+    if (mass > masscut && massPhot > 0.1) signDca = sDca;\r
+    \r
+    if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);\r
+    if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);\r
+  }\r
+\r
+  //clean up\r
+  delete param1;\r
+  delete param2;\r
+\r
+  return signDca;\r
+}\r
+\r
+\r
+//__________________________________________________________________\r
+//PhotonicPrim() removed, it's shit.\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaBtag::PhotonicV0(Int_t id) \r
+{\r
+  //This method checks to see whether a track that has been flagged as\r
+  //an electron was determined to match to a V0 candidate with\r
+  //invariant mass consistent with photon conversion\r
+\r
+  Bool_t itIS = kFALSE;\r
+\r
+  Double_t massEta = 0.547;\r
+  Double_t massRho0 = 0.770;\r
+  Double_t massOmega = 0.782;\r
+  Double_t massPhi = 1.020;\r
+  \r
+  //---Get V0s---\r
+  AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();\r
+  int nv0s = aod->GetNumberOfV0s();\r
+  for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
+    AliAODv0 *v0 = aod->GetV0(iV0);\r
+    if (!v0) continue;\r
+    double radius = v0->RadiusV0();\r
+    double mass = v0->InvMass2Prongs(0,1,11,11);\r
+    if(GetDebug() > 0) {\r
+      printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );\r
+      printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);\r
+      printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );\r
+    }\r
+    if (mass < 0.100 ||\r
+       (mass > massEta-0.05 || mass < massEta+0.05) ||\r
+       (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
+       (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
+       (mass > massPhi-0.05 || mass < massPhi+0.05)) {\r
+      if ( id == v0->GetNegID() || id == v0->GetPosID()) {\r
+       itIS=kTRUE;\r
+       if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );\r
+      }\r
+    } }\r
+  return itIS;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t AliAnaBtag::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3]) \r
+{\r
+  //Use the Event vertex and AOD track information to get\r
+  //a real impact parameter for the track\r
+\r
+\r
+  Double_t maxD = 100000.; //max transverse IP\r
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
+    AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();\r
+    AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();\r
+    AliESDtrack esdTrack(track);\r
+    Double_t bfield[3];\r
+    esdTrack.GetBxByBz(bfield);\r
+    bfield[0]=0;\r
+    bfield[1]=0;\r
+    bfield[2]=5.;\r
+    \r
+    Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov);\r
+    return gotit;\r
+  }\r
+\r
+  return kFALSE;\r
+\r
+}\r
+//__________________________________________________________________\r
+Bool_t AliAnaBtag::CheckIfBjet(const AliAODTrack* track)\r
+{\r
+  Bool_t bjet = kFALSE;\r
+  Int_t trackId = track->GetID(); //get the index in the reader\r
+  Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
+  if(GetDebug() > 3) printf("AliAnaBtag::CheckIfBjet() - aod branch entries %d\n", naod);\r
+\r
+  for(Int_t iaod = 0; iaod < naod ; iaod++){\r
+    AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
+    Int_t electronlabel = ele->GetTrackLabel(0);\r
+    if(electronlabel != trackId) continue;  //skip to the next one if they don't match\r
+    if(ele->GetBtag()>0)\r
+      bjet = kTRUE;\r
+  } \r
+\r
+  return bjet;\r
+} \r
+\r
+\r
+\r
+//__________________________________________________________________\r
+AliAODMCParticle* AliAnaBtag::GetMCParticle(Int_t ipart) \r
+{\r
+  //Get the MC particle at position ipart\r
+\r
+  AliAODMCParticle* aodprimary = 0x0;\r
+  TClonesArray * mcparticles0 = 0x0;\r
+  TClonesArray * mcparticles1 = 0x0;\r
+\r
+  if(GetReader()->ReadAODMCParticles()){\r
+    //Get the list of MC particles                                                                                                                           \r
+    mcparticles0 = GetReader()->GetAODMCParticles(0);\r
+    if(!mcparticles0 && GetDebug() > 0) {\r
+      printf("AliAnaBtag::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
+    }\r
+    if(GetReader()->GetSecondInputAODTree()){\r
+      mcparticles1 = GetReader()->GetAODMCParticles(1);\r
+      if(!mcparticles1 && GetDebug() > 0) {\r
+       printf("AliAnaBtag::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");\r
+      }\r
+    }\r
+\r
+    Int_t npart0 = mcparticles0->GetEntriesFast();\r
+    Int_t npart1 = 0;\r
+    if(mcparticles1) npart1 = mcparticles1->GetEntriesFast();\r
+    if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
+    else aodprimary = (AliAODMCParticle*)mcparticles1->At(ipart-npart0);\r
+    if(!aodprimary) {\r
+      printf("AliAnaBtag::GetMCParticle() *** no primary ***:  label %d \n", ipart);\r
+      return 0x0;\r
+    }\r
+\r
+  } else {\r
+    printf("AliAnaBtag::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
+  }\r
+  return aodprimary;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t  AliAnaBtag::IsMcBJet(Double_t jeta, Double_t jphi)\r
+{\r
+  //Check the jet eta,phi against that of the b-quark\r
+  //to decide whether it is an MC B-jet\r
+  Bool_t bjet=kFALSE;\r
+\r
+  //      printf("MTH: McStack ,nparticles=%d \n", stack->GetNtrack() );\r
+\r
+  AliStack* stack = 0x0;\r
+  \r
+  for(Int_t ipart = 0; ipart < 100; ipart++) {\r
+\r
+    Double_t pphi = -999.;\r
+    Double_t peta = -999.;\r
+    Int_t pdg = 0;\r
+    if(GetReader()->ReadStack()) {\r
+      stack = GetMCStack();\r
+      if(!stack) {\r
+       printf("AliAnaBtag::IsMCBJet() *** no stack ***: \n");\r
+       return kFALSE;\r
+      }\r
+      TParticle* primary = stack->Particle(ipart);\r
+      if (!primary) continue;\r
+      pdg = primary->GetPdgCode();\r
+      pphi = primary->Phi();\r
+      peta = primary->Eta();\r
+    } else if(GetReader()->ReadAODMCParticles()) {\r
+      AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
+      if(!aodprimary) continue;\r
+      pdg = aodprimary->GetPdgCode();\r
+      pphi = aodprimary->Phi();\r
+      peta = aodprimary->Eta();\r
+    }\r
+    if ( TMath::Abs(pdg) != 5) continue;\r
+    \r
+    //      printf("MTH: IsMcBJet : %d, pdg=%d : pt=%f \n", ipart, pdgcode, primary->Pt());\r
+    Double_t dphi = jphi - pphi;\r
+    Double_t deta = jeta - peta;\r
+    Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+    \r
+    if (dr < 0.2) {\r
+      bjet=kTRUE;\r
+      //printf("MTH: **** found matching MC-Bjet: PDG=%d, pt=%f,dr=%f \n", pdgcode, primary->Pt(),dr );\r
+      break;\r
+    }\r
+  }\r
+  return bjet;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+Bool_t  AliAnaBtag::IsMcDJet(Double_t jeta, Double_t jphi)\r
+{\r
+  //Check if this jet is a charm jet\r
+  Bool_t cjet=kFALSE;\r
+\r
+  AliStack* stack = 0x0;\r
+\r
+  for(Int_t ipart = 0; ipart < 100; ipart++) {\r
+    \r
+    Double_t pphi = -999.;\r
+    Double_t peta = -999.;\r
+    Int_t pdg = 0;\r
+    if(GetReader()->ReadStack()) {\r
+      stack = GetMCStack();\r
+      if(!stack) {\r
+       printf("AliAnaBtag::IsMCDJet() *** no stack ***: \n");\r
+       return kFALSE;\r
+      }\r
+      TParticle* primary = stack->Particle(ipart);\r
+      if (!primary) continue;\r
+      pdg = primary->GetPdgCode();\r
+      pphi = primary->Phi();\r
+      peta = primary->Eta();\r
+    } else if(GetReader()->ReadAODMCParticles()) {\r
+      AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
+      if(!aodprimary) continue;\r
+      pdg = aodprimary->GetPdgCode();\r
+      pphi = aodprimary->Phi();\r
+      peta = aodprimary->Eta();\r
+    }\r
+\r
+    if ( TMath::Abs(pdg) != 4) continue;\r
+\r
+    Double_t dphi = jphi - pphi;\r
+    Double_t deta = jeta - peta;\r
+    Double_t dr = sqrt(deta*deta + dphi*dphi);\r
+    \r
+    if (dr < 0.2) {\r
+      cjet=kTRUE;\r
+      break;\r
+    }\r
+  }\r
+\r
+  return cjet;\r
+\r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaBtag::Print(const Option_t * opt) const\r
+{\r
+  //Print some relevant parameters set for the analysis\r
+  \r
+  if(! opt)\r
+    return;\r
+  \r
+  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
+  AliAnaPartCorrBaseClass::Print(" ");\r
+\r
+  printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;\r
+  printf("pOverE range           =     %f - %f\n",fpOverEmin,fpOverEmax);\r
+  printf("residual cut           =     %f\n",fResidualCut);\r
+  printf("---DVM Btagging\n");\r
+  printf("max IP-cut (e,h)       =     %f\n",fImpactCut);\r
+  printf("min ITS-hits           =     %d\n",fITSCut);\r
+  printf("max dR (e,h)           =     %f\n",fDrCut);\r
+  printf("max pairDCA            =     %f\n",fPairDcaCut);\r
+  printf("max decaylength        =     %f\n",fDecayLenCut);\r
+  printf("min Associated Pt      =     %f\n",fAssocPtCut);\r
+  printf("---IPSig Btagging\n");\r
+  printf("min tag track          =     %d\n",fNTagTrkCut);\r
+  printf("min IP significance    =     %f\n",fIPSigCut);\r
+  printf("    \n") ;\r
+       \r
+} \r
+\r
+\r
+//__________________________________________________________________\r
+void  AliAnaBtag::Terminate(TList* outputList)\r
+{\r
\r
+  //Do some plots to end\r
+  //Recover histograms from output histograms list, needed for\r
+  //distributed analysis.                \r
+  //ReadHistograms(outputList);\r
+\r
+  printf(" AliAnaBtag::Terminate()  *** %s Report: %d outputs\n", GetName(), outputList->GetEntries()) ;\r
+\r
+}\r
+\r
diff --git a/PWG4/PartCorrDep/AliAnaBtag.h b/PWG4/PartCorrDep/AliAnaBtag.h
new file mode 100644 (file)
index 0000000..8ffce1b
--- /dev/null
@@ -0,0 +1,170 @@
+#ifndef ALIANABTAG_H\r
+#define ALIANABTAG_H\r
+/* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice     */\r
+/* $Id:  $ */\r
+\r
+//_________________________________________________________________________\r
+//\r
+// Class for the electron identification and B-tagging.\r
+// Clusters from EMCAL matched to tracks are selected \r
+// and kept in the AOD. Few histograms produced.\r
+//\r
+\r
+//-- Author T.R.P-R.Aronsson, J.Klay\r
+\r
+// --- ROOT system ---\r
+class TH2F ;\r
+class TString ;\r
+class TNtuple ;\r
+class TH3F;\r
+\r
+// --- ANALYSIS system ---\r
+#include "AliAnaPartCorrBaseClass.h"\r
+\r
+class AliAODMCParticle;\r
+class AliCaloTrackReader;\r
+class AliAODTrack;\r
+class TList ;\r
+\r
+class AliAnaBtag : public AliAnaPartCorrBaseClass {\r
+\r
+public: \r
+  AliAnaBtag() ; // default ctor\r
+  virtual ~AliAnaBtag() ; //virtual dtor\r
+private:\r
+       AliAnaBtag(const AliAnaBtag & g) ; // cpy ctor\r
+       AliAnaBtag & operator = (const AliAnaBtag & g) ;//cpy assignment\r
+       \r
+public:\r
+       \r
+  TList *  GetCreateOutputObjects();\r
+\r
+  //Main functions\r
+  void Init();\r
+  void MakeAnalysisFillAOD()  ;\r
+  void MakeAnalysisFillHistograms() ; \r
+  \r
+  //B-tagging\r
+  Int_t GetDVMBtag(AliAODTrack * tr); //returns # tracks from secvtx\r
+\r
+  //Temporary local method to get DCA because AliAODTrack is stupid\r
+  Bool_t GetDCA(const AliAODTrack* tr,Double_t imp[2], Double_t cov[3]);\r
+  Bool_t PhotonicV0(Int_t trackId); //check with V0 list\r
+\r
+  //check if track has been flagged as a non-photonic or DVM electron\r
+  //used with the jet tracks to tag bjets\r
+  Bool_t CheckIfBjet(const AliAODTrack* track);\r
+  Bool_t IsMcBJet(Double_t x, Double_t y);\r
+  Bool_t IsMcDJet(Double_t x, Double_t y);\r
+\r
+  void Print(const Option_t * opt)const;\r
+  \r
+  TString GetCalorimeter()   const {return fCalorimeter ; }\r
+  Double_t GetpOverEmin()   const {return fpOverEmin ; }\r
+  Double_t GetpOverEmax()   const {return fpOverEmax ; }\r
+  Bool_t GetWriteNtuple()   const {return fWriteNtuple ; }\r
+\r
+  Double_t GetDrCut() const { return fDrCut; }\r
+  Double_t GetPairDcaCut() const { return fPairDcaCut; }\r
+  Double_t GetDecayLenCut() const { return fDecayLenCut; }\r
+  Double_t GetImpactCut() const { return fImpactCut; }\r
+  Double_t GetAssocPtCut() const { return fAssocPtCut; }\r
+  Double_t GetMassCut() const { return fMassCut; }\r
+  Double_t GetSdcaCut() const { return fSdcaCut; }\r
+  Int_t    GetITSCut() const { return fITSCut; }\r
+  Int_t    GetNTagTrackCut() const { return fNTagTrkCut; }\r
+  Double_t GetIPSigCut() const { return fIPSigCut; }\r
+  Double_t GetMinClusEne() const { return fMinClusEne; }\r
+\r
+  void SetCalorimeter(TString det)    {fCalorimeter = det ; }\r
+  void SetpOverEmin(Double_t min)     {fpOverEmin = min ; }\r
+  void SetpOverEmax(Double_t max)     {fpOverEmax = max ; }\r
+  void SetResidualCut(Double_t cut)     {fResidualCut = cut ; }\r
+  void SetWriteNtuple(Bool_t val)     {fWriteNtuple = val ; }\r
+\r
+  void SetDrCut(Double_t dr)  { fDrCut = dr; }\r
+  void SetPairDcaCut(Double_t pdca) { fPairDcaCut = pdca; }\r
+  void SetDecayLenCut(Double_t dlen) { fDecayLenCut = dlen; }\r
+  void SetImpactCut(Double_t imp) { fImpactCut = imp; }\r
+  void SetAssocPtCut(Double_t pt) { fAssocPtCut = pt; }\r
+  void SetMassCut(Double_t mass) { fMassCut = mass; }\r
+  void SetSdcaCut(Double_t sdca) { fSdcaCut = sdca; }\r
+  void SetITSCut(Int_t its) { fITSCut = its; }\r
+  void SetNTagTrackCut(Int_t ntr) { fNTagTrkCut = ntr; }\r
+  void SetIPSigCut(Double_t ips) { fIPSigCut = ips; }\r
+  void SetMinClusEne(Double_t ene) { fMinClusEne = ene; }\r
+\r
+  void InitParameters();\r
+  void Terminate(TList * outputList);\r
+         \r
+  private:\r
+  //For DVM B-tag method\r
+  Double_t ComputeSignDca(AliAODTrack *track, AliAODTrack *track2, float cut1, double pdcacut);\r
+\r
+  //Int_t GetMCSource(Int_t mctag);\r
+  AliAODMCParticle* GetMCParticle(Int_t part);\r
+\r
+\r
+  private:\r
+  TString  fCalorimeter;  //! Which detector? EMCAL or PHOS\r
+  Double_t fpOverEmin;    //! Minimum p/E value for Electrons\r
+  Double_t fpOverEmax;    //! Maximum p/E value for Electrons\r
+  Double_t fResidualCut;  //! Track-cluster matching distance\r
+  Double_t fMinClusEne;   //! Min clus energy for matching\r
+\r
+  //DVM B-tagging\r
+  Double_t fDrCut;       //max dR\r
+  Double_t fPairDcaCut;  //max pair-DCA\r
+  Double_t fDecayLenCut; //max 3d-decaylength\r
+  Double_t fImpactCut;   //max track impact param\r
+  Double_t fAssocPtCut;  //min associated pt\r
+  Double_t fMassCut;     //min Minv cut\r
+  Double_t fSdcaCut;     //min signDca\r
+  Int_t   fITSCut;       //min ITS hits (both)\r
+  //IP Sig B-tagging\r
+  Int_t    fNTagTrkCut;  //min number of tracks required for IP sig tag\r
+  Double_t fIPSigCut;    //min IP significance cut\r
+\r
+  Double_t fJetEtaCut;   //max eta for jets\r
+  Double_t fJetPhiMin;   //min phi for jets\r
+  Double_t fJetPhiMax;   //max phi for jets\r
+\r
+  Bool_t  fWriteNtuple; //flag for filling ntuple or not\r
+  \r
+\r
+  ///////////////////////////////////////\r
+  //Output histograms and Ntuples\r
+\r
+  TNtuple * fEleNtuple;                 //Ntuple for electrons if neede\r
+\r
+\r
+  TH1F * fhEmcalElectrons;              //All electrons, as id:ed by EMCAL\r
+  TH1F * fhTRDElectrons;                //Electrons from TRD\r
+  TH1F * fhTPCElectrons;                //Electrons from TPC\r
+  TH1F * fhDVM1;                        //initial b-tag dvm1\r
+  TH1F * fhDVM2;                        //initial b-tag dvm2 \r
+  TH2F * fhJets;                        //All jets 2d\r
+  TH2F * fhJetsAllEtaPhi;               //Eta phi for all jets\r
+  TH2F * fhJetsLeadingBElectronEtaPhi;  //deta dphi btw leading jet and bele\r
+  TH2F * fhDVM1EtaPhi;                  //eta phi for b-electrons\r
+  TH2F * fhBJetElectronDetaDphi;        //eta phi for jet with b-electrons\r
+  TH1F * fhClusterEnergy;               //cluster E of EMCAL\r
+  TH1F * fhTestalle;\r
+  TH1F * fhResidual;                    //Residuals from trackmatching\r
+\r
+  //Analysis of electrons\r
+  TH2F * fhElectrons;\r
+  //Analysis for tracks from esd/aod\r
+  TH2F * fhTracks;\r
+\r
+\r
+  ClassDef(AliAnaBtag,12)\r
+\r
+} ;\r
\r
+\r
+#endif//ALIANABTAG_H\r
+\r
+\r
+\r
index c9e523b97252a0dee19e39798599d2ac62e1583a..01db9a7e42e612b8fa27d14cf08114762e0a60f8 100755 (executable)
@@ -10,6 +10,7 @@ SRCS = PartCorrDep/AliAnaCaloTrigger.cxx \
        PartCorrDep/AliAnaPhoton.cxx PartCorrDep/AliAnaPi0.cxx  \
        PartCorrDep/AliAnaPi0EbE.cxx PartCorrDep/AliAnaChargedParticles.cxx \
        PartCorrDep/AliAnaCalorimeterQA.cxx PartCorrDep/AliAnaElectron.cxx \
+       PartCorrDep/AliAnaBtag.cxx \
        PartCorrDep/AliAnalysisTaskTaggedPhotons.cxx PartCorrDep/AliAnaOmegaToPi0Gamma.cxx 
 
 HDRS:= $(SRCS:.cxx=.h) 
index db7926828d9f55d08356e5c2501ead23b2524327..4b46abac20f552981283e21d8ae042f20beb3183 100644 (file)
@@ -136,7 +136,7 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   
   
   // -------------------------------------------------
-  // --- Photon Isolation and Correlation Analysis ---
+  // --- Photon/Pi0/Omega/Electron Analysis ---
   // -------------------------------------------------
   
   AliAnaPhoton *anaphoton = new AliAnaPhoton();
@@ -201,9 +201,97 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   anapi0->SetHistoMassRangeAndNBins(0., 0.6, 200) ;
   anapi0->SetHistoAsymmetryRangeAndNBins(0., 1. , 10) ;
   if(kPrintSettings) anapi0->Print("");
-  
+       
+  //---------------------------  
+  //Pi0, event by event
+  //---------------------------  
+
+  AliNeutralMesonSelection *nms = new AliNeutralMesonSelection();
+  nms->SetInvMassCutRange(0.05, 0.2)     ;
+  nms->KeepNeutralMesonSelectionHistos(kTRUE);
+  //Set Histrograms bins and ranges
+  nms->SetHistoERangeAndNBins(0, 50, 200) ;
+  //      nms->SetHistoPtRangeAndNBins(0, 50, 100) ;
+  //      nms->SetHistoAngleRangeAndNBins(0, 0.3, 100) ;
+  //      nsm->SetHistoIMRangeAndNBins(0, 0.4, 100) ;  
+               
+  AliAnaPi0EbE *anapi0ebe = new AliAnaPi0EbE();
+  anapi0ebe->SetDebug(-1);//10 for lots of messages
+  anapi0ebe->SetAnalysisType(AliAnaPi0EbE::kIMCalo);
+  anapi0ebe->SetMinPt(0);
+  anapi0ebe->SetCalorimeter(calorimeter);
+  anapi0ebe->SetInputAODName(Form("Photons%s",calorimeter.Data()));
+  if(!data.Contains("delta")) {
+       anapi0ebe->SetOutputAODName(Form("Pi0s%s",calorimeter.Data()));
+       anapi0ebe->SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
+  }
+  else  anapi0ebe->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
+       
+  if(kUseKinematics) anapi0ebe->SwitchOnDataMC() ;//Access MC stack and fill more histograms
+  else  anapi0ebe->SwitchOffDataMC() ; 
+  anapi0ebe->SetNeutralMesonSelection(nms);
+  //Set Histrograms bins and ranges
+  anapi0ebe->SetHistoPtRangeAndNBins(0, 50, 200) ;
+  //      anapi0->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;
+  //      anapi0->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
+  if(kPrintSettings) anapi0ebe->Print("");
+       
+  //-------------------------------------
+  //*** analysis the omega->pi0+gamma ***
+  //------------------------------------
+  AliAnaOmegaToPi0Gamma *anaomegaToPi0Gamma = new AliAnaOmegaToPi0Gamma();
+  anaomegaToPi0Gamma->SetDebug(-1);//10 for lots of messages
+  anaomegaToPi0Gamma->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
+  anaomegaToPi0Gamma->SetInputAODPhotonName(Form("Photons%s",calorimeter.Data()));
+  anaomegaToPi0Gamma->SetNPID(1);
+  anaomegaToPi0Gamma->SetNVtxZ(1);
+  anaomegaToPi0Gamma->SetNEventsMixed(4);
+  if(calorimeter=="PHOS")
+       anaomegaToPi0Gamma->SetPi0MassPeakWidthCut(0.008); // PHOS
+  else if(calorimeter=="EMCAL")
+       anaomegaToPi0Gamma->SetPi0MassPeakWidthCut(0.012); // EMCAL 
+  anaomegaToPi0Gamma->SetHistoPtRangeAndNBins(0, 20, 100) ;
+  anaomegaToPi0Gamma->SetHistoMassRangeAndNBins(0, 1, 100) ;
+  anaomegaToPi0Gamma->SetPi0OverOmegaPtCut(0.8);
+  anaomegaToPi0Gamma->SetGammaOverOmegaPtCut(0.2);
+  if(kUseKinematics) anaomegaToPi0Gamma->SwitchOnDataMC() ;//Access MC stack and fill more histograms
+  else anaomegaToPi0Gamma->SwitchOffDataMC() ;//Access MC stack and fill more histograms
+  anaomegaToPi0Gamma->AddToHistogramsName(Form("AnaOmegaToPi0Gamma%s_",calorimeter.Data()));
+  if(kPrintSettings)   anaomegaToPi0Gamma->Print("");
+       
+       
+  //---------------------------------------------------------------------
+  // Electron/btag
+  //---------------------------------------------------------------------
+  if(calorimeter=="EMCAL"){
+       Double_t pOverEmin = 0.8;  //tight
+       Double_t pOverEmax = 1.2;  //tight
+       Double_t dRmax     = 0.02; //tight
+       
+       AliAnaBtag *anaelectron = new AliAnaBtag();
+       anaelectron->SetDebug(-1); //10 for lots of messages
+       anaelectron->SetCalorimeter("EMCAL");
+       if(kUseKinematics){
+               anaelectron->SwitchOffDataMC();
+               anaelectron->SetMinPt(1.);
+       }
+       anaelectron->SetOutputAODName("Electrons");
+       anaelectron->SetOutputAODClassName("AliAODPWG4Particle");
+       anaelectron->SetWriteNtuple(kFALSE);
+       //Determine which cuts to use based on enum
+       anaelectron->SetpOverEmin(pOverEmin);
+       anaelectron->SetpOverEmax(pOverEmax);
+       anaelectron->SetResidualCut(dRmax);
+       //Set Histrograms bins and ranges 
+       anaelectron->SetHistoPtRangeAndNBins(0, 100, 100) ;
+       anaelectron->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;
+       anaelectron->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
+       if(kPrintSettings)anaelectron->Print("");
+  }
+  //==================================
   // ### Isolation analysis ###        
-  
+  //=================================
+  //Photon
   AliAnaParticleIsolation *anaisol = new AliAnaParticleIsolation();
   anaisol->SetDebug(-1);
   anaisol->SetMinPt(0);
@@ -234,6 +322,39 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   anaisol->AddToHistogramsName("AnaIsolPhoton_");
   if(kPrintSettings) anaisol->Print("");
   
+  //Pi0
+  AliAnaParticleIsolation *anaisolpi0 = new AliAnaParticleIsolation();
+  anaisolpi0->SetDebug(-1);
+  anaisolpi0->SetMinPt(0);
+  anaisolpi0->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
+  anaisolpi0->AddToHistogramsName("AnaIsolPi0_");
+  anaisolpi0->SetCalorimeter(calorimeter);
+  if(kUseKinematics) anaisolpi0->SwitchOnDataMC() ;//Access MC stack and fill more histograms
+  else  anaisolpi0->SwitchOffDataMC() ;
+  //Select clusters with no pair, if both clusters with pi0 mass
+  anaisolpi0->SwitchOffInvariantMass();
+  //anaisol->SetNeutralMesonSelection(nms);
+  //Do isolation cut
+  AliIsolationCut * ic2 =  anaisolpi0->GetIsolationCut();      
+  ic2->SetConeSize(0.4);
+  ic2->SetPtThreshold(0.2);
+  ic2->SetICMethod(AliIsolationCut::kPtThresIC);
+  if(kPrintSettings) ic->Print("");
+  //Do or not do isolation with previously produced AODs.
+  //No effect if use of SwitchOnSeveralIsolation()
+  anaisolpi0->SwitchOffReIsolation();
+  //Multiple IC
+  anaisolpi0->SwitchOffSeveralIsolation() ;
+  //Set Histograms bins and ranges
+  anaisolpi0->SetHistoPtRangeAndNBins(0, 50, 200) ;
+  //      ana->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;
+  //      ana->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
+  if(kPrintSettings) anaisol->Print("");
+       
+  //===========================
+  //Correlation analysis
+  //===========================
+       
   // ### Correlation with Jet Finder AOD output
   AliAnaParticleJetFinderCorrelation *anacorrjet = new AliAnaParticleJetFinderCorrelation();
   anacorrjet->SetInputAODName(Form("Photons%s",calorimeter.Data()));
@@ -301,68 +422,6 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   //      ana->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
   if(kPrintSettings) anacorrisohadron->Print("");
   
-  // -------------------------------------------------
-  // --- Pi0    Isolation and Correlation Analysis ---
-  // -------------------------------------------------
-       
-  AliNeutralMesonSelection *nms = new AliNeutralMesonSelection();
-  nms->SetInvMassCutRange(0.05, 0.2)     ;
-  nms->KeepNeutralMesonSelectionHistos(kTRUE);
-  //Set Histrograms bins and ranges
-  nms->SetHistoERangeAndNBins(0, 50, 200) ;
-  //      nms->SetHistoPtRangeAndNBins(0, 50, 100) ;
-  //      nms->SetHistoAngleRangeAndNBins(0, 0.3, 100) ;
-  //      nsm->SetHistoIMRangeAndNBins(0, 0.4, 100) ;  
-  
-  AliAnaPi0EbE *anapi0ebe = new AliAnaPi0EbE();
-  anapi0ebe->SetDebug(-1);//10 for lots of messages
-  anapi0ebe->SetAnalysisType(AliAnaPi0EbE::kIMCalo);
-  anapi0ebe->SetMinPt(0);
-  anapi0ebe->SetCalorimeter(calorimeter);
-  anapi0ebe->SetInputAODName(Form("Photons%s",calorimeter.Data()));
-  if(!data.Contains("delta")) {
-    anapi0ebe->SetOutputAODName(Form("Pi0s%s",calorimeter.Data()));
-    anapi0ebe->SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
-  }
-  else  anapi0ebe->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
-  
-  if(kUseKinematics) anapi0ebe->SwitchOnDataMC() ;//Access MC stack and fill more histograms
-  else  anapi0ebe->SwitchOffDataMC() ; 
-  anapi0ebe->SetNeutralMesonSelection(nms);
-  //Set Histrograms bins and ranges
-  anapi0ebe->SetHistoPtRangeAndNBins(0, 50, 200) ;
-  //      anapi0->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;
-  //      anapi0->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
-  if(kPrintSettings) anapi0ebe->Print("");
-  
-  AliAnaParticleIsolation *anaisolpi0 = new AliAnaParticleIsolation();
-  anaisolpi0->SetDebug(-1);
-  anaisolpi0->SetMinPt(0);
-  anaisolpi0->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
-  anaisolpi0->AddToHistogramsName("AnaIsolPi0_");
-  anaisolpi0->SetCalorimeter(calorimeter);
-  if(kUseKinematics) anaisolpi0->SwitchOnDataMC() ;//Access MC stack and fill more histograms
-  else  anaisolpi0->SwitchOffDataMC() ;
-  //Select clusters with no pair, if both clusters with pi0 mass
-  anaisolpi0->SwitchOffInvariantMass();
-  //anaisol->SetNeutralMesonSelection(nms);
-  //Do isolation cut
-  AliIsolationCut * ic2 =  anaisolpi0->GetIsolationCut();      
-  ic2->SetConeSize(0.4);
-  ic2->SetPtThreshold(0.2);
-  ic2->SetICMethod(AliIsolationCut::kPtThresIC);
-  if(kPrintSettings) ic->Print("");
-  //Do or not do isolation with previously produced AODs.
-  //No effect if use of SwitchOnSeveralIsolation()
-  anaisolpi0->SwitchOffReIsolation();
-  //Multiple IC
-  anaisolpi0->SwitchOffSeveralIsolation() ;
-  //Set Histograms bins and ranges
-  anaisolpi0->SetHistoPtRangeAndNBins(0, 50, 200) ;
-  //      ana->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;
-  //      ana->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
-  if(kPrintSettings) anaisol->Print("");
-  
   
   // ### Pi0 Correlation with hadrons, not isolated
   AliAnaParticleHadronCorrelation *anacorrhadronpi0 = new AliAnaParticleHadronCorrelation();
@@ -414,26 +473,6 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   //      ana->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;
   if(kPrintSettings) anacorrhadronisopi0->Print("");
   
-  //*** analysis the omega->pi0+gamma ***
-  AliAnaOmegaToPi0Gamma *anaomegaToPi0Gamma = new AliAnaOmegaToPi0Gamma();
-  anaomegaToPi0Gamma->SetDebug(-1);//10 for lots of messages
-  anaomegaToPi0Gamma->SetInputAODName(Form("Pi0s%s",calorimeter.Data()));
-  anaomegaToPi0Gamma->SetInputAODPhotonName(Form("Photons%s",calorimeter.Data()));
-  anaomegaToPi0Gamma->SetNPID(1);
-  anaomegaToPi0Gamma->SetNVtxZ(1);
-  anaomegaToPi0Gamma->SetNEventsMixed(4);
-  if(calorimeter=="PHOS")
-    anaomegaToPi0Gamma->SetPi0MassPeakWidthCut(0.008); // PHOS
-  else if(calorimeter=="EMCAL")
-    anaomegaToPi0Gamma->SetPi0MassPeakWidthCut(0.012); // EMCAL 
-  anaomegaToPi0Gamma->SetHistoPtRangeAndNBins(0, 20, 100) ;
-  anaomegaToPi0Gamma->SetHistoMassRangeAndNBins(0, 1, 100) ;
-  anaomegaToPi0Gamma->SetPi0OverOmegaPtCut(0.8);
-  anaomegaToPi0Gamma->SetGammaOverOmegaPtCut(0.2);
-  if(kUseKinematics) anaomegaToPi0Gamma->SwitchOnDataMC() ;//Access MC stack and fill more histograms
-  else anaomegaToPi0Gamma->SwitchOffDataMC() ;//Access MC stack and fill more histograms
-  anaomegaToPi0Gamma->AddToHistogramsName(Form("AnaOmegaToPi0Gamma%s_",calorimeter.Data()));
-  if(kPrintSettings)   anaomegaToPi0Gamma->Print("");
   
   // #### Configure Maker ####
   AliAnaPartCorrMaker * maker = new AliAnaPartCorrMaker();
@@ -444,7 +483,8 @@ AliAnalysisTaskParticleCorrelation *AddTaskPartCorr(TString data, TString calori
   maker->AddAnalysis(anaphoton,n++);
   maker->AddAnalysis(anapi0,n++);
   maker->AddAnalysis(anapi0ebe,n++);
-  maker->AddAnalysis(anaomegaToPi0Gamma,n++);   
+  maker->AddAnalysis(anaomegaToPi0Gamma,n++);  
+  if(calorimeter=="EMCAL")maker->AddAnalysis(anaelectron,n++);   
   // Isolation analysis
   maker->AddAnalysis(anaisol,n++);
   maker->AddAnalysis(anacorrisohadron,n++);
diff --git a/PWG4/macros/ConfigAnalysisBTag.C b/PWG4/macros/ConfigAnalysisBTag.C
new file mode 100644 (file)
index 0000000..8d177f2
--- /dev/null
@@ -0,0 +1,179 @@
+//------------------------------------\r
+// Configuration macro example:\r
+//\r
+// Configure EMCal electron analysis\r
+//\r
+// Modified by: K. Read\r
+//\r
+//------------------------------------\r
+\r
+AliAnaPartCorrMaker*  ConfigAnalysis()\r
+{\r
+  //\r
+  // Configuration goes here\r
+  // \r
+  printf("======================== \n");\r
+  printf("ConfigAnalysisElectron() \n");\r
+  printf("======================== \n");\r
+  Bool_t kInputIsESD = kTRUE;     //uncomment for input ESD\r
+//Bool_t kInputIsESD = kFALSE;    //uncomment for input AODs\r
+  Bool_t kFollowsFilter = kTRUE;  //uncomment if follows ESD filter task\r
+//Bool_t kFollowsFilter = kFALSE; //uncomment if no ESD filter task\r
+  Bool_t kMC = kTRUE; //set to kFALSE for data\r
+\r
+  //enum for the different electron cut sets\r
+  //defined for dR and p/E\r
+  //kTight2 is the default standard cuts\r
+  enum kCutTypes {kTight2, kLooseTight, kTightLoose, kLoose2};\r
+  Int_t kCutSet = kTight2;\r
+  Double_t pOverEmin = 0.8;  //tight\r
+  Double_t pOverEmax = 1.2;  //tight\r
+  Double_t dRmax     = 0.02; //tight\r
+  if (gSystem->Getenv("ELECUTSET")){\r
+    kCutSet = atoi(gSystem->Getenv("ELECUTSET"));\r
+  }\r
+  if(kCutSet == kLooseTight) {\r
+    pOverEmin = 0.6;  //loose\r
+    pOverEmax = 1.4;  //loose\r
+    dRmax     = 0.02; //tight\r
+  }\r
+  if(kCutSet == kTightLoose) {\r
+    pOverEmin = 0.8;  //tight\r
+    pOverEmax = 1.2;  //tight\r
+    dRmax     = 0.05; //loose\r
+  }\r
+  if(kCutSet == kLoose2) {\r
+    pOverEmin = 0.6;  //loose\r
+    pOverEmax = 1.4;  //loose\r
+    dRmax     = 0.05; //loose\r
+  }    \r
+\r
+  //Alternatively, select input via anaInputData environment variable.\r
+  if (gSystem->Getenv("anaInputData")){\r
+    TString kInputData = gSystem->Getenv("anaInputData");\r
+    if( kInputData == "AOD" ){\r
+      kInputIsESD = kFALSE;\r
+      kFollowsFilter = kFALSE;\r
+    }\r
+  }\r
+\r
+  //Alternatively, adjust for real data based on kMC value.\r
+  if (gSystem->Getenv("anakMC")){\r
+    kMC = atoi(gSystem->Getenv("anakMC"));\r
+  }\r
+\r
+  //Detector Fiducial Cuts\r
+  AliFiducialCut * fidCut = new AliFiducialCut();\r
+  fidCut->DoCTSFiducialCut(kFALSE) ;\r
+  fidCut->DoEMCALFiducialCut(kFALSE) ;\r
+  fidCut->DoPHOSFiducialCut(kFALSE) ;\r
+\r
+  //fidCut->SetSimpleCTSFiducialCut(0.9,0.,360.);\r
+  //fidCut->SetSimpleEMCALFiducialCut(0.7,80.,190.);\r
+  //fidCut->SetSimplePHOSFiducialCut(0.13,220.,320.);\r
+\r
+  fidCut->Print("");\r
+\r
+  //-----------------------------------------------------------  \r
+  // Reader\r
+  //-----------------------------------------------------------\r
+  if(kInputIsESD && !kFollowsFilter)AliCaloTrackESDReader *reader = new AliCaloTrackESDReader();\r
+  else           AliCaloTrackAODReader *reader = new AliCaloTrackAODReader();\r
+  reader->SetDebug(-1);//10 for lots of messages\r
+\r
+  //Switch on or off the detectors information that you want\r
+  reader->SwitchOnEMCAL();\r
+  reader->SwitchOnCTS();\r
+  //reader->SwitchOffEMCALCells();     \r
+  reader->SwitchOffPHOS();\r
+  //reader->SwitchOffPHOSCells();      \r
+\r
+  //Kine\r
+  if(kMC && !kInputIsESD){\r
+    reader->SwitchOffStack();          // On  by default, remember to SwitchOnMCData() in analysis classes\r
+    reader->SwitchOnAODMCParticles();  // Off by default, remember to SwitchOnMCData() in analysis classes\r
+  }\r
+\r
+  //Select Trigger Class for real data\r
+  if(!kMC) reader->SetFiredTriggerClassName("CINT1B-ABCE-NOPF-ALL");\r
+\r
+  //Min particle pT\r
+  reader->SetCTSPtMin(0.0);   //new\r
+  reader->SetEMCALPtMin(0.0); //new\r
+  if(kFollowsFilter)reader->SetTrackStatus(0);  //to prevent automatic TPC and ITS refit\r
+\r
+  //In case of generating jet events (with PYTHIA), if pt hard bin is small\r
+  //reject events with large difference between ptHard and triggered jet\r
+  //reader->SetPtHardAndJetPtComparison(kTRUE);\r
+\r
+  reader->SetFiducialCut(fidCut);\r
+\r
+  if(!kInputIsESD){\r
+    // Analysis with tracks, select only tracks with\r
+    // following bits\r
+\r
+    //     //We want tracks fitted in the detectors:\r
+    //     ULong_t status=AliAODTrack::kTPCrefit;\r
+    //     status|=AliAODTrack::kITSrefit;\r
+   \r
+    //     We want tracks whose PID bit is set:\r
+    //     ULong_t status =AliAODTrack::kITSpid;\r
+    //     status|=AliAODTrack::kTPCpid;       \r
+\r
+    // reader->SetTrackStatus(status);\r
+  }\r
+\r
+  reader->Print("");\r
+\r
+\r
+  //Detector Fiducial Cuts\r
+  AliFiducialCut * fidCut2 = new AliFiducialCut();\r
+  fidCut2->DoEMCALFiducialCut(kTRUE) ;\r
+  fidCut2->SetSimpleEMCALFiducialCut(0.7,80.,190.);\r
+\r
+  fidCut2->DoCTSFiducialCut(kTRUE) ;\r
+  fidCut2->SetSimpleCTSFiducialCut(0.9,0.,360.); \r
+\r
+  fidCut2->Print("");\r
+\r
+  //---------------------------------------------------------------------\r
+  // Analysis algorithm\r
+  //---------------------------------------------------------------------\r
+\r
+  AliAnaBtag *anaelectron = new AliAnaBtag();\r
+  anaelectron->SetDebug(-1); //10 for lots of messages\r
+  anaelectron->SetCalorimeter("EMCAL");\r
+  if(kMC){\r
+    anaelectron->SwitchOffDataMC();\r
+    anaelectron->SetMinPt(1.);\r
+  }\r
+  anaelectron->SetOutputAODName("Electrons");\r
+  anaelectron->SetOutputAODClassName("AliAODPWG4Particle");\r
+  anaelectron->SetWriteNtuple(kFALSE);\r
+  //Determine which cuts to use based on enum\r
+  anaelectron->SetpOverEmin(pOverEmin);\r
+  anaelectron->SetpOverEmax(pOverEmax);\r
+  anaelectron->SetResidualCut(dRmax);\r
+  //Set Histrograms bins and ranges\r
+  anaelectron->SetHistoPtRangeAndNBins(0, 100, 100) ;\r
+  anaelectron->SetHistoPhiRangeAndNBins(0, TMath::TwoPi(), 100) ;\r
+  anaelectron->SetHistoEtaRangeAndNBins(-0.7, 0.7, 100) ;\r
+  anaelectron->Print("");\r
+\r
+  //---------------------------------------------------------------------\r
+  // Set  analysis algorithm and reader\r
+  //---------------------------------------------------------------------\r
+  maker = new AliAnaPartCorrMaker();\r
+  maker->SetReader(reader);//pointer to reader\r
+  maker->AddAnalysis(anaelectron,0);\r
+  maker->SetAnaDebug(-1)  ;\r
+  maker->SwitchOnHistogramsMaker()  ;\r
+  maker->SwitchOnAODsMaker()  ;\r
+\r
+  maker->Print("");\r
+  //\r
+  printf("============================ \n");\r
+  printf("END ConfigAnalysisElectron() \n");\r
+  printf("============================ \n");\r
+  return maker ;\r
+}\r