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

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "TList.h"
#include "TParticle.h"
#include "TParticlePDG.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TFile.h"
#include "TRandom.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliMCParticle.h"
#include "AliMCEvent.h"

#include "AliLog.h"
#include "AliESDVertex.h"
#include "AliESDInputHandler.h"
#include "AliESDtrackCuts.h"
#include "AliMultiplicity.h"

#include "AliAnalysisTaskHardSoft.h"
#include "AliExternalTrackParam.h"
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliGenDPMjetEventHeader.h"

// Analysis Task for Hard and Soft event characteristics

// Author: E.Sicking

ClassImp(AliAnalysisTaskHardSoft)

//________________________________________________________________________
  AliAnalysisTaskHardSoft::AliAnalysisTaskHardSoft(const char *name) 
    : AliAnalysisTaskSE(name) 
    ,fUseMc(false)
    ,fUseNeutral(false)
    ,fRadiusCut(0.7)
    ,fTriggerPtCut(0.7)
    ,fAssociatePtCut(0.4)
    ,fMap(0x0)
    ,fMapLeading(0x0)
    ,fCuts(0)
    ,fFieldOn(kTRUE)
    ,fHists(0)
    
{
  for(Int_t i = 0;i< 2;i++){
    fPt[i]                  = 0;
    fEta[i]                 = 0;
    fPhi[i]                 = 0;
    fEtaPhi[i]              = 0;
    fEtaPhiLeading[i]       = 0;
    fNch[i]                 = 0;
    fPtLeading[i]           = 0;
    fPtLeadingNch[i]        = 0;
    fPtSumNch[i]            = 0;
    fPtAvNch[i]             = 0;
    fDPhiLeading[i]         = 0;
    fRadiusLeading[i]       = 0;
    fDPhiLeadingR[i]        = 0;
    fRadiusLeadingR[i]      = 0;
    fDPhiLeadingRS[i]       = 0;
    fRadiusLeadingRS[i]     = 0;
    fNchAssInR[i]           = 0;
    fTrigger[i]             = 0;

    for(Int_t j=0;j<100;j++){
      fDPhiLeadingNchBin[i][j]   = 0;
    }

    for(Int_t j=0;j<2;j++){
      fNchHardSoft[i][j]   = 0;
      fPtHardSoft[i][j]   = 0;
      fPtAvHardSoftNch[i][j]   = 0;
    }
  }
  DefineOutput(1,  TList::Class()); 
}


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

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  fHists = new TList();

 
  TString labels[2];
  labels[0]="MC";
  labels[1]="ESD";

  for(Int_t i=0;i<2;i++){
    fPt[i]                 = new TH1F(Form("fPt%s",labels[i].Data()),
				      Form("fPt%s",labels[i].Data()) ,  
				      500, 0., 50);
    fEta[i]                = new TH1F (Form("fEta%s",labels[i].Data()),
				       Form("fEta%s",labels[i].Data()) ,  
				       100, -1., 1);
    fPhi[i]                = new TH1F(Form("fPhi%s",labels[i].Data()),
				      Form("fPhi%s",labels[i].Data()) ,  
				      360, 0.,2*TMath::Pi());
    fEtaPhi[i]             = new TH2F(Form("fEtaPhi%s",labels[i].Data()),
				      Form("fEtaPhi%s",labels[i].Data()) ,  
				      100,-1.,1.,
				      360, 0.,2*TMath::Pi());
    fEtaPhiLeading[i]      = new TH2F(Form("fEtaPhiLeading%s",labels[i].Data()),
				      Form("fEtaPhiLeading%s",labels[i].Data()) ,  
				      100,-1.,1.,
				      360, 0.,2*TMath::Pi());
    fNch[i]                = new TH1F(Form("fNch%s",labels[i].Data()),
				      Form("fNch%s",labels[i].Data()) ,  
				      250, -0.5, 249.5);
    fPtLeading[i]          = new TH1F(Form("fPtLeading%s",labels[i].Data()),
				      Form("fPtLeading%s",labels[i].Data()) ,  
				      500, 0., 50);
    fPtLeadingNch[i]       = new TProfile(Form("fPtLeadingNch%s",labels[i].Data()),
					  Form("fPtLeadingNch%s",labels[i].Data()) ,  
					  250, -0.5, 249.5);
    fPtSumNch[i]           = new TProfile(Form("fPtSumNch%s",labels[i].Data()),
					  Form("fPtSumNch%s",labels[i].Data()) ,  
					  250, -0.5, 249.5);
    fPtAvNch[i]            = new TProfile(Form("fPtAvNch%s",labels[i].Data()),
					  Form("fPtAvNch%s",labels[i].Data()) ,  
					  250, -0.5, 249.5);
    fDPhiLeading[i]        = new TH1F(Form("fDPhiLeading%s",labels[i].Data()),
				      Form("fDPhiLeading%s",labels[i].Data()) ,  
				      180, 0., TMath::Pi());
    fRadiusLeading[i]      = new TH1F(Form("fRadiusLeading%s",labels[i].Data()),
				      Form("fRadiusLeading%s",labels[i].Data()) ,  
				      180, 0., 2*TMath::Pi());
    fDPhiLeadingR[i]       = new TH1F(Form("fDPhiLeadingR%s",labels[i].Data()),
				      Form("fDPhiLeadingR%s",labels[i].Data()) ,  
				      180, 0., TMath::Pi());
    fRadiusLeadingR[i]     = new TH1F(Form("fRadiusLeadingR%s",labels[i].Data()),
				      Form("fRadiusLeadingR%s",labels[i].Data()) ,  
				      180, 0., 2*TMath::Pi());
    fDPhiLeadingRS[i]      = new TH1F(Form("fDPhiLeadingRS%s",labels[i].Data()),
				      Form("fDPhiLeadingRS%s",labels[i].Data()) ,  
				      180, 0., TMath::Pi());
    fRadiusLeadingRS[i]    = new TH1F(Form("fRadiusLeadingRS%s",labels[i].Data()),
				      Form("fRadiusLeadingRS%s",labels[i].Data()) ,  
				      180, 0., 2*TMath::Pi());
    fNchAssInR[i]          = new TProfile(Form("fNchAssInR%s",labels[i].Data()),
					  Form("fNchAssInR%s",labels[i].Data()) ,  
					  250, -0.5, 249.5);
    fTrigger[i]            = new TH1F(Form("fTrigger%s",labels[i].Data()),
				      Form("fTrigger%s",labels[i].Data()) ,  
				      250, -0.5, 249.5);
    for(Int_t j=0;j<100;j++){
      fDPhiLeadingNchBin[i][j]     = new TH1F(Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j),
					      Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j) ,  
					      180, 0., TMath::Pi());
    }
    
    for(Int_t j=0;j<2;j++){
      fNchHardSoft[i][j]     = new TH1F(Form("fNchHardSoft%s%d",labels[i].Data(),j),
					Form("fNchHardSoft%s%d",labels[i].Data(),j) ,  
					250, -0.5, 249.5);
      fPtHardSoft[i][j]     = new TH1F(Form("fPtHardSoft%s%d",labels[i].Data(),j),
				       Form("fPtHardSoft%s%d",labels[i].Data(),j) ,  
				       500, 0., 50.0);
      fPtAvHardSoftNch[i][j]     = new TProfile(Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j),
						Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j) ,  
						250, -0.5, 249.5);
    }
    
  }
  
 
  fHists->SetOwner();
  for(Int_t i=0;i<2;i++){
    fHists->Add(fPt[i]);
    fHists->Add(fEta[i]);
    fHists->Add(fPhi[i]);
    fHists->Add(fEtaPhi[i]);
    fHists->Add(fEtaPhiLeading[i]);
    fHists->Add(fNch[i]);
    fHists->Add(fPtLeading[i]);
    fHists->Add(fPtLeadingNch[i]);
    fHists->Add(fPtSumNch[i]);
    fHists->Add(fPtAvNch[i]);
    fHists->Add(fDPhiLeading[i]);
    fHists->Add(fRadiusLeading[i]);
    fHists->Add(fDPhiLeadingR[i]);
    fHists->Add(fRadiusLeadingR[i]);
    fHists->Add(fDPhiLeadingRS[i]);
    fHists->Add(fRadiusLeadingRS[i]);
    fHists->Add(fNchAssInR[i]);
    fHists->Add(fTrigger[i]);

    for(Int_t j=0;j<100;j++){
      fHists->Add(fDPhiLeadingNchBin[i][j]);
    }

    for(Int_t j=0;j<2;j++){
      fHists->Add(fNchHardSoft[i][j]);
      fHists->Add(fPtHardSoft[i][j]);
      fHists->Add(fPtAvHardSoftNch[i][j]);
    }
  }

 
  TH1::AddDirectory(oldStatus);
}

//__________________________________________________________

void AliAnalysisTaskHardSoft::UserExec(Option_t *) 
{
  Int_t nentries[2]         = {0,0};  // tracks in event (before selection)
  Int_t nTracksAll[2]       = {0,0};  // accepted tracks in event
  Int_t nTracksAssociate[2] = {0,0};  // accepted tracks in event within R=fRadiusCut around leading track
  
  Float_t pt[2]    = {0.,0.}; // pt
  Float_t eta[2]   = {0.,0.}; // eta
  Float_t phi[2]   = {0.,0.}; // phi
  Float_t ptSum[2] = {0.,0.}; // pt sum
  Float_t ptAv[2]  = {0.,0.}; // average pt


  Float_t ptLeading[2]     = {0.,0.}; // pt leading
  Float_t etaLeading[2]    = {0.,0.}; // eta leading
  Float_t phiLeading[2]    = {0.,0.}; // phi leading

  Float_t ptOthers[2]     = {0.,0.}; // pt others // for second track loop around leading track
  Float_t etaOthers[2]    = {0.,0.}; // eta others
  Float_t phiOthers[2]    = {0.,0.}; // phi others

  Double_t etaLeadingRandom[2]   = {0.,0.}; // eta leading for random particle position (from fMapLeading)
  Double_t phiLeadingRandom[2]   = {0.,0.}; // phi leading  "
  Double_t etaOthersRandom[2]    = {0.,0.}; // eta others  for random particle position (from fMap)
  Double_t phiOthersRandom[2]    = {0.,0.}; // phi others   "

  Int_t fEventType[2]={1,1}; //hard=0, soft=1 -> set to hard if trigger and associate particle 
                             //are within R=0.7 

  //get event and vertex cut :(MC and ESD)
  //---------------------
  //MC
  //---------------------
  AliStack *stack = 0x0; // needed for MC studies
  Float_t vzMC=0.;       // MC vertex position in z
  if(fUseMc==true){
    stack = MCEvent()->Stack();
    AliGenEventHeader*  header = MCEvent()->GenEventHeader();
    TArrayF mcV;
    header->PrimaryVertex(mcV);
    vzMC = mcV[2];
  }
  //---------------------
  //ESD
  //---------------------
  AliVEvent* event = InputEvent();
  if (!event) {
    Printf("ERROR: Could not retrieve event");
    return;
  }
  if(Entry()==0){
    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
    if(esd)Printf("We are reading from ESD");
  }
  const AliVVertex* vertex = event->GetPrimaryVertex();
  Float_t vz = vertex->GetZ();
  //---------------------


  //Selection of particle(mcesd=0) or esdtracks(mcesd=1)
  for(Int_t mcesd=0;mcesd<2;mcesd++){
    if(mcesd==0){
      if(fUseMc==false) //MC part can be switched off for real data by function SetUserMc(kFALSE)
	continue;
    }

    // vertex cut and number of particles/tracks per event
    //---------------------------------------
    if(mcesd==0){//mc particles
      if (TMath::Abs(vzMC) > 10.) return;
      nentries[mcesd]=MCEvent()->GetNumberOfTracks();
    }
    else{// esd tracks
      if (TMath::Abs(vz) > 10.) return;
      nentries[mcesd]=event->GetNumberOfTracks();
    }//---------------------------------------
   

    // arrays for leading track determination (done with TMath::Sort of Array)
    Float_t * ptArray = new Float_t[nentries[mcesd]];

    //array of track properties
    Float_t * etaArray = new Float_t[nentries[mcesd]];
    Float_t * phiArray = new Float_t[nentries[mcesd]];

    Int_t *pindex  = new Int_t[nentries[mcesd]];
    for (Int_t i = 0; i < nentries[mcesd]; i++) {
      ptArray[i]=0.;
      etaArray[i]=0.;
      phiArray[i]=0.;
      pindex[i]=0;
    }
  
      
    //first track loop
    for (Int_t iTrack = 0; iTrack < nentries[mcesd]; iTrack++) {

      pt[mcesd]  = 0.;
      eta[mcesd] = 0.;
      phi[mcesd] = 0;
         
      //get properties of mc particle
      if(mcesd==0){//mc
	AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
	// Primaries only
	if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
	if(!fUseNeutral)if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
	//same cuts as on ESDtracks
	if(TMath::Abs(mcP->Eta())>0.9)continue;
	if(mcP->Pt()<0.2)continue;
	if(mcP->Pt()>200)continue;

	pt[mcesd]  = mcP->Pt();
	eta[mcesd] = mcP->Eta();
	phi[mcesd] = mcP->Phi();
      } 

      //get properties of esdtracks
      else{//esd
	AliVParticle *track = event->GetTrack(iTrack);
	if (!track) {
	  Printf("ERROR: Could not receive track %d", iTrack);
	  continue;
	}
	AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
	if(!esdtrack)continue;
	if (!fCuts->AcceptTrack(esdtrack)) continue;
 
	pt[mcesd]=esdtrack->Pt();
	eta[mcesd]=esdtrack->Eta();
	phi[mcesd]=esdtrack->Phi();
      }
    
      ptArray[nTracksAll[mcesd]]    = pt[mcesd];  // fill pt array
      etaArray[nTracksAll[mcesd]]   = eta[mcesd]; // fill eta array
      phiArray[nTracksAll[mcesd]++] = phi[mcesd]; // count tracks and fill phi array

      fPt[mcesd]     -> Fill(pt[mcesd]);
      fEta[mcesd]    -> Fill(eta[mcesd]);
      fPhi[mcesd]    -> Fill(phi[mcesd]);
      fEtaPhi[mcesd] -> Fill(eta[mcesd],phi[mcesd]);
      

    }//end first track loop

    fNch[mcesd]    -> Fill(nTracksAll[mcesd]);
    

    //find leading pt tracks
    if(nentries[mcesd]>0) TMath::Sort(nentries[mcesd], ptArray, pindex, kTRUE);  
    //for(Int_t i=0;i<nTracksAll[mcesd];i++){//show just the filled entries, skip empty ones.
    //     printf("%d:  pt = %f, number %i \n",mcesd, ptArray[pindex[i]],i);
    //}
    
    
    if(nTracksAll[mcesd]>0){
      fPtLeading[mcesd]->Fill(ptArray[pindex[0]]);  //first entry in array is highest
      fPtLeadingNch[mcesd]->Fill(nTracksAll[mcesd],ptArray[pindex[0]]);

      for(Int_t i=0;i<nTracksAll[mcesd];i++){      //calculate ptsum
	ptSum[mcesd]+=ptArray[pindex[i]];
      }
      ptAv[mcesd]=ptSum[mcesd]/nTracksAll[mcesd];  //calculate <pt>

      fPtSumNch[mcesd]->Fill(nTracksAll[mcesd],ptSum[mcesd]);
      fPtAvNch[mcesd]->Fill(nTracksAll[mcesd],ptAv[mcesd]);
    }
    

    if(nTracksAll[mcesd]>1){ // require at least two tracks (leading and prob. accosicates)
      
      //Leading track properties
      ptLeading[mcesd]  = ptArray[pindex[0]];
      etaLeading[mcesd] = etaArray[pindex[0]];
      phiLeading[mcesd] = phiArray[pindex[0]];

      fEtaPhiLeading[mcesd] -> Fill(etaLeading[mcesd],phiLeading[mcesd]);

      fMapLeading->GetRandom2(etaLeadingRandom[mcesd],phiLeadingRandom[mcesd]);
           
      //second track loop for event propoerties around leading tracks with pt>triggerPtCut
      //loop only over already accepted tracks except leading track 
      if(ptLeading[mcesd]>fTriggerPtCut){
	
	fTrigger[mcesd]->Fill(nTracksAll[mcesd]); // how often is there a trigger particle at a certain Nch bin
	
	for (Int_t iTrack = 1; iTrack < nTracksAll[mcesd]; iTrack++) { // starting at second highest track

	  ptOthers[mcesd]   = ptArray[pindex[iTrack]];
	  etaOthers[mcesd]  = etaArray[pindex[iTrack]];
	  phiOthers[mcesd]  = phiArray[pindex[iTrack]];

	  fMap->GetRandom2(etaOthersRandom[mcesd],phiOthersRandom[mcesd]);
	  
	  if(ptOthers[mcesd]>fAssociatePtCut){ // only tracks which fullfill associate pt cut

	    //1. real data

	    Float_t dPhi=TMath::Abs(phiOthers[mcesd]-phiLeading[mcesd]);
	    if(dPhi>TMath::Pi())      dPhi=2*TMath::Pi()-dPhi;
	    Float_t dEta=etaOthers[mcesd]-etaLeading[mcesd];
	    
	    Float_t radius=TMath::Sqrt(dPhi*dPhi+dEta*dEta);
	    fRadiusLeading[mcesd]->Fill(radius);
	    fDPhiLeading[mcesd]->Fill(dPhi);
	    if(nTracksAll[mcesd]<100)fDPhiLeadingNchBin[mcesd][nTracksAll[mcesd]]->Fill(dPhi);
	    
	    if(radius<fRadiusCut){
	      fEventType[mcesd]=0;
	      nTracksAssociate[mcesd]++;
	    }
	    
	    //2. random position
	    Float_t dPhiR=TMath::Abs(phiOthersRandom[mcesd]-phiLeadingRandom[mcesd]);
	    if(dPhiR>TMath::Pi())      dPhiR=dPhiR-2*TMath::Pi();
	    Float_t dEtaR=etaOthersRandom[mcesd]-etaLeadingRandom[mcesd];
	    
	    Float_t radiusR=TMath::Sqrt(dPhiR*dPhiR+dEtaR*dEtaR);
	    fRadiusLeadingR[mcesd]->Fill(radiusR);
	    fDPhiLeadingR[mcesd]->Fill(dPhiR);


	    //3. random position of leading particle
	    Float_t dPhiRS=TMath::Abs(phiOthers[mcesd]-phiLeadingRandom[mcesd]);
	    if(dPhiRS>TMath::Pi())      dPhiRS=dPhiRS-2*TMath::Pi();
	    Float_t dEtaRS=etaOthers[mcesd]-etaLeadingRandom[mcesd];
	    
	    Float_t radiusRS=TMath::Sqrt(dPhiRS*dPhiRS+dEtaRS*dEtaRS);
	    fRadiusLeadingRS[mcesd]->Fill(radiusRS);
	    fDPhiLeadingRS[mcesd]->Fill(dPhiRS);
	  }
	  

	}
	//fill histogram with number of tracks (pt>fAssociatePtCut) around leading track
	fNchAssInR[mcesd]->Fill(nTracksAll[mcesd],nTracksAssociate[mcesd]);
	
      }
    }
    
    fNchHardSoft[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd]);

    for(Int_t i=0;i<nTracksAll[mcesd];i++){
      fPtHardSoft[mcesd][fEventType[mcesd]]->Fill(ptArray[i]);
    }
    
    fPtAvHardSoftNch[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd],ptAv[mcesd]);

   
  }//double loop over mcP and ESD


  // Post output data.
  PostData(1, fHists);
}      


//________________________________________________________________________
void AliAnalysisTaskHardSoft::Terminate(Option_t *) 
{


}
Commit	Line	Data
7f092c14	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TH1F.h"
	4	#include "TH2F.h"
	5	#include "TH3F.h"
	6	#include "TCanvas.h"
	7	#include "TList.h"
	8	#include "TParticle.h"
	9	#include "TParticlePDG.h"
	10	#include "TProfile.h"
	11	#include "TNtuple.h"
	12	#include "TFile.h"
	13	#include "TRandom.h"
	14
	15	#include "AliAnalysisTask.h"
	16	#include "AliAnalysisManager.h"
	17
	18	#include "AliESDEvent.h"
	19	#include "AliStack.h"
	20	#include "AliMCParticle.h"
	21	#include "AliMCEvent.h"
	22
	23	#include "AliLog.h"
	24	#include "AliESDVertex.h"
	25	#include "AliESDInputHandler.h"
	26	#include "AliESDtrackCuts.h"
	27	#include "AliMultiplicity.h"
	28
	29	#include "AliAnalysisTaskHardSoft.h"
	30	#include "AliExternalTrackParam.h"
	31	#include "AliTrackReference.h"
	32	#include "AliHeader.h"
	33	#include "AliGenEventHeader.h"
	34	#include "AliGenDPMjetEventHeader.h"
	35
	36	// Analysis Task for Hard and Soft event characteristics
	37
	38	// Author: E.Sicking
	39
	40	ClassImp(AliAnalysisTaskHardSoft)
	41
	42	//________________________________________________________________________
	43	AliAnalysisTaskHardSoft::AliAnalysisTaskHardSoft(const char *name)
	44	: AliAnalysisTaskSE(name)
	45	,fUseMc(false)
	46	,fUseNeutral(false)
	47	,fRadiusCut(0.7)
	48	,fTriggerPtCut(0.7)
	49	,fAssociatePtCut(0.4)
	50	,fMap(0x0)
	51	,fMapLeading(0x0)
	52	,fCuts(0)
	53	,fFieldOn(kTRUE)
	54	,fHists(0)
	55
	56	{
	57	for(Int_t i = 0;i< 2;i++){
	58	fPt[i] = 0;
	59	fEta[i] = 0;
	60	fPhi[i] = 0;
	61	fEtaPhi[i] = 0;
	62	fEtaPhiLeading[i] = 0;
	63	fNch[i] = 0;
	64	fPtLeading[i] = 0;
65	fPtLeadingNch[i] = 0;
66	fPtSumNch[i] = 0;
67	fPtAvNch[i] = 0;
68	fDPhiLeading[i] = 0;
69	fRadiusLeading[i] = 0;
70	fDPhiLeadingR[i] = 0;
71	fRadiusLeadingR[i] = 0;
72	fDPhiLeadingRS[i] = 0;
73	fRadiusLeadingRS[i] = 0;
74	fNchAssInR[i] = 0;
75	fTrigger[i] = 0;
76
77	for(Int_t j=0;j<100;j++){
78	fDPhiLeadingNchBin[i][j] = 0;
79	}
a2b9ced6	80
	81	for(Int_t j=0;j<2;j++){
	82	fNchHardSoft[i][j] = 0;
	83	fPtHardSoft[i][j] = 0;
	84	fPtAvHardSoftNch[i][j] = 0;
	85	}
7f092c14	86	}
	87	DefineOutput(1, TList::Class());
	88	}
	89
	90
	91	//________________________________________________________________________
	92	void AliAnalysisTaskHardSoft::UserCreateOutputObjects()
	93	{
	94	// Create histograms
	95	// Called once
	96
	97	Bool_t oldStatus = TH1::AddDirectoryStatus();
	98	TH1::AddDirectory(kFALSE);
	99
	100	fHists = new TList();
	101
	102
	103
	104	TString labels[2];
	105	labels[0]="MC";
	106	labels[1]="ESD";
	107
	108	for(Int_t i=0;i<2;i++){
	109	fPt[i] = new TH1F(Form("fPt%s",labels[i].Data()),
a2b9ced6	110	Form("fPt%s",labels[i].Data()) ,
a2b9ced6	111	500, 0., 50);
7f092c14	112	fEta[i] = new TH1F (Form("fEta%s",labels[i].Data()),
	113	Form("fEta%s",labels[i].Data()) ,
	114	100, -1., 1);
	115	fPhi[i] = new TH1F(Form("fPhi%s",labels[i].Data()),
	116	Form("fPhi%s",labels[i].Data()) ,
	117	360, 0.,2*TMath::Pi());
	118	fEtaPhi[i] = new TH2F(Form("fEtaPhi%s",labels[i].Data()),
	119	Form("fEtaPhi%s",labels[i].Data()) ,
	120	100,-1.,1.,
	121	360, 0.,2*TMath::Pi());
	122	fEtaPhiLeading[i] = new TH2F(Form("fEtaPhiLeading%s",labels[i].Data()),
	123	Form("fEtaPhiLeading%s",labels[i].Data()) ,
	124	100,-1.,1.,
	125	360, 0.,2*TMath::Pi());
	126	fNch[i] = new TH1F(Form("fNch%s",labels[i].Data()),
	127	Form("fNch%s",labels[i].Data()) ,
	128	250, -0.5, 249.5);
	129	fPtLeading[i] = new TH1F(Form("fPtLeading%s",labels[i].Data()),
	130	Form("fPtLeading%s",labels[i].Data()) ,
	131	500, 0., 50);
	132	fPtLeadingNch[i] = new TProfile(Form("fPtLeadingNch%s",labels[i].Data()),
	133	Form("fPtLeadingNch%s",labels[i].Data()) ,
	134	250, -0.5, 249.5);
	135	fPtSumNch[i] = new TProfile(Form("fPtSumNch%s",labels[i].Data()),
	136	Form("fPtSumNch%s",labels[i].Data()) ,
	137	250, -0.5, 249.5);
	138	fPtAvNch[i] = new TProfile(Form("fPtAvNch%s",labels[i].Data()),
a2b9ced6	139	Form("fPtAvNch%s",labels[i].Data()) ,
a2b9ced6	140	250, -0.5, 249.5);
7f092c14	141	fDPhiLeading[i] = new TH1F(Form("fDPhiLeading%s",labels[i].Data()),
	142	Form("fDPhiLeading%s",labels[i].Data()) ,
	143	180, 0., TMath::Pi());
	144	fRadiusLeading[i] = new TH1F(Form("fRadiusLeading%s",labels[i].Data()),
	145	Form("fRadiusLeading%s",labels[i].Data()) ,
	146	180, 0., 2*TMath::Pi());
	147	fDPhiLeadingR[i] = new TH1F(Form("fDPhiLeadingR%s",labels[i].Data()),
	148	Form("fDPhiLeadingR%s",labels[i].Data()) ,
	149	180, 0., TMath::Pi());
	150	fRadiusLeadingR[i] = new TH1F(Form("fRadiusLeadingR%s",labels[i].Data()),
	151	Form("fRadiusLeadingR%s",labels[i].Data()) ,
	152	180, 0., 2*TMath::Pi());
	153	fDPhiLeadingRS[i] = new TH1F(Form("fDPhiLeadingRS%s",labels[i].Data()),
	154	Form("fDPhiLeadingRS%s",labels[i].Data()) ,
	155	180, 0., TMath::Pi());
	156	fRadiusLeadingRS[i] = new TH1F(Form("fRadiusLeadingRS%s",labels[i].Data()),
	157	Form("fRadiusLeadingRS%s",labels[i].Data()) ,
	158	180, 0., 2*TMath::Pi());
	159	fNchAssInR[i] = new TProfile(Form("fNchAssInR%s",labels[i].Data()),
	160	Form("fNchAssInR%s",labels[i].Data()) ,
	161	250, -0.5, 249.5);
	162	fTrigger[i] = new TH1F(Form("fTrigger%s",labels[i].Data()),
	163	Form("fTrigger%s",labels[i].Data()) ,
	164	250, -0.5, 249.5);
	165	for(Int_t j=0;j<100;j++){
	166	fDPhiLeadingNchBin[i][j] = new TH1F(Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j),
	167	Form("fDPhiLeadingNchBin%s%02d",labels[i].Data(),j) ,
	168	180, 0., TMath::Pi());
	169	}
	170
a2b9ced6	171	for(Int_t j=0;j<2;j++){
	172	fNchHardSoft[i][j] = new TH1F(Form("fNchHardSoft%s%d",labels[i].Data(),j),
	173	Form("fNchHardSoft%s%d",labels[i].Data(),j) ,
	174	250, -0.5, 249.5);
	175	fPtHardSoft[i][j] = new TH1F(Form("fPtHardSoft%s%d",labels[i].Data(),j),
	176	Form("fPtHardSoft%s%d",labels[i].Data(),j) ,
	177	500, 0., 50.0);
	178	fPtAvHardSoftNch[i][j] = new TProfile(Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j),
	179	Form("fPtAvHardSoftNch%s%d",labels[i].Data(),j) ,
	180	250, -0.5, 249.5);
	181	}
	182
7f092c14	183	}
	184
	185
	186	fHists->SetOwner();
	187	for(Int_t i=0;i<2;i++){
	188	fHists->Add(fPt[i]);
	189	fHists->Add(fEta[i]);
	190	fHists->Add(fPhi[i]);
	191	fHists->Add(fEtaPhi[i]);
	192	fHists->Add(fEtaPhiLeading[i]);
	193	fHists->Add(fNch[i]);
	194	fHists->Add(fPtLeading[i]);
	195	fHists->Add(fPtLeadingNch[i]);
	196	fHists->Add(fPtSumNch[i]);
	197	fHists->Add(fPtAvNch[i]);
	198	fHists->Add(fDPhiLeading[i]);
	199	fHists->Add(fRadiusLeading[i]);
	200	fHists->Add(fDPhiLeadingR[i]);
	201	fHists->Add(fRadiusLeadingR[i]);
	202	fHists->Add(fDPhiLeadingRS[i]);
	203	fHists->Add(fRadiusLeadingRS[i]);
	204	fHists->Add(fNchAssInR[i]);
	205	fHists->Add(fTrigger[i]);
	206
	207	for(Int_t j=0;j<100;j++){
	208	fHists->Add(fDPhiLeadingNchBin[i][j]);
	209	}
a2b9ced6	210
	211	for(Int_t j=0;j<2;j++){
	212	fHists->Add(fNchHardSoft[i][j]);
	213	fHists->Add(fPtHardSoft[i][j]);
	214	fHists->Add(fPtAvHardSoftNch[i][j]);
	215	}
7f092c14	216	}
	217
	218
	219	TH1::AddDirectory(oldStatus);
	220	}
	221
	222	//__________________________________________________________
	223
	224	void AliAnalysisTaskHardSoft::UserExec(Option_t *)
	225	{
	226	Int_t nentries[2] = {0,0}; // tracks in event (before selection)
	227	Int_t nTracksAll[2] = {0,0}; // accepted tracks in event
	228	Int_t nTracksAssociate[2] = {0,0}; // accepted tracks in event within R=fRadiusCut around leading track
	229
	230	Float_t pt[2] = {0.,0.}; // pt
	231	Float_t eta[2] = {0.,0.}; // eta
	232	Float_t phi[2] = {0.,0.}; // phi
	233	Float_t ptSum[2] = {0.,0.}; // pt sum
	234	Float_t ptAv[2] = {0.,0.}; // average pt
	235
	236
	237	Float_t ptLeading[2] = {0.,0.}; // pt leading
	238	Float_t etaLeading[2] = {0.,0.}; // eta leading
	239	Float_t phiLeading[2] = {0.,0.}; // phi leading
	240
	241	Float_t ptOthers[2] = {0.,0.}; // pt others // for second track loop around leading track
	242	Float_t etaOthers[2] = {0.,0.}; // eta others
	243	Float_t phiOthers[2] = {0.,0.}; // phi others
	244
	245	Double_t etaLeadingRandom[2] = {0.,0.}; // eta leading for random particle position (from fMapLeading)
	246	Double_t phiLeadingRandom[2] = {0.,0.}; // phi leading "
	247	Double_t etaOthersRandom[2] = {0.,0.}; // eta others for random particle position (from fMap)
	248	Double_t phiOthersRandom[2] = {0.,0.}; // phi others "
	249
a2b9ced6	250	Int_t fEventType[2]={1,1}; //hard=0, soft=1 -> set to hard if trigger and associate particle
a2b9ced6	251	//are within R=0.7
7f092c14	252
	253	//get event and vertex cut :(MC and ESD)
	254	//---------------------
	255	//MC
	256	//---------------------
a2b9ced6	257	AliStack *stack = 0x0; // needed for MC studies
a2b9ced6	258	Float_t vzMC=0.; // MC vertex position in z
7f092c14	259	if(fUseMc==true){
	260	stack = MCEvent()->Stack();
	261	AliGenEventHeader* header = MCEvent()->GenEventHeader();
	262	TArrayF mcV;
	263	header->PrimaryVertex(mcV);
	264	vzMC = mcV[2];
	265	}
	266	//---------------------
	267	//ESD
	268	//---------------------
	269	AliVEvent* event = InputEvent();
	270	if (!event) {
	271	Printf("ERROR: Could not retrieve event");
	272	return;
	273	}
	274	if(Entry()==0){
	275	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
	276	if(esd)Printf("We are reading from ESD");
	277	}
	278	const AliVVertex* vertex = event->GetPrimaryVertex();
	279	Float_t vz = vertex->GetZ();
	280	//---------------------
	281
	282
	283
	284
	285	//Selection of particle(mcesd=0) or esdtracks(mcesd=1)
	286	for(Int_t mcesd=0;mcesd<2;mcesd++){
	287	if(mcesd==0){
	288	if(fUseMc==false) //MC part can be switched off for real data by function SetUserMc(kFALSE)
	289	continue;
	290	}
	291
	292	// vertex cut and number of particles/tracks per event
	293	//---------------------------------------
	294	if(mcesd==0){//mc particles
	295	if (TMath::Abs(vzMC) > 10.) return;
	296	nentries[mcesd]=MCEvent()->GetNumberOfTracks();
	297	}
	298	else{// esd tracks
	299	if (TMath::Abs(vz) > 10.) return;
	300	nentries[mcesd]=event->GetNumberOfTracks();
	301	}//---------------------------------------
	302
	303
	304
	305
	306	// arrays for leading track determination (done with TMath::Sort of Array)
	307	Float_t * ptArray = new Float_t[nentries[mcesd]];
	308
	309	//array of track properties
	310	Float_t * etaArray = new Float_t[nentries[mcesd]];
	311	Float_t * phiArray = new Float_t[nentries[mcesd]];
	312
	313	Int_t *pindex = new Int_t[nentries[mcesd]];
	314	for (Int_t i = 0; i < nentries[mcesd]; i++) {
	315	ptArray[i]=0.;
	316	etaArray[i]=0.;
	317	phiArray[i]=0.;
	318	pindex[i]=0;
	319	}
	320
	321
	322
323
324	//first track loop
325	for (Int_t iTrack = 0; iTrack < nentries[mcesd]; iTrack++) {
326
327	pt[mcesd] = 0.;
328	eta[mcesd] = 0.;
329	phi[mcesd] = 0;
330
331	//get properties of mc particle
332	if(mcesd==0){//mc
333	AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
334	// Primaries only
335	if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
336	if(!fUseNeutral)if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
337	//same cuts as on ESDtracks
338	if(TMath::Abs(mcP->Eta())>0.9)continue;
339	if(mcP->Pt()<0.2)continue;
340	if(mcP->Pt()>200)continue;
341
342	pt[mcesd] = mcP->Pt();
343	eta[mcesd] = mcP->Eta();
344	phi[mcesd] = mcP->Phi();
345	}
346
347	//get properties of esdtracks
348	else{//esd
349	AliVParticle *track = event->GetTrack(iTrack);
350	if (!track) {
351	Printf("ERROR: Could not receive track %d", iTrack);
352	continue;
353	}
354	AliESDtrack esdtrack = dynamic_cast<AliESDtrack>(track);
355	if(!esdtrack)continue;
356	if (!fCuts->AcceptTrack(esdtrack)) continue;
357
358	pt[mcesd]=esdtrack->Pt();
359	eta[mcesd]=esdtrack->Eta();
360	phi[mcesd]=esdtrack->Phi();
361	}
362
363	ptArray[nTracksAll[mcesd]] = pt[mcesd]; // fill pt array
364	etaArray[nTracksAll[mcesd]] = eta[mcesd]; // fill eta array
365	phiArray[nTracksAll[mcesd]++] = phi[mcesd]; // count tracks and fill phi array
366
367	fPt[mcesd] -> Fill(pt[mcesd]);
368	fEta[mcesd] -> Fill(eta[mcesd]);
369	fPhi[mcesd] -> Fill(phi[mcesd]);
370	fEtaPhi[mcesd] -> Fill(eta[mcesd],phi[mcesd]);
371
372
373	}//end first track loop
374
375	fNch[mcesd] -> Fill(nTracksAll[mcesd]);
376
377
378	//find leading pt tracks
379	if(nentries[mcesd]>0) TMath::Sort(nentries[mcesd], ptArray, pindex, kTRUE);
380	//for(Int_t i=0;i<nTracksAll[mcesd];i++){//show just the filled entries, skip empty ones.
381	// printf("%d: pt = %f, number %i \n",mcesd, ptArray[pindex[i]],i);
382	//}
383
384
385	if(nTracksAll[mcesd]>0){
386	fPtLeading[mcesd]->Fill(ptArray[pindex[0]]); //first entry in array is highest
387	fPtLeadingNch[mcesd]->Fill(nTracksAll[mcesd],ptArray[pindex[0]]);
388
389	for(Int_t i=0;i<nTracksAll[mcesd];i++){ //calculate ptsum
390	ptSum[mcesd]+=ptArray[pindex[i]];
391	}
392	ptAv[mcesd]=ptSum[mcesd]/nTracksAll[mcesd]; //calculate <pt>
393
394	fPtSumNch[mcesd]->Fill(nTracksAll[mcesd],ptSum[mcesd]);
395	fPtAvNch[mcesd]->Fill(nTracksAll[mcesd],ptAv[mcesd]);
396	}
397
398
399	if(nTracksAll[mcesd]>1){ // require at least two tracks (leading and prob. accosicates)
400
401	//Leading track properties
402	ptLeading[mcesd] = ptArray[pindex[0]];
403	etaLeading[mcesd] = etaArray[pindex[0]];
404	phiLeading[mcesd] = phiArray[pindex[0]];
405
406	fEtaPhiLeading[mcesd] -> Fill(etaLeading[mcesd],phiLeading[mcesd]);
407
408	fMapLeading->GetRandom2(etaLeadingRandom[mcesd],phiLeadingRandom[mcesd]);
409
410	//second track loop for event propoerties around leading tracks with pt>triggerPtCut
411	//loop only over already accepted tracks except leading track
412	if(ptLeading[mcesd]>fTriggerPtCut){
413
414	fTrigger[mcesd]->Fill(nTracksAll[mcesd]); // how often is there a trigger particle at a certain Nch bin
415
416	for (Int_t iTrack = 1; iTrack < nTracksAll[mcesd]; iTrack++) { // starting at second highest track
417
418	ptOthers[mcesd] = ptArray[pindex[iTrack]];
419	etaOthers[mcesd] = etaArray[pindex[iTrack]];
420	phiOthers[mcesd] = phiArray[pindex[iTrack]];
421
422	fMap->GetRandom2(etaOthersRandom[mcesd],phiOthersRandom[mcesd]);
423
424	if(ptOthers[mcesd]>fAssociatePtCut){ // only tracks which fullfill associate pt cut
425
426	//1. real data
427
428	Float_t dPhi=TMath::Abs(phiOthers[mcesd]-phiLeading[mcesd]);
429	if(dPhi>TMath::Pi()) dPhi=2*TMath::Pi()-dPhi;
430	Float_t dEta=etaOthers[mcesd]-etaLeading[mcesd];
431
432	Float_t radius=TMath::Sqrt(dPhidPhi+dEtadEta);
433	fRadiusLeading[mcesd]->Fill(radius);
434	fDPhiLeading[mcesd]->Fill(dPhi);
435	if(nTracksAll[mcesd]<100)fDPhiLeadingNchBin[mcesd][nTracksAll[mcesd]]->Fill(dPhi);
436
437	if(radius<fRadiusCut){
a2b9ced6	438	fEventType[mcesd]=0;
7f092c14	439	nTracksAssociate[mcesd]++;
	440	}
	441
	442	//2. random position
	443	Float_t dPhiR=TMath::Abs(phiOthersRandom[mcesd]-phiLeadingRandom[mcesd]);
	444	if(dPhiR>TMath::Pi()) dPhiR=dPhiR-2*TMath::Pi();
	445	Float_t dEtaR=etaOthersRandom[mcesd]-etaLeadingRandom[mcesd];
	446
	447	Float_t radiusR=TMath::Sqrt(dPhiRdPhiR+dEtaRdEtaR);
	448	fRadiusLeadingR[mcesd]->Fill(radiusR);
	449	fDPhiLeadingR[mcesd]->Fill(dPhiR);
	450
	451
	452
	453	//3. random position of leading particle
	454	Float_t dPhiRS=TMath::Abs(phiOthers[mcesd]-phiLeadingRandom[mcesd]);
	455	if(dPhiRS>TMath::Pi()) dPhiRS=dPhiRS-2*TMath::Pi();
	456	Float_t dEtaRS=etaOthers[mcesd]-etaLeadingRandom[mcesd];
	457
	458	Float_t radiusRS=TMath::Sqrt(dPhiRSdPhiRS+dEtaRSdEtaRS);
	459	fRadiusLeadingRS[mcesd]->Fill(radiusRS);
	460	fDPhiLeadingRS[mcesd]->Fill(dPhiRS);
	461	}
	462
	463
	464
	465	}
	466	//fill histogram with number of tracks (pt>fAssociatePtCut) around leading track
	467	fNchAssInR[mcesd]->Fill(nTracksAll[mcesd],nTracksAssociate[mcesd]);
a2b9ced6	468
7f092c14	469	}
	470	}
	471
a2b9ced6	472	fNchHardSoft[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd]);
	473
	474	for(Int_t i=0;i<nTracksAll[mcesd];i++){
	475	fPtHardSoft[mcesd][fEventType[mcesd]]->Fill(ptArray[i]);
	476	}
	477
	478	fPtAvHardSoftNch[mcesd][fEventType[mcesd]]->Fill(nTracksAll[mcesd],ptAv[mcesd]);
	479
7f092c14	480
	481	}//double loop over mcP and ESD
	482
	483
	484
	485	// Post output data.
	486	PostData(1, fHists);
	487	}
	488
	489
	490
	491
	492
	493	//________________________________________________________________________
	494	void AliAnalysisTaskHardSoft::Terminate(Option_t *)
	495	{
	496
	497
	498	}
	499
	500
	501
	502
	503