PWG4/PartCorrDep/AliAnaBtag.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4 * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes hereby granted      *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 /* $Id: $ */\r
  16 \r
  17 //_________________________________________________________________________\r
  18 //\r
  19 // Class for the electron identification and B-tagging.\r
  20 // Clusters from EMCAL matched to tracks to id electrons.\r
  21 // Btagger is run on all electrons, then jets are tagged as well.\r
  22 // \r
  23 //\r
  24 // -- Author: T.R.P.Aronsson (Yale), M. Heinz (Yale)\r
  25 //////////////////////////////////////////////////////////////////////////////\r
  26   \r
  27 #include <TH2F.h>\r
  28 #include <TH3F.h>\r
  29 #include <TParticle.h>\r
  30 #include <TNtuple.h>\r
  31 #include <TClonesArray.h>\r
  32 \r
  33 #include "AliAnaBtag.h" \r
  34 #include "AliCaloTrackReader.h"\r
  35 #include "AliMCAnalysisUtils.h"\r
  36 #include "AliVCluster.h"\r
  37 #include "AliFiducialCut.h"\r
  38 #include "AliAODTrack.h"\r
  39 #include "AliAODPid.h"\r
  40 #include "AliCaloPID.h"\r
  41 #include "AliAODMCParticle.h"\r
  42 #include "AliStack.h"\r
  43 #include "AliExternalTrackParam.h"\r
  44 #include "AliESDv0.h"\r
  45 #include "AliESDtrack.h"\r
  46 #include "AliAODJet.h"\r
  47 #include "AliAODEvent.h"\r
  48 #include "AliGenPythiaEventHeader.h"\r
  49 ClassImp(AliAnaBtag)\r
  50   \r
  51 //____________________________________________________________________________\r
  52 AliAnaBtag::AliAnaBtag() \r
  53 : AliAnaPartCorrBaseClass(),\r
  54   fWriteNtuple(0),electrons(0),pairs(0),events(0),fEventNumber(0),fNElec(0),fNElecEv(0),fNPair(0),fDrCut(0.),fPairDcaCut(0.),fDecayLenCut(0.),fImpactCut(0.),\r
  55   fAssocPtCut(0.),fMassCut(0.),fSdcaCut(0.),fITSCut(0),\r
  56   fNTagTrkCut(0),fIPSigCut(0.),fJetEtaCut(0.3),fJetPhiMin(1.8),fJetPhiMax(2.9),\r
  57   fhEmcalElectrons(0),fhTRDElectrons(0),fhTPCElectrons(0),fhEmcalMCE(0),fhTRDMCE(0),fhTPCMCE(0),fhEmcalMCEFake(0),fhTRDMCEFake(0),fhTPCMCEFake(0),fhEmcalMCP(0),fhTRDMCP(0),fhTPCMCP(0),fhEmcalMCK(0),fhTRDMCK(0),fhTPCMCK(0),fhSpecies(0),fhDVM1(0),fhDVM2(0),fhNVTX(0),fhNVTXMC(0),fhJets(0),fhJetsAllEtaPhi(0),fhJetsLeadingBElectronEtaPhi(0),fhDVM1EtaPhi(0),fhBJetElectronDetaDphi(0),fhClusterMap(0),fhClusterEnergy(0),fhTestalle(0),fhResidual(0),fhPairPt(0),fhElectrons(0),fhTracks(0)\r
  58 {\r
  59   //default ctor\r
  60   //Initialize parameters\r
  61   InitParameters();\r
  62 \r
  63 }\r
  64 \r
  65 //____________________________________________________________________________\r
  66 AliAnaBtag::~AliAnaBtag() \r
  67 {\r
  68   //dtor\r
  69 \r
  70 }\r
  71 \r
  72 \r
  73 //________________________________________________________________________\r
  74 TList *  AliAnaBtag::GetCreateOutputObjects()\r
  75 {  \r
  76   // Create histograms to be saved in output file and \r
  77   // store them in outputContainer\r
  78   TList * outputContainer = new TList() ; \r
  79   outputContainer->SetName("ElectronHistos") ; \r
  80 \r
  81   if(fWriteNtuple){\r
  82     electrons = new TNtuple("electrons","Electron Ntuple","electronnumber:eventnumber:electronineventnumber:pairs:start:stop:tpc:trd:emc:mcpdg:parentbit:btag:pt");\r
  83     outputContainer->Add(electrons);\r
  84 \r
  85     pairs = new TNtuple("pairs","Pair Ntuple","pairnumber:electronnumber:eventnumber:pdca:sdca:minv:pth:massphoton:decaylength");\r
  86     outputContainer->Add(pairs);\r
  87 \r
  88     events = new TNtuple("events","Event NTuple","event");\r
  89     outputContainer->Add(events);\r
  90   }\r
  91   \r
  92   fhEmcalElectrons = new TH1F("fhEmcalElectrons","",400,0,400);\r
  93   outputContainer->Add(fhEmcalElectrons);\r
  94   \r
  95   fhTRDElectrons = new TH1F("fhTRDElectrons","",400,0,400);\r
  96   outputContainer->Add(fhTRDElectrons);\r
  97   \r
  98   fhTPCElectrons = new TH1F("fhTPCElectrons","",400,0,400);\r
  99   outputContainer->Add(fhTPCElectrons);\r
 100   \r
 101   fhNVTX = new TH1F("fhNVTX","",20,0,20);\r
 102   outputContainer->Add(fhNVTX);\r
 103 \r
 104   fhDVM1 = new TH1F("fhDVM1","",400,0,400);\r
 105   outputContainer->Add(fhDVM1);\r
 106   \r
 107   fhDVM2 = new TH1F("fhDVM2","",400,0,400);\r
 108   outputContainer->Add(fhDVM2);\r
 109   \r
 110   if(IsDataMC()){\r
 111     fhEmcalMCE = new TH1F("fhEmcalMCE","",400,0,400);\r
 112     outputContainer->Add(fhEmcalMCE);\r
 113     \r
 114     fhTRDMCE = new TH1F("fhTRDMCE","",400,0,400);\r
 115     outputContainer->Add(fhTRDMCE);\r
 116     \r
 117     fhTPCMCE = new TH1F("fhTPCMCE","",400,0,400);\r
 118     outputContainer->Add(fhTPCMCE);\r
 119     \r
 120     fhEmcalMCEFake = new TH1F("fhEmcalMCEFake","",400,0,400);\r
 121     outputContainer->Add(fhEmcalMCEFake);\r
 122     \r
 123     fhTRDMCEFake = new TH1F("fhTRDMCEFake","",400,0,400);\r
 124     outputContainer->Add(fhTRDMCEFake);\r
 125     \r
 126     fhTPCMCEFake = new TH1F("fhTPCMCEFake","",400,0,400);\r
 127     outputContainer->Add(fhTPCMCEFake);\r
 128     \r
 129     fhEmcalMCP = new TH1F("fhEmcalMCP","",400,0,400);\r
 130     outputContainer->Add(fhEmcalMCP);\r
 131     \r
 132     fhTRDMCP = new TH1F("fhTRDMCP","",400,0,400);\r
 133     outputContainer->Add(fhTRDMCP);\r
 134     \r
 135     fhTPCMCP = new TH1F("fhTPCMCP","",400,0,400);\r
 136     outputContainer->Add(fhTPCMCP);\r
 137     \r
 138     fhEmcalMCK = new TH1F("fhEmcalMCK","",400,0,400);\r
 139     outputContainer->Add(fhEmcalMCK);\r
 140     \r
 141     fhTRDMCK = new TH1F("fhTRDMCK","",400,0,400);\r
 142     outputContainer->Add(fhTRDMCK);\r
 143     \r
 144     fhTPCMCK = new TH1F("fhTPCMCK","",400,0,400);\r
 145     outputContainer->Add(fhTPCMCK);\r
 146     \r
 147     fhSpecies  = new TH1F("fhSpecies","",1000,0,1000);\r
 148     outputContainer->Add(fhSpecies);\r
 149     \r
 150     fhNVTXMC = new TH1F("fhNVTXMC","",20,0,20);\r
 151     outputContainer->Add(fhNVTXMC);\r
 152   }\r
 153   \r
 154 \r
 155   fhJets = new TH2F("fhJets","",400,0,400,20,0,20);\r
 156   outputContainer->Add(fhJets);\r
 157   \r
 158   fhJetsAllEtaPhi = new TH2F("fhJetsAllEtaPhi","",100,-2,2,100,-2,8);\r
 159   outputContainer->Add(fhJetsAllEtaPhi);\r
 160   \r
 161   fhJetsLeadingBElectronEtaPhi = new TH2F("fhJetsLeadingBElectronEtaPhi","",100,-5,5,200,-10,10);\r
 162   outputContainer->Add(fhJetsLeadingBElectronEtaPhi);\r
 163   \r
 164   fhDVM1EtaPhi = new TH2F("fhDVM1EtaPhi","",100,-2,2,100,-2,8);\r
 165   outputContainer->Add(fhDVM1EtaPhi);\r
 166   \r
 167   fhBJetElectronDetaDphi = new TH2F("fhBJetElectronDetaDphi","",100,-5,5,200,-10,10);\r
 168   outputContainer->Add(fhBJetElectronDetaDphi);\r
 169   \r
 170   fhClusterMap = new TH2F("fhClusterMap","",100,-2,2,100,-2,8);\r
 171   outputContainer->Add(fhClusterMap);\r
 172   \r
 173   fhClusterEnergy = new TH1F("fhClusterEnergy","",100,0,10);\r
 174   outputContainer->Add(fhClusterEnergy);\r
 175   \r
 176   fhTestalle = new TH1F("fhTestalle","",400,0,400);\r
 177   outputContainer->Add(fhTestalle);\r
 178   \r
 179   fhResidual = new TH1F("fhResidual","",500,0,5);\r
 180   outputContainer->Add(fhResidual);\r
 181 \r
 182   fhPairPt = new TH1F("fhPairPt","",400,0,400);\r
 183   outputContainer->Add(fhPairPt);\r
 184   \r
 185   fhElectrons = new TH2F("fhElectrons","",200,0,100,20,0,20);\r
 186   outputContainer->Add(fhElectrons);\r
 187   \r
 188   fhTracks = new TH2F("fhTracks","",200,0,100,20,0,20);\r
 189   outputContainer->Add(fhTracks);\r
 190     \r
 191   return outputContainer ; \r
 192 }\r
 193 \r
 194 //____________________________________________________________________________\r
 195 void AliAnaBtag::Init()\r
 196 {\r
 197   //do some initialization\r
 198     printf("Rubbish init step AliAnaBtag::Init()");\r
 199 }\r
 200 \r
 201 \r
 202 //____________________________________________________________________________\r
 203 void AliAnaBtag::InitParameters()\r
 204 { \r
 205   //Initialize the parameters of the analysis.\r
 206   SetOutputAODClassName("AliAODPWG4Particle");\r
 207   SetOutputAODName("PWG4Particle");\r
 208   //AddToHistogramsName("");\r
 209   //DVM B-tagging\r
 210   fDrCut       = 1.0; \r
 211   fPairDcaCut  = 0.02;\r
 212   fDecayLenCut = 1.0;\r
 213   fImpactCut   = 0.5;\r
 214   fAssocPtCut  = 1.0;\r
 215   fMassCut     = 1.5;\r
 216   fSdcaCut     = 0.1;\r
 217   fITSCut      = 4;\r
 218   //IPSig B-tagging\r
 219   fNTagTrkCut  = 2;\r
 220   fIPSigCut    = 3.0;\r
 221   //Jet fiducial cuts\r
 222   fJetEtaCut = 0.3;\r
 223   fJetPhiMin = 1.8;\r
 224   fJetPhiMax = 2.9;\r
 225 }\r
 226 \r
 227 //__________________________________________________________________\r
 228 void  AliAnaBtag::MakeAnalysisFillAOD() \r
 229 {\r
 230   fEventNumber++;\r
 231   fNElecEv=0;\r
 232   if(fWriteNtuple)\r
 233     events->Fill(fEventNumber);\r
 234   \r
 235   //This reads in tracks, extrapolates to EMCAL, does p/E selectrons, identifies electron candidates\r
 236   //After candidates are obtained, btagging and saving into AOD.\r
 237   AliStack *stack =0x0;\r
 238   //Double_t bfield = 0.;\r
 239   if(GetDebug()>0) printf("AliAnaBtag::MakeAnalysisFillAOD() - Write ntuple flag is %d \n",fWriteNtuple);\r
 240   //if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) bfield = GetReader()->GetBField();\r
 241   TObjArray *cl = GetEMCALClusters();\r
 242   \r
 243   if(!GetCTSTracks() || GetCTSTracks()->GetEntriesFast() == 0) return ;\r
 244   Int_t ntracks = GetCTSTracks()->GetEntriesFast();\r
 245   if(GetDebug() > 0)\r
 246     printf("AliAnaBtag::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);\r
 247   \r
 248   Int_t iCluster = -999;\r
 249   Int_t ntot = cl->GetEntriesFast();\r
 250   \r
 251   //CLUSTER STUFF \r
 252   for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
 253     AliVCluster * clus = (AliVCluster*) (cl->At(iclus));\r
 254     if(!clus) continue;\r
 255     fhClusterEnergy->Fill(clus->E());\r
 256     Float_t xclus[3];\r
 257     clus->GetPosition(xclus);\r
 258     TVector3 cluspos(xclus[0],xclus[1],xclus[2]);\r
 259     \r
 260     fhClusterMap->Fill(cluspos.Eta(),cluspos.Phi());\r
 261   }\r
 262   \r
 263   \r
 264   \r
 265   for (Int_t itrk =  0; itrk <  ntracks; itrk++) {////////////// track loop\r
 266     iCluster = -999; //start with no match\r
 267     AliAODTrack * track = (AliAODTrack*) (GetCTSTracks()->At(itrk)) ;\r
 268     if(track->GetLabel()<0){\r
 269       if(GetDebug()>0)\r
 270         printf("Negative track label, aborting!\n");\r
 271       continue;\r
 272     }\r
 273     Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
 274     Bool_t dcaOkay = GetDCA(track,imp,cov);  //homegrown dca calculation until AOD is fixed\r
 275     if(!dcaOkay&&GetDebug()>0) printf("AliAnaBtag::FillAOD - Problem computing DCA to primary vertex for track %d.  Skipping it...\n",itrk);\r
 276     fhTracks->Fill(track->Pt(),1);\r
 277     \r
 278     if(track->Pt()<0)\r
 279       continue;\r
 280     \r
 281     AliAODPid* pid = (AliAODPid*) track->GetDetPid();\r
 282     if(pid == 0) {\r
 283       if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD() - No PID object - skipping track %d",itrk);\r
 284       continue;\r
 285     } \r
 286     else {\r
 287       Double_t emcpos[3];\r
 288       pid->GetEMCALPosition(emcpos);\r
 289       Double_t emcmom[3];\r
 290       pid->GetEMCALMomentum(emcmom);\r
 291       \r
 292       TVector3 pos(emcpos[0],emcpos[1],emcpos[2]);\r
 293       TVector3 mom(emcmom[0],emcmom[1],emcmom[2]);\r
 294       Double_t tphi = pos.Phi();\r
 295       Double_t teta = pos.Eta();\r
 296       Double_t tmom = mom.Mag();\r
 297       \r
 298       Bool_t in = kFALSE;\r
 299       if(track->Phi()*180./TMath::Pi() > 80. && track->Phi()*180./TMath::Pi() < 190. &&\r
 300          track->Eta() > -0.7 && track->Eta() < 0.7) in = kTRUE;\r
 301       \r
 302       \r
 303       Double_t dEdx = pid->GetTPCsignal();\r
 304       Int_t pidProb = track->GetMostProbablePID();\r
 305       Bool_t tpcEle = kFALSE; if(dEdx > 70.) tpcEle = kTRUE;\r
 306       Bool_t trkEle = kFALSE; if(pidProb == AliAODTrack::kElectron) trkEle = kTRUE;\r
 307       Bool_t emcEle = kFALSE;   \r
 308       \r
 309       \r
 310       \r
 311       \r
 312       ////////////////////////////////////////////////EMCAL//////////////////////\r
 313       if(mom.Pt() > 1.0 && in) {\r
 314         fhTracks->Fill(track->Pt(),3);\r
 315         Double_t res = 999.;\r
 316         Double_t pOverE = -999.;\r
 317         \r
 318         //Track Matching parameters\r
 319         double minRes=100.;\r
 320         Double_t minR  = 99;\r
 321         Double_t minPe =-1;\r
 322         Double_t minEp =-1;\r
 323         Double_t minMult = -1;\r
 324         Double_t minPt   = -1;\r
 325         \r
 326         for(Int_t iclus = 0; iclus < ntot; iclus++) {\r
 327           AliVCluster * clus = (AliVCluster*) (cl->At(iclus));\r
 328           if(!clus) continue;\r
 329           \r
 330           // new optimized from ben. 2010May\r
 331           if (clus->GetNCells()       < 2    ) continue;\r
 332           if (clus->GetNCells()       > 35   ) continue;\r
 333           if (clus->E()               < 0 ) continue;\r
 334           if (clus->GetDispersion()   > 1.08    ) continue;\r
 335           if (clus->GetM20()          > 0.42  ) continue;\r
 336           if (clus->GetM02()          > 0.4  ) continue;\r
 337           if (clus->GetM20()          < 0 ) continue;\r
 338           if (clus->GetM02()          < 0.06 ) continue;\r
 339           \r
 340           \r
 341           Float_t x[3];\r
 342           clus->GetPosition(x);\r
 343           TVector3 cluspos(x[0],x[1],x[2]);\r
 344           Double_t deta = teta - cluspos.Eta();\r
 345           Double_t dphi = tphi - cluspos.Phi();\r
 346           if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
 347           if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
 348           \r
 349           res = sqrt(dphi*dphi + deta*deta);\r
 350           if(res<minRes)  minRes=res;       \r
 351           \r
 352           if(res < 0.0275) { // { //Optimized from Ben\r
 353             iCluster = iclus;\r
 354             Double_t energy = clus->E(); \r
 355             if(energy > 0) pOverE = tmom/energy;\r
 356             if (res< minR) {\r
 357               minR  = res;\r
 358               minPe = pOverE;\r
 359               minEp = energy/tmom;\r
 360               minMult = clus->GetNCells() ;\r
 361               minPt = track->Pt();\r
 362             }\r
 363           } else {\r
 364             //unmatched\r
 365           }//res cut\r
 366           \r
 367         }//calo cluster loop\r
 368         fhResidual->Fill(minRes);\r
 369         \r
 370         if(minPe > 0.9 && minPe < 1.08) emcEle = kTRUE;//       if(minPe > fpOverEmin && minPe < fpOverEmax) emcEle = kTRUE;\r
 371         \r
 372       }//pt, fiducial selection = EMCAL ////////////////////END EMCAL/////////////////////////\r
 373       \r
 374       \r
 375       \r
 376       \r
 377       ////////////////////////////////////////////////////Electrons/////////////////////     \r
 378       if(emcEle ||tpcEle || trkEle) { //Obsolete (kinda...)\r
 379         fhTestalle->Fill(track->Pt());\r
 380         \r
 381         //B-tagging\r
 382         if(GetDebug() > 1) printf("Found Electron - do b-tagging\n");\r
 383         Int_t pairs1=0,start=0,stop=0;\r
 384         Int_t dvmbtag = GetDVMBtag(track,pairs1,start,stop); //add: get back #pairs, start stop in pair-Ntuple.\r
 385         if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillAOD -  Analyze, got back result from dvm: pair counts: pairs: %d, start %d, stop %d. \n",pairs1,start,stop);\r
 386         fhNVTX->Fill(dvmbtag);\r
 387         \r
 388         if(dvmbtag>0){\r
 389           fhDVM1->Fill(track->Pt());\r
 390           fhDVM1EtaPhi->Fill(track->Eta(),track->Phi());            \r
 391         }\r
 392         if(dvmbtag>1)\r
 393           fhDVM2->Fill(track->Pt());\r
 394         \r
 395         //Create particle to save in AOD///////////// Purpose of this is the AODJets needs to check against this.\r
 396         Double_t eMass = 0.511/1000; //mass in GeV\r
 397         Double_t eleE = sqrt(track->P()*track->P() + eMass*eMass);\r
 398         AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),eleE);\r
 399         tr.SetLabel(track->GetLabel());\r
 400         tr.SetCaloLabel(iCluster,-1); //sets the indices of the original caloclusters\r
 401         tr.SetTrackLabel(track->GetID(),-1); //sets the indices of the original tracks tr.SetTrackLabel(track->GetID(),-1) instead of itrk;\r
 402         tr.SetBtag(dvmbtag);\r
 403         if(track->Charge() < 0) tr.SetPdg(11); //electron is 11\r
 404         else  tr.SetPdg(-11); //positron is -11 \r
 405         \r
 406         //Set detector flags\r
 407         Int_t emcflag=0,tpcflag=0,trdflag=0;\r
 408         if(emcEle){\r
 409           fhEmcalElectrons->Fill(track->Pt());\r
 410           emcflag=1;}\r
 411         if(trkEle){\r
 412           fhTRDElectrons->Fill(track->Pt());\r
 413           trdflag=1;}\r
 414         if(tpcEle){\r
 415           fhTPCElectrons->Fill(track->Pt());\r
 416           tpcflag=1;}\r
 417         \r
 418         if(emcEle) {//PID determined by EMCAL\r
 419           tr.SetDetector("EMCAL");\r
 420         } else {\r
 421           tr.SetDetector("CTS"); //PID determined by CTS\r
 422         }\r
 423         \r
 424         \r
 425         /////////////////////////MC stuff////////////////////////////////////////////////\r
 426         Int_t pdg = 0;\r
 427         if(IsDataMC()){\r
 428           stack=GetMCStack();\r
 429           if(!stack) {printf("AliAnaBtag::MakeAnalysisFillHistograms() - Crap, no stack: \n");\r
 430           }\r
 431           else{\r
 432             //Is it really an electron?\r
 433             TParticle *partX = stack->Particle(TMath::Abs(track->GetLabel()));\r
 434             pdg = partX->GetPdgCode();\r
 435             fhSpecies->Fill(TMath::Abs(pdg));\r
 436             if(TMath::Abs(pdg)==11){ //Check MC electrons\r
 437               if(emcEle)\r
 438                 fhEmcalMCE->Fill(track->Pt());\r
 439               if(trkEle)\r
 440                 fhTRDMCE->Fill(track->Pt());\r
 441               if(tpcEle)\r
 442                 fhTPCMCE->Fill(track->Pt());\r
 443             }else{ //Fake histos!\r
 444               if(emcEle)\r
 445                 fhEmcalMCEFake->Fill(track->Pt());\r
 446               if(trkEle)\r
 447                 fhTRDMCEFake->Fill(track->Pt());\r
 448               if(tpcEle)\r
 449                 fhTPCMCEFake->Fill(track->Pt());\r
 450             }\r
 451             if(TMath::Abs(pdg)==211){ //Check MC pions\r
 452               if(emcEle)\r
 453                 fhEmcalMCP->Fill(track->Pt());\r
 454               if(trkEle)\r
 455                 fhTRDMCP->Fill(track->Pt());\r
 456               if(tpcEle)\r
 457                 fhTPCMCP->Fill(track->Pt());\r
 458             }\r
 459             if(TMath::Abs(pdg)==321){ //Check MC Kaons\r
 460               if(emcEle)\r
 461                 fhEmcalMCK->Fill(track->Pt());\r
 462               if(trkEle)\r
 463                 fhTRDMCK->Fill(track->Pt());\r
 464               if(tpcEle)\r
 465                 fhTPCMCK->Fill(track->Pt());\r
 466             }\r
 467           }\r
 468           \r
 469           //Take care of where it came from (parent bit)\r
 470           tr.SetTag(GetMCAnalysisUtils()->CheckOrigin(tr.GetLabel(),GetReader(),tr.GetInputFileIndex())); //Gets a tag bit which contains info about super grandfather particle. Use (tag&(1<<11)), 11 for direct b, and 9 for B->C\r
 471           \r
 472           if(tr.GetTag()&(1<<9)||tr.GetTag()&(1<<11)) //MC particle from b-decay\r
 473             fhNVTXMC->Fill(dvmbtag);\r
 474           \r
 475           if(fWriteNtuple){\r
 476             fNElec++;\r
 477             fNElecEv++;\r
 478             electrons->Fill(fNElec,fEventNumber,fNElecEv,pairs1,start,stop,tpcflag,trdflag,emcflag,pdg,tr.GetTag(),tr.GetBtag(),tr.Pt());\r
 479             \r
 480             \r
 481           }\r
 482           \r
 483           if(GetDebug() > 0) \r
 484             printf("AliAnaBtag::MakeAnalysisFillAOD() - Origin of candidate (parent bit) %d\n",tr.GetTag());\r
 485         }//MonteCarlo MC done\r
 486         \r
 487         AddAODParticle(tr);             \r
 488       }//electron\r
 489       \r
 490     } //pid check\r
 491   }//track loop                         \r
 492   if(GetDebug() > 1) printf("AliAnaBtag::MakeAnalysisFillAOD()  End fill AODs \n");  \r
 493   \r
 494 }\r
 495 \r
 496 //__________________________________________________________________\r
 497 void  AliAnaBtag::MakeAnalysisFillHistograms() \r
 498 {\r
 499   //Do analysis and fill histograms\r
 500   AliStack * stack = 0x0;\r
 501   //   TParticle * primary = 0x0;\r
 502   if(IsDataMC()) {\r
 503     if(GetReader()->ReadStack()){\r
 504       stack =  GetMCStack() ;      \r
 505       if(!stack)\r
 506         printf("AliAnaBtag::MakeAnalysisFillHistograms() *** no stack ***: \n");   \r
 507     }\r
 508   }// is data and MC\r
 509   \r
 510 \r
 511   double maxjetEta=-4.;\r
 512   double maxjetPhi=-4.;\r
 513   \r
 514   Int_t njets = 0; \r
 515   if(GetReader()->GetOutputEvent()) njets = (GetReader()->GetOutputEvent())->GetNJets();\r
 516   if(njets > 0) {\r
 517     if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - Jet AOD branch has %d jets.  Performing b-jet tag analysis\n",njets);\r
 518     \r
 519     ///////////////////////////////////Jet loop//////////////////////////////////////////////\r
 520     for(Int_t ijet = 0; ijet < njets ; ijet++) {\r
 521       AliAODJet * jet = (AliAODJet*)(GetReader()->GetOutputEvent())->GetJet(ijet) ;\r
 522       fhJets->Fill(jet->Pt(),1);                                                                                ////////////////FILL\r
 523       fhJetsAllEtaPhi->Fill(jet->Eta(),jet->Phi());\r
 524       if(jet->Pt() < 0.) continue; //This has to be adjusted depending on pp or AA!\r
 525       \r
 526       //Geometric EMCAL cut\r
 527       if(TMath::Abs(jet->Eta()) > 0.3) continue;\r
 528       if(jet->Phi() < 1.8 || jet->Phi() > 2.9) continue; //This is BAD FIXME \r
 529       fhJets->Fill(jet->Pt(),4); //All jets after geometric cut                                                 ////////////////FILL\r
 530       \r
 531       Bool_t leadJet  = kFALSE;\r
 532       if (ijet==0) {leadJet= kTRUE; fhJets->Fill(jet->Pt(),5);} //Leading jets                                   ////////////////FILL\r
 533       \r
 534       /////////////////////////Track loop in Jet////////////////////////////////////\r
 535       Bool_t dvmJet = kFALSE; \r
 536       Bool_t dvmMCJet = kFALSE; \r
 537       TRefArray* rt = jet->GetRefTracks();\r
 538       Int_t ntrk = rt->GetEntries();\r
 539       \r
 540       for(Int_t itrk = 0; itrk < ntrk; itrk++) {\r
 541         AliAODTrack* jetTrack = (AliAODTrack*)jet->GetTrack(itrk);\r
 542         \r
 543         Int_t trackId = jetTrack->GetID(); //get the index in the reader\r
 544         Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
 545         if(GetDebug() > 3) printf("AliAnaBtag::CheckIfBjet() - aod branch entries %d\n", naod);\r
 546         \r
 547         for(Int_t iaod = 0; iaod < naod ; iaod++){\r
 548           AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
 549           Int_t electronlabel = ele->GetTrackLabel(0);\r
 550           if(electronlabel != trackId) continue;  //skip to the next one if they don't match\r
 551           if(ele->GetBtag()>0){\r
 552             dvmJet = kTRUE;\r
 553             if(ele->GetTag()&(1<<9) || ele->GetTag()&(1<<11) )\r
 554               dvmMCJet = kTRUE;\r
 555           }\r
 556         } //Electron check\r
 557         \r
 558       }//Track loop of jet tracks\r
 559       \r
 560       if(dvmJet) fhJets->Fill(jet->Pt(),6);                                                                      ////////////////FILL\r
 561       //MC stuff\r
 562       if(IsDataMC()) {\r
 563         Bool_t isTrueBjet = IsMcBJet(jet->Eta(), jet->Phi());\r
 564         //Bool_t isTrueDjet = IsMcDJet(jet->Eta(), jet->Phi());\r
 565         if(dvmJet){\r
 566           if(dvmMCJet){\r
 567             fhJets->Fill(jet->Pt(),8);   //True                                                                   ////////////////FILL\r
 568           }else{\r
 569             fhJets->Fill(jet->Pt(),9);   //False                                                                 ////////////////FILL\r
 570           }\r
 571           if(isTrueBjet){ \r
 572             fhJets->Fill(jet->Pt(),10);     //True                                                                 ////////////////FILL\r
 573           }else{\r
 574             fhJets->Fill(jet->Pt(),11);     //False                                                                ////////////////FILL\r
 575           }\r
 576         }\r
 577         if(isTrueBjet){ \r
 578           fhJets->Fill(jet->Pt(),12);     //True                                                                 ////////////////FILL\r
 579         }else{\r
 580           fhJets->Fill(jet->Pt(),13);     //False                                                                ////////////////FILL\r
 581         }\r
 582 \r
 583       }//MC stuff\r
 584 \r
 585     }//jet loop\r
 586   } //jets exist\r
 587   \r
 588   //Electron loop, read back electrons, fill histos; mostly photonic shit.\r
 589   Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
 590   if(GetDebug() > 0) printf("AliAnaBtag::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);\r
 591   \r
 592   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
 593     AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
 594     Int_t pdg = ele->GetPdg();\r
 595     \r
 596     \r
 597     if(TMath::Abs(pdg) != AliCaloPID::kElectron) continue; //not necessary..\r
 598     \r
 599     //MC tag of this electron\r
 600     //    Int_t mctag = ele->GetTag();\r
 601     \r
 602     \r
 603     fhElectrons->Fill(ele->Pt(),1); //All electrons\r
 604     Bool_t photonic = kFALSE;\r
 605     photonic = PhotonicV0(ele->GetTrackLabel(0)); //check against V0s\r
 606     if(!photonic) fhElectrons->Fill(ele->Pt(),3); //nonphotonic electrons\r
 607     if(photonic) fhElectrons->Fill(ele->Pt(),4);  //photonic electrons\r
 608     \r
 609     //Fill electron histograms \r
 610     Float_t phiele = ele->Phi();\r
 611     Float_t etaele = ele->Eta();\r
 612     \r
 613     \r
 614     if(ele->GetBtag()>0){ // removed bit tag shit\r
 615       fhElectrons->Fill(ele->Pt(),5);\r
 616       if(!photonic) fhElectrons->Fill(ele->Pt(),6);\r
 617       if(photonic) fhElectrons->Fill(ele->Pt(),7);\r
 618       fhJetsLeadingBElectronEtaPhi->Fill(maxjetEta,maxjetPhi); \r
 619       double deta=etaele-maxjetEta;\r
 620       double dphi=phiele-maxjetPhi;\r
 621       \r
 622       //double r = sqrt((dphi*dphi)+(deta*deta));\r
 623       fhBJetElectronDetaDphi->Fill(deta,dphi);\r
 624       \r
 625     }\r
 626     \r
 627   }//electron aod loop\r
 628   \r
 629 }\r
 630 \r
 631 //__________________________________________________________________\r
 632 Int_t AliAnaBtag::GetDVMBtag(AliAODTrack * tr, Int_t &pair, Int_t &start, Int_t &stop)\r
 633 {\r
 634   //This method uses the Displaced Vertex between electron-hadron\r
 635   //pairs and the primary vertex to determine whether an electron is\r
 636   //likely from a B hadron.\r
 637   Int_t pairstart = fNPair;\r
 638   Int_t pairstop = 0;\r
 639   Int_t pairn = 0;\r
 640   Int_t ncls1 = 0;\r
 641   for(Int_t l = 0; l < 6; l++) if(TESTBIT(tr->GetITSClusterMap(),l)) ncls1++;\r
 642 \r
 643   if (ncls1 < fITSCut) return 0;\r
 644 \r
 645   Double_t imp[2] = {-999.,-999.}; Double_t cov[3] = {-999.,-999.,-999.};\r
 646   Bool_t dcaOkay = GetDCA(tr,imp,cov);  //homegrown dca calculation until AOD is fixed                  \r
 647   if(!dcaOkay) {\r
 648     printf("AliAnaBtag::GetDVMBtag - Problem computing DCA to primary vertex for track %d",tr->GetID());\r
 649     return 0;\r
 650   }\r
 651 \r
 652   if (TMath::Abs(imp[0])   > fImpactCut ) return 0;\r
 653   if (TMath::Abs(imp[1])   > fImpactCut ) return 0;\r
 654 \r
 655   Int_t nvtx = 0;\r
 656 //   Int_t nvtx2 = 0;\r
 657 //   Int_t nvtx3 = 0;\r
 658 \r
 659   for (Int_t k2 =0; k2 < GetCTSTracks()->GetEntriesFast() ; k2++) {\r
 660 \r
 661     //loop over assoc\r
 662     AliAODTrack* track2 = (AliAODTrack*)GetCTSTracks()->At(k2);\r
 663     Int_t id1 = tr->GetID();\r
 664     Int_t id2 = track2->GetID();\r
 665     if(id1 == id2) continue;\r
 666 \r
 667     Int_t ncls2 = 0;\r
 668     for(Int_t l = 0; l < 6; l++) if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;\r
 669     if (ncls2 < fITSCut) continue;\r
 670 \r
 671 \r
 672     Double_t sdca=0,pdca=0,minv=0,pth=0,massphoton=0,decaylength=0;\r
 673 \r
 674     sdca = ComputeSignDca(tr,track2,minv,pdca,massphoton,decaylength);\r
 675     pth=track2->Pt();\r
 676 \r
 677 \r
 678 \r
 679 \r
 680     Double_t dphi = tr->Phi() - track2->Phi();\r
 681     if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();\r
 682     if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();\r
 683     Double_t deta = tr->Eta() - track2->Eta();\r
 684     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
 685 \r
 686     if(dr > fDrCut) continue;\r
 687     fNPair++;\r
 688     pairn++;\r
 689     if(GetDebug() > 0) \r
 690       printf("pairs: %d !!!!!!!!!!! \n",fNPair);\r
 691     if(fWriteNtuple){\r
 692       pairs->Fill(fNPair,fNElec,fEventNumber,pdca,sdca,minv,pth,massphoton,decaylength);\r
 693     }\r
 694     pairstop=fNPair;\r
 695     //pairs->Fill(1,1,1,1,1,1,1);\r
 696     fhPairPt->Fill(pth);\r
 697     if(decaylength>1.0) continue;\r
 698     if(tr->Pt()<6. && pth>0.4 && minv>1.4 && pdca<0.025 && sdca>0.06 && massphoton>0.1) nvtx++; \r
 699     if(tr->Pt()>6.&&tr->Pt()<10. && pth>0.2 && minv>1.7 && pdca<0.012 && sdca>0.06 && massphoton>0.1) nvtx++;\r
 700     if(tr->Pt()>10.&& pth>0.6 && minv>1.5 && pdca<0.14 && sdca>0.04 && massphoton>0.1) nvtx++;\r
 701 \r
 702 \r
 703 \r
 704   } //loop over hadrons\r
 705 \r
 706   if(GetDebug() > 0) printf("AliAnaBtag::GetDVMBtag - result of btagging: %d \n",nvtx);\r
 707 \r
 708 \r
 709   pair=pairn;\r
 710   start=pairstart+1;\r
 711   stop=pairstop;\r
 712   if(GetDebug() > 0) printf("End of DVM, pair counts: pairs: %d, start %d, stop %d. \n",pair,start,stop);\r
 713   return nvtx;\r
 714 }\r
 715 \r
 716 //__________________________________________________________________\r
 717 Double_t AliAnaBtag::ComputeSignDca(AliAODTrack *tr, AliAODTrack *tr2 , Double_t &masscut, Double_t &pdcacut, Double_t &massphoton, Double_t &decay)\r
 718 {\r
 719   //Compute the signed dca between two tracks\r
 720   //and return the result\r
 721 \r
 722   Double_t signDca=-999.;\r
 723   if(GetDebug() > 2 ) printf(">>ComputeSdca:: track1 %d, track2 %d, masscut %f \n", tr->GetLabel(), tr2->GetLabel(), masscut);\r
 724 \r
 725   //=====Now calculate DCA between both tracks=======  \r
 726   Double_t massE = 0.000511;\r
 727   Double_t massK = 0.493677;\r
 728   Double_t vertex[3] = {-999.,-999.,-999}; //vertex\r
 729 \r
 730   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
 731     GetVertex(vertex); //If only one file, get the vertex from there  \r
 732   }\r
 733   \r
 734   TVector3 primV(vertex[0],vertex[1],vertex[2]) ; \r
 735   if(GetDebug()>0) printf(">>ComputeSdca:: primary vertex = %2.2f,%2.2f,%2.2f \n",vertex[0],vertex[1],vertex[2]) ;\r
 736 \r
 737   AliExternalTrackParam *param1 = new AliExternalTrackParam(tr);\r
 738   AliExternalTrackParam *param2 = new AliExternalTrackParam(tr2);\r
 739 \r
 740   Double_t bfield[3];\r
 741   param1->GetBxByBz(bfield);\r
 742   bfield[0]=0;\r
 743   bfield[1]=0;\r
 744   bfield[2]=5.;\r
 745   Double_t bz = 5.; // = param1->GetBz();\r
 746   Double_t xplane1 = 0.; Double_t xplane2 = 0.;\r
 747   Double_t pairdca = param1->GetDCA(param2,bz,xplane1,xplane2);\r
 748 \r
 749   param1->PropagateToBxByBz(xplane1,bfield);\r
 750   param2->PropagateToBxByBz(xplane2,bfield);\r
 751 \r
 752   Int_t id1 = 0, id2 = 0;\r
 753   AliESDv0 bvertex(*param1,id1,*param2,id2);\r
 754   Double_t vx,vy,vz;\r
 755   bvertex.GetXYZ(vx,vy,vz);\r
 756 \r
 757   Double_t emom[3];\r
 758   Double_t hmom[3];\r
 759   param1->PxPyPz(emom);\r
 760   param2->PxPyPz(hmom);\r
 761   TVector3 emomAtB(emom[0],emom[1],emom[2]);\r
 762   TVector3 hmomAtB(hmom[0],hmom[1],hmom[2]);\r
 763   TVector3 secvtxpt(vx,vy,vz);\r
 764   TVector3 decayvector(0,0,0);\r
 765   decayvector = secvtxpt - primV; //decay vector from PrimVtx\r
 766   Double_t decaylength = decayvector.Mag();\r
 767   decay=decaylength;\r
 768 \r
 769   if(GetDebug() > 0) {\r
 770     printf(">>ComputeSdca:: mom1=%f, mom2=%f \n", emomAtB.Perp(), hmomAtB.Perp() );\r
 771     printf(">>ComputeSdca:: pairDCA=%f, length=%f \n", pairdca,decaylength );\r
 772   }\r
 773 \r
 774   if (emomAtB.Mag()>0 /*&& decaylength < fDecayLenCut*/ ) {\r
 775     TVector3 sumMom = emomAtB+hmomAtB;\r
 776     Double_t ener1 = sqrt(pow(emomAtB.Mag(),2) + massE*massE);\r
 777     Double_t ener2 = sqrt(pow(hmomAtB.Mag(),2) + massK*massK);\r
 778     Double_t ener3 = sqrt(pow(hmomAtB.Mag(),2) + massE*massE);\r
 779     Double_t mass = sqrt(pow((ener1+ener2),2) - pow(sumMom.Mag(),2));\r
 780     Double_t massPhot = sqrt(pow((ener1+ener3),2) - pow(sumMom.Mag(),2));\r
 781     Double_t sDca = decayvector.Dot(emomAtB)/emomAtB.Mag();\r
 782     pdcacut=pairdca;//\r
 783     masscut=mass; // \r
 784     massphoton=massPhot;//                                                                     Send it back!\r
 785     signDca = sDca;\r
 786     if(GetDebug() > 0) printf("\t>>ComputeSdca:: mass=%f \n", mass);\r
 787     if(GetDebug() > 0) printf("\t>>ComputeSdca:: sec vtx-signdca :%f\n",signDca);\r
 788   }\r
 789   //clean up\r
 790   delete param1;\r
 791   delete param2;\r
 792   return signDca;\r
 793 }\r
 794 \r
 795 \r
 796 //__________________________________________________________________\r
 797 Bool_t AliAnaBtag::PhotonicV0(Int_t id) \r
 798 {\r
 799   //This method checks to see whether a track that has been flagged as\r
 800   //an electron was determined to match to a V0 candidate with\r
 801   //invariant mass consistent with photon conversion\r
 802   Bool_t itIS = kFALSE;\r
 803   Double_t massEta = 0.547;\r
 804   Double_t massRho0 = 0.770;\r
 805   Double_t massOmega = 0.782;\r
 806   Double_t massPhi = 1.020;\r
 807   \r
 808   //---Get V0s---\r
 809   AliAODEvent *aod = (AliAODEvent*) GetReader()->GetInputEvent();\r
 810   int nv0s = aod->GetNumberOfV0s();\r
 811   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
 812     AliAODv0 *v0 = aod->GetV0(iV0);\r
 813     if (!v0) continue;\r
 814     double radius = v0->RadiusV0();\r
 815     double mass = v0->InvMass2Prongs(0,1,11,11);\r
 816     if(GetDebug() > 0) {\r
 817       printf("## PhotonicV0() :: v0: %d, radius: %f \n", iV0 , radius );\r
 818       printf("## PhotonicV0() :: neg-id: %d, pos-id: %d, THIS id: %d\n", v0->GetNegID(), v0->GetPosID(), id);\r
 819       printf("## PhotonicV0() :: Minv(e,e): %f \n", v0->InvMass2Prongs(0,1,11,11) );\r
 820     }\r
 821     if (mass < 0.100 ||\r
 822         (mass > massEta-0.05 || mass < massEta+0.05) ||\r
 823         (mass > massRho0-0.05 || mass < massRho0+0.05) ||\r
 824         (mass > massOmega-0.05 || mass < massOmega+0.05) ||\r
 825         (mass > massPhi-0.05 || mass < massPhi+0.05)) {\r
 826       if ( id == v0->GetNegID() || id == v0->GetPosID()) {\r
 827         itIS=kTRUE;\r
 828         if(GetDebug() > 0) printf("## PhotonicV0() :: It's a conversion electron!!! \n" );\r
 829       }\r
 830     } }\r
 831   return itIS;\r
 832 }\r
 833 \r
 834 //__________________________________________________________________\r
 835 Bool_t AliAnaBtag::GetDCA(const AliAODTrack* track,Double_t impPar[2], Double_t cov[3]) \r
 836 {\r
 837   //Use the Event vertex and AOD track information to get\r
 838   //a real impact parameter for the track\r
 839   Double_t maxD = 100000.; //max transverse IP\r
 840   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) {\r
 841     AliVEvent* ve = (AliVEvent*)GetReader()->GetInputEvent();\r
 842     AliVVertex *vv = (AliVVertex*)ve->GetPrimaryVertex();\r
 843     AliESDtrack esdTrack(track);\r
 844     Double_t bfield[3];\r
 845     esdTrack.GetBxByBz(bfield);\r
 846     bfield[0]=0;\r
 847     bfield[1]=0;\r
 848     bfield[2]=5.;\r
 849     \r
 850     Bool_t gotit = esdTrack.PropagateToDCABxByBz(vv,bfield,maxD,impPar,cov);\r
 851     return gotit;\r
 852   }\r
 853   return kFALSE;\r
 854 }\r
 855 //__________________________________________________________________\r
 856 Bool_t AliAnaBtag::CheckIfBjet(const AliAODTrack* track)\r
 857 {\r
 858   Bool_t bjet = kFALSE;\r
 859   Int_t trackId = track->GetID(); //get the index in the reader\r
 860   Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
 861   if(GetDebug() > 3) printf("AliAnaBtag::CheckIfBjet() - aod branch entries %d\n", naod);\r
 862 \r
 863   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
 864     AliAODPWG4Particle* ele =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
 865     Int_t electronlabel = ele->GetTrackLabel(0);\r
 866     if(electronlabel != trackId) continue;  //skip to the next one if they don't match\r
 867     if(ele->GetBtag()>0)\r
 868       bjet = kTRUE;\r
 869   } \r
 870 \r
 871   return bjet;\r
 872 } \r
 873 \r
 874 //__________________________________________________________________\r
 875 AliAODMCParticle* AliAnaBtag::GetMCParticle(Int_t ipart) \r
 876 {\r
 877   //Get the MC particle at position ipart\r
 878   \r
 879   AliAODMCParticle* aodprimary = 0x0;\r
 880   TClonesArray * mcparticles0 = 0x0;\r
 881   \r
 882   if(GetReader()->ReadAODMCParticles()){\r
 883     //Get the list of MC particles                                                                                                                           \r
 884     mcparticles0 = GetReader()->GetAODMCParticles(0);\r
 885     if(!mcparticles0) {\r
 886       if(GetDebug() > 0)printf("AliAnaBtag::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");\r
 887     }\r
 888     else{\r
 889       Int_t npart0 = mcparticles0->GetEntriesFast();\r
 890       if(ipart < npart0) aodprimary = (AliAODMCParticle*)mcparticles0->At(ipart);\r
 891       \r
 892       if(!aodprimary) {\r
 893         printf("AliAnaBtag::GetMCParticle() *** no primary ***:  label %d \n", ipart);\r
 894         return 0x0;\r
 895       }\r
 896     }\r
 897   } else {\r
 898     printf("AliAnaBtag::GetMCParticle() - Asked for AliAODMCParticle but we have a stack reader.\n");\r
 899   }\r
 900   return aodprimary;\r
 901 }\r
 902 \r
 903 //__________________________________________________________________\r
 904 Bool_t  AliAnaBtag::IsMcBJet(Double_t jeta, Double_t jphi)\r
 905 {\r
 906   //Check the jet eta,phi against that of the b-quark\r
 907   //to decide whether it is an MC B-jet\r
 908   Bool_t bjet=kFALSE;\r
 909   AliStack* stack = 0x0;\r
 910   \r
 911   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
 912 \r
 913     Double_t pphi = -999.;\r
 914     Double_t peta = -999.;\r
 915     Int_t pdg = 0;\r
 916     if(GetReader()->ReadStack()) {\r
 917       stack = GetMCStack();\r
 918       if(!stack) {\r
 919         printf("AliAnaBtag::IsMCBJet() *** no stack ***: \n");\r
 920         return kFALSE;\r
 921       }\r
 922       TParticle* primary = stack->Particle(ipart);\r
 923       if (!primary) continue;\r
 924       pdg = primary->GetPdgCode();\r
 925       pphi = primary->Phi();\r
 926       peta = primary->Eta();\r
 927     } else if(GetReader()->ReadAODMCParticles()) {\r
 928       AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
 929       if(!aodprimary) continue;\r
 930       pdg = aodprimary->GetPdgCode();\r
 931       pphi = aodprimary->Phi();\r
 932       peta = aodprimary->Eta();\r
 933     }\r
 934     if ( TMath::Abs(pdg) != 5) continue;\r
 935     Double_t dphi = jphi - pphi;\r
 936     Double_t deta = jeta - peta;\r
 937     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
 938     \r
 939     if (dr < 0.2) {\r
 940       bjet=kTRUE;\r
 941       break;\r
 942     }\r
 943   }\r
 944   return bjet;\r
 945 }\r
 946 \r
 947 //__________________________________________________________________\r
 948 Bool_t  AliAnaBtag::IsMcDJet(Double_t jeta, Double_t jphi)\r
 949 {\r
 950   //Check if this jet is a charm jet\r
 951   Bool_t cjet=kFALSE;\r
 952 \r
 953   AliStack* stack = 0x0;\r
 954 \r
 955   for(Int_t ipart = 0; ipart < 100; ipart++) {\r
 956     \r
 957     Double_t pphi = -999.;\r
 958     Double_t peta = -999.;\r
 959     Int_t pdg = 0;\r
 960     if(GetReader()->ReadStack()) {\r
 961       stack = GetMCStack();\r
 962       if(!stack) {\r
 963         printf("AliAnaBtag::IsMCDJet() *** no stack ***: \n");\r
 964         return kFALSE;\r
 965       }\r
 966       TParticle* primary = stack->Particle(ipart);\r
 967       if (!primary) continue;\r
 968       pdg = primary->GetPdgCode();\r
 969       pphi = primary->Phi();\r
 970       peta = primary->Eta();\r
 971     } else if(GetReader()->ReadAODMCParticles()) {\r
 972       AliAODMCParticle* aodprimary = GetMCParticle(ipart);\r
 973       if(!aodprimary) continue;\r
 974       pdg = aodprimary->GetPdgCode();\r
 975       pphi = aodprimary->Phi();\r
 976       peta = aodprimary->Eta();\r
 977     }\r
 978 \r
 979     if ( TMath::Abs(pdg) != 4) continue;\r
 980     Double_t dphi = jphi - pphi;\r
 981     Double_t deta = jeta - peta;\r
 982     Double_t dr = sqrt(deta*deta + dphi*dphi);\r
 983    \r
 984     if (dr < 0.2) {\r
 985       cjet=kTRUE;\r
 986       break;\r
 987     }\r
 988   }\r
 989   return cjet;\r
 990 }\r
 991 \r
 992 \r
 993 //__________________________________________________________________\r
 994 void  AliAnaBtag::Terminate(TList* outputList)\r
 995 {\r
 996   printf(" AliAnaBtag::Terminate()  *** %s Report: %d outputs/histograms \n", GetName(), outputList->GetEntries()) ;\r
 997 }\r