[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskQGSep.cxx

#include <TChain.h>
#include <TList.h>

#include <TTree.h>
#include <TFile.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TCanvas.h>

#include "AliVEvent.h"
#include "AliVParticle.h"

#include "AliESDEvent.h"
#include "AliESDtrack.h"

#include "AliAODEvent.h"
#include "AliAODTrack.h"

#include "AliMCEvent.h"
#include "AliMCParticle.h"
#include "AliAnalysisManager.h"
#include "AliMCEvent.h"
#include "TParticle.h"
#include "AliStack.h"

#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODHandler.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
#include "AliAODMCParticle.h"

#include "AliAnalysisTaskQGSep.h"
#include "AliAnalysisHelperJetTasks.h"
ClassImp(AliAnalysisTaskQGSep)

//________________________________________________________________________
AliAnalysisTaskQGSep::AliAnalysisTaskQGSep(const char *name)
  : AliAnalysisTaskSE(name),
    fBranchRec("jets"),
    fUseMC(kFALSE),
    fUseAOD(kFALSE),
    fXsection(0),
    fWeight(0),
    fMyAODEvent(0),
    fOutputList(0),
    fpHistPtAvEQ(0),
    fpHistPtAvEG(0),
    fpHistDrEQ(0),
    fpHistDrEG(0),
    fpHistDrE(0),
    fpHistPtAvE(0),
    fpHistDrE3(0),
    fpHistPtAvE3(0)

{
  // Constructor

  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskQGSep::UserCreateOutputObjects()
{
  // Create histograms
  // Called once

  fOutputList = new TList();

  if(fUseMC){

    //histos for quarks
    fpHistPtAvEQ = new TProfile("fpHistPtAvEQ", "", 100, 0, 100, 0, 10);
    fOutputList->Add(fpHistPtAvEQ);
    fpHistDrEQ = new TProfile("fpHistDrEQ", "", 100, 0, 100, 0, 0.5);
    fOutputList->Add(fpHistDrEQ);

    //histos for gluons
    fpHistPtAvEG = new TProfile("fpHistPtAvEG", "", 100, 0, 100, 0., 10.);
    fOutputList->Add(fpHistPtAvEG);
    fpHistDrEG = new TProfile("fpHistDrEG", "", 100, 0, 100, 0, 0.5);
    fOutputList->Add(fpHistDrEG);
  }
  
  //histos for full spetra, no separation
  fpHistDrE = new TProfile("fpHistDrE", "", 100, 0, 100, 0, 0.5);
  fOutputList->Add(fpHistDrE);
  fpHistPtAvE = new TProfile("fpHistPtAvE", "", 100, 0, 100, 0., 10.);
  fOutputList->Add(fpHistPtAvE);
  fpHistDrE3 = new TProfile("fpHistDrE3", "", 100, 0, 100, 0, 0.5);
  fOutputList->Add(fpHistDrE3);
  fpHistPtAvE3 = new TProfile("fpHistPtAvE3", "", 100, 0, 100, 0., 10.);
  fOutputList->Add(fpHistPtAvE3);

  for(Int_t i = 0; i < fOutputList->GetEntries(); i++){
    TH1 * h = dynamic_cast<TH1*>(fOutputList->At(i));
    if(h) h->Sumw2();
  }

  if(fDebug)Printf("~# QGSep User objects created");
}
  
//__________________________________________________________________
Bool_t AliAnalysisTaskQGSep::Notify()
{
  
  if(fUseMC){
    //
    // Implemented Notify() to read the cross sections
    // and number of trials from pyxsec.root
    // 
    
    
    TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
    Float_t xsection = 0;
    Float_t ftrials  = 1;
    
    Float_t fAvgTrials = 1;
    if(tree){
      TFile *curfile = tree->GetCurrentFile();
      if (!curfile) {
	Error("Notify","No current file");
	return kFALSE;
      }
      AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
      // construct a poor man average trials 
      Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
      if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
    }  
  
    if(xsection>0){
      fXsection  = xsection;
      fWeight = fXsection/fAvgTrials;
    }
    else fWeight = 1;
    return kTRUE;
  }
  
  else{
    fWeight = 1;
    return kTRUE;
  }
  
}

//________________________________________________________________________
void AliAnalysisTaskQGSep::UserExec(Option_t *)
{
  // Main loop
  // Called for each event
  
  if(fUseAOD){
    AliVEvent *event = InputEvent();
    if (!event) {
      Error("UserExec", "Could not retrieve event");
      return;
    }
  
  fMyAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
  }

  else{
        //  assume that the AOD is in the general output...
    fMyAODEvent  = AODEvent();
    if(!fMyAODEvent){
      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
      return;
    }    
  }

  if (fMyAODEvent){
    if(fUseMC)
      LoopAODMC();
    else
      LoopAOD();
  }
	
  // Post output data.
  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskQGSep::Terminate(Option_t*)
{
  // Draw result to the screen
  // Called once at the end of the query
  
  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    Error("Terminate","fOutputList not available");
    return;
  }
  
}

//________________________________________________________________________
void AliAnalysisTaskQGSep::LoopAOD(){
  AliAODJet recJets[4];
  AliAODJet jets[4];
  AliAODJet rJets[4];
  Int_t nRecJets = 0;
  
  //array of reconstructed jets from the AOD input
  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(fBranchRec.Data()));
  if(!aodRecJets){
    PostData(1, fOutputList);
    return;
  }
   
  // reconstructed jets
  nRecJets = aodRecJets->GetEntries(); 
  if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
  nRecJets = TMath::Min(nRecJets, 4);
  for(int ir = 0;ir < nRecJets;++ir)
    {
      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
      if(!tmp)continue;
      jets[ir] = *tmp;
    }
  
  Int_t counter = 0;
  Int_t tag = 0;
  Int_t nGenSel = 0;
  
  TLorentzVector v[4];
  Double_t eSum = 0.;
  Double_t pxSum = 0.;
  Double_t pySum = 0.;
  Double_t pzSum = 0.;
  
  for(Int_t i = 0; i < nRecJets; i++)
    {
      if(nRecJets == 1)
	{
	  rJets[nGenSel] = jets[i];
	  v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
	  eSum += jets[i].E();
	  pxSum += jets[i].Px();
	  pySum += jets[i].Py();
	  pzSum += jets[i].Pz();		    
	  nGenSel++;
	}
      else
	{
	  counter = 0;
	  tag = 0;
	  for(Int_t j = 0; j < nRecJets; j++)
	    {
	      if(i!=j)
		{
		  Double_t dRij = jets[i].DeltaR(&jets[j]);
		  counter++;
		  if(dRij > 2*0.4) tag++;
		}
	    }
	  if(counter!=0)
	    {
	      if(tag/counter == 1)
		{
		  rJets[nGenSel] = jets[i];
		  v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
		  eSum += jets[i].E();
		  pxSum += jets[i].Px();
		  pySum += jets[i].Py();
		  pzSum += jets[i].Pz();		    
		  nGenSel++;
		}
	    }
	}
    }
   
  nRecJets = nGenSel;
  
  if(nRecJets == 0){
    PostData(1, fOutputList);
    return;
  }

  TLorentzVector vB;
  vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
  Double_t e[4];
  for(Int_t i = 0; i < nRecJets; i++){
    v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
    e[i] = v[i].E();
  }

  Int_t idxj[4];
  TMath::Sort(TMath::Min(nRecJets,4), e, idxj);
  for(Int_t i = 0; i < nRecJets; i++){
    recJets[i] = rJets[idxj[i]];
  }
  
   
  for(Int_t iJ = 0; iJ < nRecJets; iJ++){
    Double_t pTsum = 0.;
    TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
    if(!tra) continue;
    Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
    Double_t dR[1000];
    for(Int_t iT = 0; iT < nAODtracks; iT++){
      AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
      if(!jetTrack) continue;
      pTsum += jetTrack->Pt();
      dR[iT] = recJets[iJ].DeltaR(jetTrack);
    }
     
 
    fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
           
    if(iJ > 1) //fill mulit-jet histo
      fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
      
    Int_t idxAOD[1000];
    TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);
     
    Double_t pTsum90Inv=0.;
    for(Int_t iT = 0; iT < nAODtracks; iT++){
      AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
      if(!track) continue;
      pTsum90Inv += track->Pt();
      if(pTsum90Inv >= 0.9*pTsum){
	Double_t deltaR = recJets[iJ].DeltaR(track);

	fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);

	if(iJ > 1)
	  fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);

	break;
      }
    }    
  }
}
 
//__________________________________________________________________
void AliAnalysisTaskQGSep::LoopAODMC(){
  
  AliAODJet recJets[4];
  AliAODJet jets[4];
  AliAODJet rJets[4];
  Int_t nRecJets = 0;
  
  //array of reconstructed jets from the AOD input
  TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(fBranchRec.Data()));
  if(!aodRecJets){
    PostData(1, fOutputList);
    return;
  }
  
  // reconstructed jets
  nRecJets = aodRecJets->GetEntries(); 
  if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
  nRecJets = TMath::Min(nRecJets, 4);
  for(int ir = 0;ir < nRecJets;++ir)
    {
      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
      if(!tmp)continue;
      jets[ir] = *tmp;
    }
  
  Int_t counter = 0;
  Int_t tag = 0;
  Int_t nGenSel = 0;
  
  TLorentzVector v[4];
  Double_t eSum = 0.;
  Double_t pxSum = 0.;
  Double_t pySum = 0.;
  Double_t pzSum = 0.;
  
  for(Int_t i = 0; i < nRecJets; i++)
    {
      if(nRecJets == 1)
	{
	  rJets[nGenSel] = jets[i];
	  v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
	  eSum += jets[i].E();
	  pxSum += jets[i].Px();
	  pySum += jets[i].Py();
	  pzSum += jets[i].Pz();		    
	  nGenSel++;
	}
      else
	{
	  counter = 0;
	  tag = 0;
	  for(Int_t j = 0; j < nRecJets; j++)
	    {
	      if(i!=j)
		{
		  Double_t dRij = jets[i].DeltaR(&jets[j]);
		  counter++;
		  if(dRij > 2*0.4) tag++;
		}
	    }
	  if(counter!=0)
	    {
	      if(tag/counter == 1)
		{
		  rJets[nGenSel] = jets[i];
		  v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
		  eSum += jets[i].E();
		  pxSum += jets[i].Px();
		  pySum += jets[i].Py();
		  pzSum += jets[i].Pz();		    
		  nGenSel++;
		}
	    }
	}
    }
  
  nRecJets = nGenSel;
  
  if(nRecJets == 0){
    PostData(1, fOutputList);
    return;
   }
   
  TLorentzVector vB;
  vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
  Double_t e[4];
  for(Int_t i = 0; i < nRecJets; i++){
    v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
    e[i] = v[i].E();
  }

  Int_t idxj[4];
  TMath::Sort(TMath::Min(nRecJets,4), e, idxj);
  for(Int_t i = 0; i < nRecJets; i++){
    recJets[i] = rJets[idxj[i]];
  }
  
  Int_t nMCtracks = 0;
  TClonesArray *tca = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
  if(!tca) return;
  nMCtracks = TMath::Min(tca->GetEntries(), 10000);
  
  //sort AliAODMCParticles in pT
  Double_t pTMC[10000];
  for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
    AliAODMCParticle *partMC = dynamic_cast<AliAODMCParticle*>(tca->At(iMC));
    if(!partMC) continue;
    pTMC[iMC]=partMC->Pt();
  }
  
  Int_t idxMC[10000];
  TMath::Sort(nMCtracks, pTMC, idxMC);
  
  
  Int_t flagQ[4] = {0}, flagG[4] = {0};   
  for(Int_t iJ = 0; iJ < nRecJets; iJ++){
    //flag jet as q/g
    
    //look for highest momentum parton in the jet cone
    for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
      AliAODMCParticle *partMC = dynamic_cast<AliAODMCParticle*>(tca->At(idxMC[iMC]));
      if(!partMC) continue;
      Double_t r = recJets[iJ].DeltaR(partMC);
      if(r < 0.4){
	if(TMath::Abs(partMC->GetPdgCode()) < 9)
	  flagQ[iJ] = 1;
	if(TMath::Abs(partMC->GetPdgCode()) == 21)
	  flagG[iJ] = 1;
	break;
      }
    }

    Double_t pTsum = 0.;
    TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
    if(!tra) continue;
    Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
    Double_t dR[1000];
    for(Int_t iT = 0; iT < nAODtracks; iT++){
      AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
      if(!jetTrack) continue;
      pTsum += jetTrack->Pt();
      dR[iT] = recJets[iJ].DeltaR(jetTrack);
    }

    fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
    if(flagQ[iJ] == 1) fpHistPtAvEQ->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
    if(flagG[iJ] == 1) fpHistPtAvEG->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
        
    if(iJ > 1) //fill mulit-jet histo
      fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
      
    Int_t idxAOD[1000];
    TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);

    Double_t pTsum90Inv=0.;
    for(Int_t iT = 0; iT < nAODtracks; iT++){
      AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
      if(!track) continue;
      pTsum90Inv += track->Pt();
      if(pTsum90Inv >= 0.9*pTsum){
	Double_t deltaR = recJets[iJ].DeltaR(track);

	fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);
	if(flagQ[iJ] == 1) fpHistDrEQ->Fill(recJets[iJ].E(), deltaR, fWeight);
	if(flagG[iJ] == 1) fpHistDrEG->Fill(recJets[iJ].E(), deltaR, fWeight);

	if(iJ > 1)
	  fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);

	break;
      }
    }    
  }
}
Commit	Line	Data
cc636095	1	#include <TChain.h>
	2	#include <TList.h>
	3
	4	#include <TTree.h>
	5	#include <TFile.h>
	6	#include <TH1.h>
	7	#include <TH1F.h>
	8	#include <TH2F.h>
	9	#include <TProfile.h>
	10	#include <TCanvas.h>
	11
	12	#include "AliVEvent.h"
	13	#include "AliVParticle.h"
	14
	15	#include "AliESDEvent.h"
	16	#include "AliESDtrack.h"
	17
	18	#include "AliAODEvent.h"
	19	#include "AliAODTrack.h"
	20
	21	#include "AliMCEvent.h"
	22	#include "AliMCParticle.h"
	23	#include "AliAnalysisManager.h"
	24	#include "AliMCEvent.h"
	25	#include "TParticle.h"
	26	#include "AliStack.h"
	27
	28	#include "AliAODEvent.h"
	29	#include "AliAODVertex.h"
	30	#include "AliAODHandler.h"
	31	#include "AliAODTrack.h"
	32	#include "AliAODJet.h"
	33	#include "AliAODMCParticle.h"
	34
	35	#include "AliAnalysisTaskQGSep.h"
	36	#include "AliAnalysisHelperJetTasks.h"
	37	ClassImp(AliAnalysisTaskQGSep)
	38
	39	//________________________________________________________________________
	40	AliAnalysisTaskQGSep::AliAnalysisTaskQGSep(const char *name)
	41	: AliAnalysisTaskSE(name),
	42	fBranchRec("jets"),
	43	fUseMC(kFALSE),
	44	fUseAOD(kFALSE),
	45	fXsection(0),
	46	fWeight(0),
	47	fMyAODEvent(0),
	48	fOutputList(0),
	49	fpHistPtAvEQ(0),
	50	fpHistPtAvEG(0),
	51	fpHistDrEQ(0),
	52	fpHistDrEG(0),
	53	fpHistDrE(0),
	54	fpHistPtAvE(0),
	55	fpHistDrE3(0),
	56	fpHistPtAvE3(0)
	57
	58	{
	59	// Constructor
	60
	61	DefineOutput(1, TList::Class());
	62	}
	63
	64	//________________________________________________________________________
65	void AliAnalysisTaskQGSep::UserCreateOutputObjects()
66	{
67	// Create histograms
68	// Called once
69
70	fOutputList = new TList();
71
72	if(fUseMC){
73
74	//histos for quarks
75	fpHistPtAvEQ = new TProfile("fpHistPtAvEQ", "", 100, 0, 100, 0, 10);
76	fOutputList->Add(fpHistPtAvEQ);
77	fpHistDrEQ = new TProfile("fpHistDrEQ", "", 100, 0, 100, 0, 0.5);
78	fOutputList->Add(fpHistDrEQ);
79
80	//histos for gluons
81	fpHistPtAvEG = new TProfile("fpHistPtAvEG", "", 100, 0, 100, 0., 10.);
82	fOutputList->Add(fpHistPtAvEG);
83	fpHistDrEG = new TProfile("fpHistDrEG", "", 100, 0, 100, 0, 0.5);
84	fOutputList->Add(fpHistDrEG);
85	}
86
87	//histos for full spetra, no separation
88	fpHistDrE = new TProfile("fpHistDrE", "", 100, 0, 100, 0, 0.5);
89	fOutputList->Add(fpHistDrE);
90	fpHistPtAvE = new TProfile("fpHistPtAvE", "", 100, 0, 100, 0., 10.);
91	fOutputList->Add(fpHistPtAvE);
92	fpHistDrE3 = new TProfile("fpHistDrE3", "", 100, 0, 100, 0, 0.5);
93	fOutputList->Add(fpHistDrE3);
94	fpHistPtAvE3 = new TProfile("fpHistPtAvE3", "", 100, 0, 100, 0., 10.);
95	fOutputList->Add(fpHistPtAvE3);
96
97	for(Int_t i = 0; i < fOutputList->GetEntries(); i++){
98	TH1 * h = dynamic_cast<TH1*>(fOutputList->At(i));
99	if(h) h->Sumw2();
100	}
101
102	if(fDebug)Printf("~# QGSep User objects created");
103	}
104
105	//__________________________________________________________________
106	Bool_t AliAnalysisTaskQGSep::Notify()
107	{
108
109	if(fUseMC){
110	//
111	// Implemented Notify() to read the cross sections
112	// and number of trials from pyxsec.root
113	//
114
115
116	TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
117	Float_t xsection = 0;
118	Float_t ftrials = 1;
119
120	Float_t fAvgTrials = 1;
121	if(tree){
122	TFile *curfile = tree->GetCurrentFile();
123	if (!curfile) {
124	Error("Notify","No current file");
125	return kFALSE;
126	}
127	AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
128	// construct a poor man average trials
129	Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
130	if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; // CKB take this into account for normalisation
131	}
132
133	if(xsection>0){
134	fXsection = xsection;
135	fWeight = fXsection/fAvgTrials;
136	}
137	else fWeight = 1;
138	return kTRUE;
139	}
140
141	else{
142	fWeight = 1;
143	return kTRUE;
144	}
145
146	}
147
148	//________________________________________________________________________
149	void AliAnalysisTaskQGSep::UserExec(Option_t *)
150	{
151	// Main loop
152	// Called for each event
153
154	if(fUseAOD){
155	AliVEvent *event = InputEvent();
156	if (!event) {
157	Error("UserExec", "Could not retrieve event");
158	return;
159	}
160
161	fMyAODEvent = dynamic_cast<AliAODEvent*> (InputEvent());
162	}
163
164	else{
165	// assume that the AOD is in the general output...
166	fMyAODEvent = AODEvent();
167	if(!fMyAODEvent){
168	Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
169	return;
170	}
171	}
172
173	if (fMyAODEvent){
174	if(fUseMC)
175	LoopAODMC();
176	else
177	LoopAOD();
178	}
179
180	// Post output data.
181	PostData(1, fOutputList);
182	}
183
184	//________________________________________________________________________
185	void AliAnalysisTaskQGSep::Terminate(Option_t*)
186	{
187	// Draw result to the screen
188	// Called once at the end of the query
189
190	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
191	if (!fOutputList) {
192	Error("Terminate","fOutputList not available");
193	return;
194	}
195
196	}
197
198	//________________________________________________________________________
199	void AliAnalysisTaskQGSep::LoopAOD(){
200	AliAODJet recJets[4];
201	AliAODJet jets[4];
202	AliAODJet rJets[4];
203	Int_t nRecJets = 0;
204
205	//array of reconstructed jets from the AOD input
206	TClonesArray aodRecJets = dynamic_cast<TClonesArray>(fMyAODEvent->FindListObject(fBranchRec.Data()));
207	if(!aodRecJets){
208	PostData(1, fOutputList);
209	return;
210	}
211
212	// reconstructed jets
213	nRecJets = aodRecJets->GetEntries();
214	if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
215	nRecJets = TMath::Min(nRecJets, 4);
216	for(int ir = 0;ir < nRecJets;++ir)
217	{
218	AliAODJet tmp = dynamic_cast<AliAODJet>(aodRecJets->At(ir));
219	if(!tmp)continue;
220	jets[ir] = *tmp;
221	}
222
223	Int_t counter = 0;
224	Int_t tag = 0;
225	Int_t nGenSel = 0;
226
227	TLorentzVector v[4];
228	Double_t eSum = 0.;
229	Double_t pxSum = 0.;
230	Double_t pySum = 0.;
231	Double_t pzSum = 0.;
232
233	for(Int_t i = 0; i < nRecJets; i++)
234	{
235	if(nRecJets == 1)
236	{
237	rJets[nGenSel] = jets[i];
238	v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
239	eSum += jets[i].E();
240	pxSum += jets[i].Px();
241	pySum += jets[i].Py();
242	pzSum += jets[i].Pz();
243	nGenSel++;
244	}
245	else
246	{
247	counter = 0;
248	tag = 0;
249	for(Int_t j = 0; j < nRecJets; j++)
250	{
251	if(i!=j)
252	{
253	Double_t dRij = jets[i].DeltaR(&jets[j]);
254	counter++;
255	if(dRij > 2*0.4) tag++;
256	}
257	}
258	if(counter!=0)
259	{
260	if(tag/counter == 1)
261	{
262	rJets[nGenSel] = jets[i];
263	v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
264	eSum += jets[i].E();
265	pxSum += jets[i].Px();
266	pySum += jets[i].Py();
267	pzSum += jets[i].Pz();
268	nGenSel++;
269	}
270	}
271	}
272	}
273
274	nRecJets = nGenSel;
275
276	if(nRecJets == 0){
277	PostData(1, fOutputList);
278	return;
279	}
280
281	TLorentzVector vB;
282	vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
283	Double_t e[4];
284	for(Int_t i = 0; i < nRecJets; i++){
285	v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
286	e[i] = v[i].E();
287	}
288
289	Int_t idxj[4];
e855f5c5	290	TMath::Sort(TMath::Min(nRecJets,4), e, idxj);
cc636095	291	for(Int_t i = 0; i < nRecJets; i++){
	292	recJets[i] = rJets[idxj[i]];
	293	}
	294
	295
	296	for(Int_t iJ = 0; iJ < nRecJets; iJ++){
	297	Double_t pTsum = 0.;
	298	TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
	299	if(!tra) continue;
	300	Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
	301	Double_t dR[1000];
	302	for(Int_t iT = 0; iT < nAODtracks; iT++){
	303	AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
	304	if(!jetTrack) continue;
	305	pTsum += jetTrack->Pt();
	306	dR[iT] = recJets[iJ].DeltaR(jetTrack);
	307	}
	308
	309
	310	fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	311
	312	if(iJ > 1) //fill mulit-jet histo
	313	fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	314
	315	Int_t idxAOD[1000];
	316	TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);
	317
	318	Double_t pTsum90Inv=0.;
	319	for(Int_t iT = 0; iT < nAODtracks; iT++){
	320	AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
	321	if(!track) continue;
	322	pTsum90Inv += track->Pt();
	323	if(pTsum90Inv >= 0.9*pTsum){
	324	Double_t deltaR = recJets[iJ].DeltaR(track);
	325
	326	fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);
	327
	328	if(iJ > 1)
	329	fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);
	330
	331	break;
	332	}
	333	}
	334	}
	335	}
	336
	337	//__________________________________________________________________
	338	void AliAnalysisTaskQGSep::LoopAODMC(){
	339
	340	AliAODJet recJets[4];
	341	AliAODJet jets[4];
	342	AliAODJet rJets[4];
	343	Int_t nRecJets = 0;
	344
	345	//array of reconstructed jets from the AOD input
	346	TClonesArray aodRecJets = dynamic_cast<TClonesArray>(fMyAODEvent->FindListObject(fBranchRec.Data()));
	347	if(!aodRecJets){
	348	PostData(1, fOutputList);
	349	return;
	350	}
	351
	352	// reconstructed jets
	353	nRecJets = aodRecJets->GetEntries();
	354	if(fDebug)Printf("--- Jets found in bRec: %d", nRecJets);
355	nRecJets = TMath::Min(nRecJets, 4);
356	for(int ir = 0;ir < nRecJets;++ir)
357	{
358	AliAODJet tmp = dynamic_cast<AliAODJet>(aodRecJets->At(ir));
359	if(!tmp)continue;
360	jets[ir] = *tmp;
361	}
362
363	Int_t counter = 0;
364	Int_t tag = 0;
365	Int_t nGenSel = 0;
366
367	TLorentzVector v[4];
368	Double_t eSum = 0.;
369	Double_t pxSum = 0.;
370	Double_t pySum = 0.;
371	Double_t pzSum = 0.;
372
373	for(Int_t i = 0; i < nRecJets; i++)
374	{
375	if(nRecJets == 1)
376	{
377	rJets[nGenSel] = jets[i];
378	v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
379	eSum += jets[i].E();
380	pxSum += jets[i].Px();
381	pySum += jets[i].Py();
382	pzSum += jets[i].Pz();
383	nGenSel++;
384	}
385	else
386	{
387	counter = 0;
388	tag = 0;
389	for(Int_t j = 0; j < nRecJets; j++)
390	{
391	if(i!=j)
392	{
393	Double_t dRij = jets[i].DeltaR(&jets[j]);
394	counter++;
395	if(dRij > 2*0.4) tag++;
396	}
397	}
398	if(counter!=0)
399	{
400	if(tag/counter == 1)
401	{
402	rJets[nGenSel] = jets[i];
403	v[nGenSel].SetPxPyPzE(jets[i].Px(), jets[i].Py(), jets[i].Pz(), jets[i].E());
404	eSum += jets[i].E();
405	pxSum += jets[i].Px();
406	pySum += jets[i].Py();
407	pzSum += jets[i].Pz();
408	nGenSel++;
409	}
410	}
411	}
412	}
413
414	nRecJets = nGenSel;
415
416	if(nRecJets == 0){
417	PostData(1, fOutputList);
418	return;
419	}
420
421	TLorentzVector vB;
422	vB.SetPxPyPzE(pxSum, pySum, pzSum, eSum);
423	Double_t e[4];
424	for(Int_t i = 0; i < nRecJets; i++){
425	v[i].Boost(-vB.Px()/vB.E(),-vB.Py()/vB.E(),-vB.Pz()/vB.E());
426	e[i] = v[i].E();
427	}
428
429	Int_t idxj[4];
e855f5c5	430	TMath::Sort(TMath::Min(nRecJets,4), e, idxj);
cc636095	431	for(Int_t i = 0; i < nRecJets; i++){
	432	recJets[i] = rJets[idxj[i]];
	433	}
	434
	435	Int_t nMCtracks = 0;
	436	TClonesArray tca = dynamic_cast<TClonesArray>(fMyAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
	437	if(!tca) return;
	438	nMCtracks = TMath::Min(tca->GetEntries(), 10000);
	439
	440	//sort AliAODMCParticles in pT
	441	Double_t pTMC[10000];
	442	for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
	443	AliAODMCParticle partMC = dynamic_cast<AliAODMCParticle>(tca->At(iMC));
	444	if(!partMC) continue;
	445	pTMC[iMC]=partMC->Pt();
	446	}
	447
	448	Int_t idxMC[10000];
	449	TMath::Sort(nMCtracks, pTMC, idxMC);
	450
	451
e855f5c5	452	Int_t flagQ[4] = {0}, flagG[4] = {0};
cc636095	453	for(Int_t iJ = 0; iJ < nRecJets; iJ++){
	454	//flag jet as q/g
	455
	456	//look for highest momentum parton in the jet cone
	457	for(Int_t iMC = 0; iMC < nMCtracks; iMC++){
	458	AliAODMCParticle partMC = dynamic_cast<AliAODMCParticle>(tca->At(idxMC[iMC]));
	459	if(!partMC) continue;
	460	Double_t r = recJets[iJ].DeltaR(partMC);
	461	if(r < 0.4){
	462	if(TMath::Abs(partMC->GetPdgCode()) < 9)
	463	flagQ[iJ] = 1;
	464	if(TMath::Abs(partMC->GetPdgCode()) == 21)
	465	flagG[iJ] = 1;
	466	break;
	467	}
	468	}
	469
	470	Double_t pTsum = 0.;
	471	TRefArray * tra = dynamic_cast<TRefArray*>(recJets[iJ].GetRefTracks());
	472	if(!tra) continue;
	473	Int_t nAODtracks = TMath::Min(1000, tra->GetEntries());
	474	Double_t dR[1000];
	475	for(Int_t iT = 0; iT < nAODtracks; iT++){
	476	AliAODTrack * jetTrack = dynamic_cast<AliAODTrack*>(tra->At(iT));
	477	if(!jetTrack) continue;
	478	pTsum += jetTrack->Pt();
	479	dR[iT] = recJets[iJ].DeltaR(jetTrack);
	480	}
	481
	482	fpHistPtAvE->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	483	if(flagQ[iJ] == 1) fpHistPtAvEQ->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	484	if(flagG[iJ] == 1) fpHistPtAvEG->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	485
	486	if(iJ > 1) //fill mulit-jet histo
	487	fpHistPtAvE3->Fill(recJets[iJ].E(), (Double_t)pTsum/(Double_t)nAODtracks, fWeight);
	488
	489	Int_t idxAOD[1000];
	490	TMath::Sort(nAODtracks, dR, idxAOD, kFALSE);
	491
	492	Double_t pTsum90Inv=0.;
	493	for(Int_t iT = 0; iT < nAODtracks; iT++){
	494	AliAODTrack * track = dynamic_cast<AliAODTrack*>(tra->At(idxAOD[iT]));
	495	if(!track) continue;
	496	pTsum90Inv += track->Pt();
	497	if(pTsum90Inv >= 0.9*pTsum){
	498	Double_t deltaR = recJets[iJ].DeltaR(track);
	499
	500	fpHistDrE->Fill(recJets[iJ].E(), deltaR, fWeight);
	501	if(flagQ[iJ] == 1) fpHistDrEQ->Fill(recJets[iJ].E(), deltaR, fWeight);
	502	if(flagG[iJ] == 1) fpHistDrEG->Fill(recJets[iJ].E(), deltaR, fWeight);
	503
	504	if(iJ > 1)
	505	fpHistDrE3->Fill(recJets[iJ].E(), deltaR, fWeight);
	506
	507	break;
	508	}
	509	}
	510	}
	511	}