[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskEmcalJetFlavourTagExample.cxx

/**************************************************************************
* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

/*
*   This sample task demostrates the basics of tagging a jet. The PID response
*   is retrived for both the TPC.  A Hadronic tag is implemented for
*   clusters matched to tracks in a jet with a E/P < 0.2
*
*   Author: Andrew Castro (UTK) and Joel Mazer (UTK)
*/

#include "AliAnalysisTaskEmcalJetFlavourTagExample.h"

// general ROOT includes                                                                                                                                                  
#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THnSparse.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TParameter.h>
#include <TParticle.h>
#include <TTree.h>
#include <TVector3.h>
#include <TObjArray.h>

// AliROOT includes                                                                                                                         
#include "AliAODEvent.h"
#include "AliESDEvent.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTask.h"
#include "AliCentrality.h"
#include "AliEmcalJet.h"
#include "AliAODJet.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include <AliVEvent.h>
#include <AliVParticle.h>
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"
#include "AliEmcalParticle.h"
#include "AliESDCaloCluster.h"
#include <AliESDtrackCuts.h>
#include "AliPID.h"

// event handler (and pico's) includes                                                                                                      
#include <AliInputEventHandler.h>
#include <AliVEventHandler.h>
#include "AliESDInputHandler.h"
#include "AliPicoTrack.h"
#include "AliEventPoolManager.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"

// PID includes                                                                                                                             
#include "AliPIDResponse.h"
#include "AliTPCPIDResponse.h"
#include "AliESDpid.h"

// magnetic field includes
#include "TGeoGlobalMagField.h"
#include "AliMagF.h"

using std::cout;
using std::endl;

ClassImp(AliAnalysisTaskEmcalJetFlavourTagExample)

//________________________________________________________________________
AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample() : 
  AliAnalysisTaskEmcalJet("heavyF",kFALSE), 
  event(0),
  fCuts(0),
  fPhimin(-10), fPhimax(10),
  fEtamin(-0.9), fEtamax(0.9),
  fAreacut(0.0),
  fJetHIpt(50.0),
  fTrackPtCut(2.0),
  fTrackEta(0.9),
  fesdTrackCuts(0),
  fPIDResponse(0x0), fTPCResponse(),
  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
  fESD(0), fAOD(0),
  fHistJetPhi(0),
  fHistCorJetPt(0), fHistJetPt(0),
  fHistHighJetPt(0),
  fHistnSigElecPt(0),
  fHistdEdXvsPt(0),
  fHistnJetTrackvnJetClusters(0)
{
  // Default constructor.


  SetMakeGeneralHistograms(kTRUE);

}

//________________________________________________________________________
AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample(const char *name) :
  AliAnalysisTaskEmcalJet(name,kTRUE),
  event(0),
  fCuts(0),
  fPhimin(-10), fPhimax(10),
  fEtamin(-0.9), fEtamax(0.9),
  fAreacut(0.0),
  fJetHIpt(50.0),
  fTrackPtCut(2.0),
  fTrackEta(0.9),
  fesdTrackCuts(0),
  fPIDResponse(0x0), fTPCResponse(),
  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
  fESD(0), fAOD(0),
  fHistJetPhi(0),
  fHistCorJetPt(0), fHistJetPt(0),
  fHistHighJetPt(0),
  fHistnSigElecPt(0),
  fHistdEdXvsPt(0),
  fHistnJetTrackvnJetClusters(0)
{ 

   SetMakeGeneralHistograms(kTRUE);
 
   DefineInput(0,TChain::Class());
   DefineOutput(1, TList::Class());
}

//_______________________________________________________________________
AliAnalysisTaskEmcalJetFlavourTagExample::~AliAnalysisTaskEmcalJetFlavourTagExample()
{
  // destructor
  //
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetFlavourTagExample::UserCreateOutputObjects()
{
  if (! fCreateHisto)
    return;
  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  //fJetsCont           = GetJetContainer(0);
  if(fJetsCont) { //get particles and clusters connected to jets
    fTracksJetCont       = fJetsCont->GetParticleContainer();
    fCaloClustersJetCont = fJetsCont->GetClusterContainer();
  }
 else {        //no jets, just analysis tracks and clusters
  fTracksCont       = GetParticleContainer(0);
  fCaloClustersCont = GetClusterContainer(0);
}
fTracksCont->SetClassName("AliVTrack");
fCaloClustersCont->SetClassName("AliVCluster");

  fHistJetPhi                = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
  fHistJetPt                 = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
  fOutput->Add(fHistJetPhi);
  fOutput->Add(fHistJetPt);

 
  TString histname;

  fHistCorJetPt		            = new TH1F("CorrJetPt", "CorrJetPt", 300, -100, 200);
  fHistnSigElecPt             = new TH2F("nsigvsPt(TPC)","nSigmaElectronTPC_v_Pt",60,0,60,40,-10,10);
  fHistdEdXvsPt               = new TH2F("dEdXvsTrackPt","dEdXvsTrackPt",60,0,60,80,0,80);
  fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
  fHistHighJetPt              = new TH1F("HighestPtJetPerEvent","HighJetPt",80,0,80);
    
  TString name;
  TString title;

  fOutput->Add(fHistCorJetPt);
  fOutput->Add(fHistnSigElecPt);
  fOutput->Add(fHistdEdXvsPt);
  fOutput->Add(fHistnJetTrackvnJetClusters);
  fOutput->Add(fHistHighJetPt);


  // ****************************** PID *****************************************************                                               
  // set up PID handler                                                                                                                     
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
  if(!inputHandler) {
    AliFatal("Input handler needed");
  }

  // PID response object                                                                                                                    
  //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();                                                                         
  //  inputHandler->CreatePIDResponse(fIsMC);         // needed to create object, why though?                                                 
  fPIDResponse = inputHandler->GetPIDResponse();
  if (!fPIDResponse) {
    AliError("PIDResponse object was not created");
  }
  // ***************************************************************************************

  PostData(1, fOutput);

}

//________________________________________________________
void AliAnalysisTaskEmcalJetFlavourTagExample::ExecOnce()
{
  //  Initialize the analysis
  AliAnalysisTaskEmcalJet::ExecOnce();
  
  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;


} // end of ExecOnce

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetFlavourTagExample::Run()
{
  // check to see if we have any tracks
  if (!fTracks)  return kTRUE;
  if (!fJets)  return kTRUE;

  // what kind of event do we have: AOD or ESD?
  Bool_t useAOD;
  if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
  else useAOD = kFALSE;
  
  // if we have ESD event, set up ESD object
  if(!useAOD){
    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    if (!fESD) {
      AliError(Form("ERROR: fESD not available\n"));
      return kTRUE;
    }
  }

  // if we have AOD event, set up AOD object
  if(useAOD){
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
    if(!fAOD) {
      AliError(Form("ERROR: fAOD not available\n"));
      return kTRUE;
    }
  }
  
  // get magnetic field info for DCA
  Double_t  MagF = fESD->GetMagneticField();
  Double_t MagSign = 1.0;
  if(MagF<0)MagSign = -1.0;
  // set magnetic field
  if (!TGeoGlobalMagField::Instance()->GetField()) {
    AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
    TGeoGlobalMagField::Instance()->SetField(field);
  }

  // get vertex information
  Double_t fvertex[3]={0,0,0};
  InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
  //Double_t zVtx=fvertex[2];

  // create pointer to list of input event                                                                                                  
  TList *list = InputEvent()->GetList();
  if(!list) {
    AliError(Form("ERROR: list not attached\n"));
    return kTRUE;
  }

  // background density                                                                                                                                                                                                                               
  fRhoVal = fRho->GetVal();

  // initialize TClonesArray pointers to jets and tracks                                                                                    
  TClonesArray *jets = 0;
  //TClonesArray *tracks = 0;
  //TClonesArray *clusters = 0;
  
  // get Jets object                                                                                                                        
  jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
  if(!jets){
    AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
    return kTRUE;
  } // verify existence of jets
  
  // get number of jets and tracks                                                                                                          
  const Int_t Njets = jets->GetEntries();
  if(Njets<1)     return kTRUE;
  
  if (fTracksCont) {
    AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
    while(track) {
      track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
    }
  }
  if (fCaloClustersCont) {
    AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
    while(cluster) {
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, fVertex);
      cluster = fCaloClustersCont->GetNextAcceptCluster();
    }
  }
  
    //  Start Jet Analysis
    // initialize jet parameters
    Int_t ijethi=-1;
    Double_t highestjetpt=0.0;
  
  // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
  for (Int_t ijet = 0; ijet < Njets; ijet++){
    // get our jets
    AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
    if (!jet) continue;
    
    // apply jet cuts
    if(!AcceptMyJet(jet)) continue;
    
    if(highestjetpt<jet->Pt()){
      ijethi=ijet;
      highestjetpt=jet->Pt();
    }
  } // end of looping over jets
  
    fHistHighJetPt->Fill(ijethi);
  
 // **********************************************************************
 //                JET LOOP
 // **********************************************************************
  
    // loop over jets in the event and make appropriate cuts
    for (Int_t iJets = 0; iJets < Njets; ++iJets) {
      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
      if (!jet)  // see if we have a jet
        continue;
      
      // apply jet cuts
      if(!AcceptMyJet(jet)) continue;
      
      //AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
      //jet->AddFlavourTag(tag);
    
 
      Int_t JetClusters = jet->GetNumberOfClusters();
      Int_t JetTracks = jet -> GetNumberOfTracks();
      fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
      // Initializations and Calculations
      Double_t jetPt = -500;                                     // initialize corr jet pt LOCAL
      Double_t jetarea = -500;                                   // initialize jet area
      jetPt = jet->Pt() - jetarea*fRhoVal;                       // semi-corrected pT of jet from GLOBAL rho value
      fHistCorJetPt->Fill(jetPt);
      
      Bool_t bkgrnd1 = kFALSE;
      
      if(jet->Pt() > fJetHIpt) {
        if(!fTracksCont || !fCaloClustersCont) continue;
          
        //******************************Cluster Matched To Closest Track
        //**************************************************************
        Int_t NumbTrackContainer = -999;
        NumbTrackContainer = fTracksCont->GetNParticles();
        for(int iTracks = 0; iTracks <= NumbTrackContainer; iTracks++){
          AliVTrack *AcceptedTrack =static_cast<AliVTrack*>(fTracksCont->GetParticle(iTracks));
          if(!AcceptedTrack){
            AliError(Form("Couldn't get AliVTrack Container %d\n", iTracks));
            continue;
          }
          if(!IsJetTrack(jet,iTracks,kFALSE))continue;
          //Get matched cluster
          Int_t emc1 = AcceptedTrack->GetEMCALcluster();
        
          Double_t acceptTrackP = AcceptedTrack->P();
          Double_t acceptTrackPt = AcceptedTrack->Pt();
          Double_t nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(AcceptedTrack,AliPID::kElectron);
          fHistnSigElecPt->Fill(acceptTrackPt,nSigmaElectron_TPC);

          AliESDtrack *ESDacceptedTrack = static_cast<AliESDtrack*>(AcceptedTrack);
          if(!ESDacceptedTrack){
            AliError(Form("Couldn't get AliESDTrack %d\n", iTracks));
            continue;
          }
          Double_t dEdx = AcceptedTrack->GetTPCsignal();
                              fHistdEdXvsPt->Fill(acceptTrackPt,dEdx);
          if(fCaloClustersCont && emc1>=0) {
            AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
            if(!clusMatch){
              AliError(Form("Couldn't get matched AliVCluster %d\n", emc1));
              continue;
            }
 
            Double_t mClusterE = clusMatch->E();
            Double_t EovP_mc = -999;
            EovP_mc = mClusterE/acceptTrackP;
          
            if(EovP_mc < 0.2) bkgrnd1 = kTRUE;
          }
          
      } //loop over tracks for matching to closest cluster

     
    } // highest pt jet cut
    
      
      if(bkgrnd1 == kTRUE) {
        AliEmcalJet::EFlavourTag tag=AliEmcalJet::kBckgrd1;
        jet->AddFlavourTag(tag);
    }
        
  } // LOOP over JETS in event


  return kTRUE;
  
}
//________________________________________________________________________
void AliAnalysisTaskEmcalJetFlavourTagExample::Terminate(Option_t *)
{
  cout<<"###########################"<<endl;
  cout<<"####   Task Finished   ####"<<endl;
  cout<<"###########################"<<endl;
  cout<<"###########################"<<endl;
} // end of terminate

//________________________________________________________________________
Int_t AliAnalysisTaskEmcalJetFlavourTagExample::AcceptMyJet(AliEmcalJet *jet) {
  //applies all jet cuts except pt
  if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
  if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
  if (jet->Area()<fAreacut) return 0;
  // prevents 0 area jets from sneaking by when area cut == 0
  if (jet->Area()==0) return 0;
  //exclude jets with extremely high pt tracks which are likely misreconstructed
  if(jet->MaxTrackPt()>100) return 0;
  
  //passed all above cuts
  return 1;
}
Commit	Line	Data
9382bd6e	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	* This sample task demostrates the basics of tagging a jet. The PID response
	18	* is retrived for both the TPC. A Hadronic tag is implemented for
	19	* clusters matched to tracks in a jet with a E/P < 0.2
	20	*
	21	* Author: Andrew Castro (UTK) and Joel Mazer (UTK)
	22	*/
	23
	24	#include "AliAnalysisTaskEmcalJetFlavourTagExample.h"
	25
	26	// general ROOT includes
	27	#include <TCanvas.h>
	28	#include <TChain.h>
	29	#include <TClonesArray.h>
	30	#include <TH1F.h>
	31	#include <TH2F.h>
	32	#include <TH3F.h>
	33	#include <THnSparse.h>
	34	#include <TList.h>
	35	#include <TLorentzVector.h>
	36	#include <TParameter.h>
	37	#include <TParticle.h>
	38	#include <TTree.h>
	39	#include <TVector3.h>
	40	#include <TObjArray.h>
	41
	42	// AliROOT includes
	43	#include "AliAODEvent.h"
	44	#include "AliESDEvent.h"
	45	#include "AliAnalysisManager.h"
	46	#include "AliAnalysisTask.h"
	47	#include "AliCentrality.h"
	48	#include "AliEmcalJet.h"
	49	#include "AliAODJet.h"
	50	#include "AliVCluster.h"
	51	#include "AliVTrack.h"
	52	#include <AliVEvent.h>
	53	#include <AliVParticle.h>
	54	#include "AliRhoParameter.h"
	55	#include "AliLog.h"
	56	#include "AliJetContainer.h"
	57	#include "AliParticleContainer.h"
	58	#include "AliClusterContainer.h"
	59	#include "AliEmcalParticle.h"
	60	#include "AliESDCaloCluster.h"
	61	#include <AliESDtrackCuts.h>
	62	#include "AliPID.h"
	63
	64	// event handler (and pico's) includes
65	#include <AliInputEventHandler.h>
66	#include <AliVEventHandler.h>
67	#include "AliESDInputHandler.h"
68	#include "AliPicoTrack.h"
69	#include "AliEventPoolManager.h"
70	#include "AliAODTrack.h"
71	#include "AliESDtrack.h"
72
73	// PID includes
74	#include "AliPIDResponse.h"
75	#include "AliTPCPIDResponse.h"
76	#include "AliESDpid.h"
77
78	// magnetic field includes
79	#include "TGeoGlobalMagField.h"
80	#include "AliMagF.h"
81
a9951dde	82	using std::cout;
	83	using std::endl;
	84
9382bd6e	85	ClassImp(AliAnalysisTaskEmcalJetFlavourTagExample)
	86
	87	//________________________________________________________________________
	88	AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample() :
	89	AliAnalysisTaskEmcalJet("heavyF",kFALSE),
	90	event(0),
	91	fCuts(0),
	92	fPhimin(-10), fPhimax(10),
	93	fEtamin(-0.9), fEtamax(0.9),
	94	fAreacut(0.0),
	95	fJetHIpt(50.0),
	96	fTrackPtCut(2.0),
	97	fTrackEta(0.9),
	98	fesdTrackCuts(0),
	99	fPIDResponse(0x0), fTPCResponse(),
	100	fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
	101	fESD(0), fAOD(0),
	102	fHistJetPhi(0),
	103	fHistCorJetPt(0), fHistJetPt(0),
	104	fHistHighJetPt(0),
	105	fHistnSigElecPt(0),
	106	fHistdEdXvsPt(0),
	107	fHistnJetTrackvnJetClusters(0)
	108	{
	109	// Default constructor.
	110
	111
	112	SetMakeGeneralHistograms(kTRUE);
	113
	114	}
	115
	116	//________________________________________________________________________
	117	AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample(const char *name) :
	118	AliAnalysisTaskEmcalJet(name,kTRUE),
	119	event(0),
	120	fCuts(0),
	121	fPhimin(-10), fPhimax(10),
	122	fEtamin(-0.9), fEtamax(0.9),
	123	fAreacut(0.0),
	124	fJetHIpt(50.0),
	125	fTrackPtCut(2.0),
	126	fTrackEta(0.9),
	127	fesdTrackCuts(0),
	128	fPIDResponse(0x0), fTPCResponse(),
	129	fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
	130	fESD(0), fAOD(0),
	131	fHistJetPhi(0),
	132	fHistCorJetPt(0), fHistJetPt(0),
	133	fHistHighJetPt(0),
	134	fHistnSigElecPt(0),
	135	fHistdEdXvsPt(0),
	136	fHistnJetTrackvnJetClusters(0)
	137	{
	138
	139	SetMakeGeneralHistograms(kTRUE);
	140
	141	DefineInput(0,TChain::Class());
	142	DefineOutput(1, TList::Class());
	143	}
	144
	145	//_______________________________________________________________________
	146	AliAnalysisTaskEmcalJetFlavourTagExample::~AliAnalysisTaskEmcalJetFlavourTagExample()
	147	{
	148	// destructor
149	//
150	if (fOutput) {
151	delete fOutput;
152	fOutput = 0;
153	}
154	}
155
156	//________________________________________________________________________
157	void AliAnalysisTaskEmcalJetFlavourTagExample::UserCreateOutputObjects()
158	{
159	if (! fCreateHisto)
160	return;
161	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
162
163	//fJetsCont = GetJetContainer(0);
164	if(fJetsCont) { //get particles and clusters connected to jets
165	fTracksJetCont = fJetsCont->GetParticleContainer();
166	fCaloClustersJetCont = fJetsCont->GetClusterContainer();
167	}
168	else { //no jets, just analysis tracks and clusters
169	fTracksCont = GetParticleContainer(0);
170	fCaloClustersCont = GetClusterContainer(0);
171	}
172	fTracksCont->SetClassName("AliVTrack");
173	fCaloClustersCont->SetClassName("AliVCluster");
174
175	fHistJetPhi = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2TMath::Pi(),2TMath::Pi());
176	fHistJetPt = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
177	fOutput->Add(fHistJetPhi);
178	fOutput->Add(fHistJetPt);
179
180
181	TString histname;
182
183	fHistCorJetPt = new TH1F("CorrJetPt", "CorrJetPt", 300, -100, 200);
184	fHistnSigElecPt = new TH2F("nsigvsPt(TPC)","nSigmaElectronTPC_v_Pt",60,0,60,40,-10,10);
185	fHistdEdXvsPt = new TH2F("dEdXvsTrackPt","dEdXvsTrackPt",60,0,60,80,0,80);
186	fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
187	fHistHighJetPt = new TH1F("HighestPtJetPerEvent","HighJetPt",80,0,80);
188
189	TString name;
190	TString title;
191
192	fOutput->Add(fHistCorJetPt);
193	fOutput->Add(fHistnSigElecPt);
194	fOutput->Add(fHistdEdXvsPt);
195	fOutput->Add(fHistnJetTrackvnJetClusters);
196	fOutput->Add(fHistHighJetPt);
197
198
199	// **************************** PID ***************************************************
200	// set up PID handler
201	AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
202	AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
203	if(!inputHandler) {
204	AliFatal("Input handler needed");
205	}
206
207	// PID response object
208	//fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
209	// inputHandler->CreatePIDResponse(fIsMC); // needed to create object, why though?
210	fPIDResponse = inputHandler->GetPIDResponse();
211	if (!fPIDResponse) {
212	AliError("PIDResponse object was not created");
213	}
214	// ***************************************************************************************
215
216	PostData(1, fOutput);
217
218	}
219
220	//________________________________________________________
221	void AliAnalysisTaskEmcalJetFlavourTagExample::ExecOnce()
222	{
223	// Initialize the analysis
224	AliAnalysisTaskEmcalJet::ExecOnce();
225
226	if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
227	if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
228	if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
229
230
231	} // end of ExecOnce
232
233	//________________________________________________________________________
234	Bool_t AliAnalysisTaskEmcalJetFlavourTagExample::Run()
235	{
236	// check to see if we have any tracks
237	if (!fTracks) return kTRUE;
238	if (!fJets) return kTRUE;
239
240	// what kind of event do we have: AOD or ESD?
241	Bool_t useAOD;
242	if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
243	else useAOD = kFALSE;
244
245	// if we have ESD event, set up ESD object
246	if(!useAOD){
247	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
248	if (!fESD) {
249	AliError(Form("ERROR: fESD not available\n"));
250	return kTRUE;
251	}
252	}
253
254	// if we have AOD event, set up AOD object
255	if(useAOD){
256	fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
257	if(!fAOD) {
258	AliError(Form("ERROR: fAOD not available\n"));
259	return kTRUE;
260	}
261	}
262
263	// get magnetic field info for DCA
264	Double_t MagF = fESD->GetMagneticField();
265	Double_t MagSign = 1.0;
266	if(MagF<0)MagSign = -1.0;
267	// set magnetic field
268	if (!TGeoGlobalMagField::Instance()->GetField()) {
269	AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
270	TGeoGlobalMagField::Instance()->SetField(field);
271	}
272
273	// get vertex information
274	Double_t fvertex[3]={0,0,0};
275	InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
276	//Double_t zVtx=fvertex[2];
277
278	// create pointer to list of input event
279	TList *list = InputEvent()->GetList();
280	if(!list) {
281	AliError(Form("ERROR: list not attached\n"));
282	return kTRUE;
283	}
284
285	// background density
286	fRhoVal = fRho->GetVal();
287
288	// initialize TClonesArray pointers to jets and tracks
289	TClonesArray *jets = 0;
290	//TClonesArray *tracks = 0;
291	//TClonesArray *clusters = 0;
292
293	// get Jets object
294	jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
295	if(!jets){
296	AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
297	return kTRUE;
298	} // verify existence of jets
299
300	// get number of jets and tracks
301	const Int_t Njets = jets->GetEntries();
302	if(Njets<1) return kTRUE;
303
304	if (fTracksCont) {
305	AliVTrack track = static_cast<AliVTrack>(fTracksCont->GetNextAcceptParticle(0));
306	while(track) {
307	track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
308	}
309	}
310	if (fCaloClustersCont) {
311	AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
312	while(cluster) {
313	TLorentzVector nPart;
314	cluster->GetMomentum(nPart, fVertex);
315	cluster = fCaloClustersCont->GetNextAcceptCluster();
316	}
317	}
318
319	// Start Jet Analysis
320	// initialize jet parameters
321	Int_t ijethi=-1;
322	Double_t highestjetpt=0.0;
323
324	// loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
325	for (Int_t ijet = 0; ijet < Njets; ijet++){
326	// get our jets
327	AliEmcalJet jet = static_cast<AliEmcalJet>(jets->At(ijet));
328	if (!jet) continue;
329
330	// apply jet cuts
331	if(!AcceptMyJet(jet)) continue;
332
333	if(highestjetpt<jet->Pt()){
334	ijethi=ijet;
335	highestjetpt=jet->Pt();
336	}
337	} // end of looping over jets
338
339	fHistHighJetPt->Fill(ijethi);
340
341	// **********************************************************************
342	// JET LOOP
343	// **********************************************************************
344
345	// loop over jets in the event and make appropriate cuts
346	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
347	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(iJets));
348	if (!jet) // see if we have a jet
349	continue;
350
351	// apply jet cuts
352	if(!AcceptMyJet(jet)) continue;
353
354	//AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
355	//jet->AddFlavourTag(tag);
356
357
358	Int_t JetClusters = jet->GetNumberOfClusters();
359	Int_t JetTracks = jet -> GetNumberOfTracks();
360	fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
361	// Initializations and Calculations
362	Double_t jetPt = -500; // initialize corr jet pt LOCAL
363	Double_t jetarea = -500; // initialize jet area
364	jetPt = jet->Pt() - jetarea*fRhoVal; // semi-corrected pT of jet from GLOBAL rho value
365	fHistCorJetPt->Fill(jetPt);
366
367	Bool_t bkgrnd1 = kFALSE;
368
369	if(jet->Pt() > fJetHIpt) {
370	if(!fTracksCont \|\| !fCaloClustersCont) continue;
371
372	//******************************Cluster Matched To Closest Track
373	//**************************************************************
374	Int_t NumbTrackContainer = -999;
375	NumbTrackContainer = fTracksCont->GetNParticles();
376	for(int iTracks = 0; iTracks <= NumbTrackContainer; iTracks++){
377	AliVTrack AcceptedTrack =static_cast<AliVTrack>(fTracksCont->GetParticle(iTracks));
378	if(!AcceptedTrack){
379	AliError(Form("Couldn't get AliVTrack Container %d\n", iTracks));
380	continue;
381	}
382	if(!IsJetTrack(jet,iTracks,kFALSE))continue;
383	//Get matched cluster
384	Int_t emc1 = AcceptedTrack->GetEMCALcluster();
385
386	Double_t acceptTrackP = AcceptedTrack->P();
387	Double_t acceptTrackPt = AcceptedTrack->Pt();
388	Double_t nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(AcceptedTrack,AliPID::kElectron);
389	fHistnSigElecPt->Fill(acceptTrackPt,nSigmaElectron_TPC);
390
391	AliESDtrack ESDacceptedTrack = static_cast<AliESDtrack>(AcceptedTrack);
392	if(!ESDacceptedTrack){
393	AliError(Form("Couldn't get AliESDTrack %d\n", iTracks));
394	continue;
395	}
396	Double_t dEdx = AcceptedTrack->GetTPCsignal();
397	fHistdEdXvsPt->Fill(acceptTrackPt,dEdx);
398	if(fCaloClustersCont && emc1>=0) {
399	AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
400	if(!clusMatch){
401	AliError(Form("Couldn't get matched AliVCluster %d\n", emc1));
402	continue;
403	}
404
405	Double_t mClusterE = clusMatch->E();
406	Double_t EovP_mc = -999;
407	EovP_mc = mClusterE/acceptTrackP;
408
409	if(EovP_mc < 0.2) bkgrnd1 = kTRUE;
410	}
411
412	} //loop over tracks for matching to closest cluster
413
414
415
416
417	} // highest pt jet cut
418
419
420	if(bkgrnd1 == kTRUE) {
421	AliEmcalJet::EFlavourTag tag=AliEmcalJet::kBckgrd1;
422	jet->AddFlavourTag(tag);
423	}
424
425	} // LOOP over JETS in event
426
427
428	return kTRUE;
429
430	}
431	//________________________________________________________________________
432	void AliAnalysisTaskEmcalJetFlavourTagExample::Terminate(Option_t *)
433	{
434	cout<<"###########################"<<endl;
435	cout<<"#### Task Finished ####"<<endl;
436	cout<<"###########################"<<endl;
437	cout<<"###########################"<<endl;
438	} // end of terminate
439
440	//________________________________________________________________________
441	Int_t AliAnalysisTaskEmcalJetFlavourTagExample::AcceptMyJet(AliEmcalJet *jet) {
442	//applies all jet cuts except pt
443	if ((jet->Phi()<fPhimin)\|\|(jet->Phi()>fPhimax)) return 0;
444	if ((jet->Eta()<fEtamin)\|\|(jet->Eta()>fEtamax)) return 0;
445	if (jet->Area()<fAreacut) return 0;
446	// prevents 0 area jets from sneaking by when area cut == 0
447	if (jet->Area()==0) return 0;
448	//exclude jets with extremely high pt tracks which are likely misreconstructed
449	if(jet->MaxTrackPt()>100) return 0;
450
451	//passed all above cuts
452	return 1;
453	}
454
455
456
457
458