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