[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskDndeta.cxx

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "AliFMDAnalysisTaskDndeta.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliLog.h"
#include "AliESDVertex.h"
#include "TMath.h"
#include "AliFMDAnaParameters.h"
//#include "AliFMDGeometry.h"
#include "AliGenEventHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliHeader.h"
//#include "TDatabasePDG.h"
//#include "TParticlePDG.h"
#include "AliFMDStripIndex.h"
#include "AliESDInputHandler.h"
ClassImp(AliFMDAnalysisTaskDndeta)


AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
: fDebug(0),
  fOutputList(0),
  fInputList(0),
  fVertexString(0x0),
  fNevents(),
  fNMCevents(),
  fStandalone(kTRUE),
  fLastTrackByStrip(0)
{
  // Default constructor
  DefineInput (0, TList::Class());
  DefineOutput(0, TList::Class());
}
//_____________________________________________________________________
AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
    AliAnalysisTask(name, "Density"),
    fDebug(0),
    fOutputList(0),
    fInputList(0),
    fVertexString(0x0),
    fNevents(),
    fNMCevents(),
    fStandalone(kTRUE),
    fLastTrackByStrip(0)
{
  fStandalone = SE;
  if(fStandalone) {
    DefineInput (0, TList::Class());
    DefineInput(1, TObjString::Class());
    DefineOutput(0, TList::Class());
    
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  if(!fOutputList)
    fOutputList = new TList();
  fOutputList->SetName("BackgroundCorrected");
  
  
  TH2F* hMult = 0;
  TH1F* hHits = 0;
  TH2F* hMultTrVtx = 0;
  TH1F* hPrimVertexBin = 0;
  
  
  TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
  TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
			    hBgTmp->GetNbinsX(),
			    hBgTmp->GetXaxis()->GetXmin(),
			    hBgTmp->GetXaxis()->GetXmax());
  hPrimary->Sumw2();
  fOutputList->Add(hPrimary);
  Int_t nVtxbins = pars->GetNvtxBins();
  TH2F* hBg = 0;
  for(Int_t i = 0; i< nVtxbins; i++) {
     
    for(Int_t det =1; det<=3;det++)
      {
	Int_t nRings = (det==1 ? 1 : 2);
	for(Int_t ring = 0;ring<nRings;ring++)
	  {
	    Char_t ringChar = (ring == 0 ? 'I' : 'O');
	    Int_t  nSec     = (ring == 0 ? 20 : 40);
	    
	    
	    hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	    hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax(),
			      nSec, 0, 2*TMath::Pi());
	    hMultTrVtx  = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
				   hBg->GetNbinsX(),
				   hBg->GetXaxis()->GetXmin(),
				   hBg->GetXaxis()->GetXmax(),
				   nSec, 0, 2*TMath::Pi());
	    hHits  = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax());
	    hHits->Sumw2();
	    fOutputList->Add(hHits);
	    
	    hMult->Sumw2();
	    fOutputList->Add(hMult);

	    hMultTrVtx->Sumw2();
	    fOutputList->Add(hMultTrVtx);
	    
	  }
      } 
  }
  
  for(Int_t i = 0; i< nVtxbins; i++) {
   
    hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
			      Form("primmult_vtxbin%d",i),
			      hBgTmp->GetNbinsX(),
			      hBgTmp->GetXaxis()->GetXmin(),
			      hBgTmp->GetXaxis()->GetXmax());
    hPrimVertexBin->Sumw2();
    fOutputList->Add(hPrimVertexBin);
    //SPD part
    TH2F* hSPDMult  = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i),
			 hBgTmp->GetNbinsX(),
			 hBgTmp->GetXaxis()->GetXmin(),
			 hBgTmp->GetXaxis()->GetXmax(),
			 20, 0, 2*TMath::Pi());
    hSPDMult->Sumw2();
    fOutputList->Add(hSPDMult);
    TH2F* hSPDMultTrVtx  = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i),
			 hBgTmp->GetNbinsX(),
			 hBgTmp->GetXaxis()->GetXmin(),
			 hBgTmp->GetXaxis()->GetXmax(),
			 20, 0, 2*TMath::Pi());
    hSPDMultTrVtx->Sumw2();
    fOutputList->Add(hSPDMultTrVtx);
  }
  
  fNevents.SetBins(nVtxbins,0,nVtxbins);
  fNevents.SetName("nEvents");
  fNMCevents.SetBins(nVtxbins,0,nVtxbins);
  fNMCevents.SetName("nMCEvents");
  fOutputList->Add(&fNevents);
  fOutputList->Add(&fNMCevents);
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
{
  if(fStandalone) {
    fInputList   = (TList*)GetInputData(0);
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  fVertexString = (TObjString*)fInputList->At(0);
  
  Int_t vtxbin   = fVertexString->GetString().Atoi();
  
  fNevents.Fill(vtxbin);
  
  for(UShort_t det=1;det<=3;det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t ringChar = (ir == 0 ? 'I' : 'O');
            
      TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      
      TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      
      hMultTotal->Add(hMultInput);
      hMultTotalTrVtx->Add(hMultInputTrVtx);
            
    }
  }
  
  //SPD code
  TH2F* hMultSPDTotal      = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin));
  TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin));
  TH2F* hMultSPDInput      = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
  TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
  hMultSPDTotal->Add(hMultSPDInput);
  hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx);
  
  
  if(pars->GetProcessPrimary())
    ProcessPrimary();
  
  if(fStandalone) {
    PostData(0, fOutputList); 
  }
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  Int_t nVtxbins = pars->GetNvtxBins();
  
  for(UShort_t det=1;det<=3;det++) {
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t ringChar = (ir == 0 ? 'I' : 'O');
      for(Int_t i =0; i<nVtxbins; i++) {
	
	TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
	//TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
	TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i));
	
	TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
	TH1D* hMultProjTrVtx   = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
	//fOutputList->Add(hMultTrVtx);
	fOutputList->Add(hMultProj);
	fOutputList->Add(hMultProjTrVtx);
      }
    }
  }
  
  for(Int_t i =0; i<nVtxbins; i++) {
    
    TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i));
    TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i));
    
    TH1D* hMultProj   = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY());
    TH1D* hMultProjTrVtx   = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY());
   
    fOutputList->Add(hMultProj);
    fOutputList->Add(hMultProjTrVtx);
  
  }
  
  std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" events"<<std::endl;
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
  
  AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  AliMCEvent* mcEvent = eventHandler->MCEvent();
  if(!mcEvent)
    return;
  
  fLastTrackByStrip.Reset(-1);
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  AliMCParticle* particle = 0;
  AliStack* stack = mcEvent->Stack();
  
  TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
  AliHeader* header            = mcEvent->Header();
  AliGenEventHeader* genHeader = header->GenEventHeader();
  
  TArrayF vertex;
  genHeader->PrimaryVertex(vertex);
  if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
    return;
  Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
  Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
  Int_t    vertexBin       = (Int_t)vertexBinDouble;
    
  Bool_t firstTrack = kTRUE;
  
  // we loop over the primaries only unless we need the hits (diagnostics running slowly)
  Int_t nTracks = stack->GetNprimary();
  if(pars->GetProcessHits())
    nTracks = stack->GetNtrack();
  
  for(Int_t i = 0 ;i<nTracks;i++) {
    particle = (AliMCParticle*) mcEvent->GetTrack(i);
    if(!particle)
      continue;
   
    if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
      hPrimary->Fill(particle->Eta());
      

      TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
      hPrimVtxBin->Fill(particle->Eta());
      if(firstTrack) {
	fNMCevents.Fill(vertexBin);
	firstTrack = kFALSE;
      }
    
    }
    if(pars->GetProcessHits()) {
           
      for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
	
	AliTrackReference* ref = particle->GetTrackReference(j);
	UShort_t det,sec,strip;
	Char_t   ring;
	if(ref->DetectorId() != AliTrackReference::kFMD)
	  continue;
	AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
	Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
	if(particle->Charge() != 0 && i != thisStripTrack ) {
	  //Double_t x,y,z;
	  
	  Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
	  TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
	  
	
	  hHits->Fill(eta);
	  
	  Float_t nstrips = (ring =='O' ? 256 : 512);
	  
	  fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
	
	  if(strip >0)
	    fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
	  if(strip < (nstrips - 1))
	    fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
	  
	}
      
	
      }
      
      
    }

      
  }
}
//_____________________________________________________________________
//
//
// EOF
Commit	Line	Data
7c3e5162	1
	2	#include <TROOT.h>
	3	#include <TSystem.h>
	4	#include <TInterpreter.h>
	5	#include <TChain.h>
	6	#include <TFile.h>
	7	#include <TList.h>
	8	#include <iostream>
bb8a464f	9	#include "TH1F.h"
7c3e5162	10	#include "TH2F.h"
	11	#include "AliFMDAnalysisTaskDndeta.h"
	12	#include "AliAnalysisManager.h"
	13	#include "AliESDFMD.h"
	14	#include "AliESDEvent.h"
	15	#include "AliAODEvent.h"
	16	#include "AliAODHandler.h"
	17	#include "AliMCEventHandler.h"
	18	#include "AliStack.h"
	19	#include "AliLog.h"
	20	#include "AliESDVertex.h"
	21	#include "TMath.h"
	22	#include "AliFMDAnaParameters.h"
78f6f750	23	//#include "AliFMDGeometry.h"
bb8a464f	24	#include "AliGenEventHeader.h"
5c7816c0	25	#include "AliGenPythiaEventHeader.h"
cc066cb9	26	#include "AliHeader.h"
78f6f750	27	//#include "TDatabasePDG.h"
78f6f750	28	//#include "TParticlePDG.h"
cc066cb9	29	#include "AliFMDStripIndex.h"
5c7816c0	30	#include "AliESDInputHandler.h"
7c3e5162	31	ClassImp(AliFMDAnalysisTaskDndeta)
	32
	33
	34	AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
	35	: fDebug(0),
	36	fOutputList(0),
	37	fInputList(0),
7c3e5162	38	fVertexString(0x0),
7c3e5162	39	fNevents(),
bb8a464f	40	fNMCevents(),
bb8a464f	41	fStandalone(kTRUE),
b64db9b1	42	fLastTrackByStrip(0)
7c3e5162	43	{
	44	// Default constructor
	45	DefineInput (0, TList::Class());
	46	DefineOutput(0, TList::Class());
	47	}
	48	//_____________________________________________________________________
	49	AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
	50	AliAnalysisTask(name, "Density"),
	51	fDebug(0),
	52	fOutputList(0),
	53	fInputList(0),
7c3e5162	54	fVertexString(0x0),
7c3e5162	55	fNevents(),
bb8a464f	56	fNMCevents(),
bb8a464f	57	fStandalone(kTRUE),
b64db9b1	58	fLastTrackByStrip(0)
7c3e5162	59	{
	60	fStandalone = SE;
	61	if(fStandalone) {
	62	DefineInput (0, TList::Class());
	63	DefineInput(1, TObjString::Class());
	64	DefineOutput(0, TList::Class());
	65
	66	}
	67	}
	68	//_____________________________________________________________________
	69	void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
	70	{
	71	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	72
7c3e5162	73	if(!fOutputList)
	74	fOutputList = new TList();
	75	fOutputList->SetName("BackgroundCorrected");
	76
	77
	78	TH2F* hMult = 0;
cc066cb9	79	TH1F* hHits = 0;
6c63da78	80	TH2F* hMultTrVtx = 0;
bb8a464f	81	TH1F* hPrimVertexBin = 0;
7c3e5162	82
bb8a464f	83
5ea1e0a9	84	TH2F* hBgTmp = pars->GetBackgroundCorrection(1, 'I', 0);
bb8a464f	85	TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
5ea1e0a9	86	hBgTmp->GetNbinsX(),
	87	hBgTmp->GetXaxis()->GetXmin(),
	88	hBgTmp->GetXaxis()->GetXmax());
bb8a464f	89	hPrimary->Sumw2();
bb8a464f	90	fOutputList->Add(hPrimary);
7c3e5162	91	Int_t nVtxbins = pars->GetNvtxBins();
b85ea106	92	TH2F* hBg = 0;
	93	for(Int_t i = 0; i< nVtxbins; i++) {
	94
	95	for(Int_t det =1; det<=3;det++)
	96	{
	97	Int_t nRings = (det==1 ? 1 : 2);
	98	for(Int_t ring = 0;ring<nRings;ring++)
	99	{
	100	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	101	Int_t nSec = (ring == 0 ? 20 : 40);
	102
	103
	104
	105	hBg = pars->GetBackgroundCorrection(det, ringChar, i);
7c3e5162	106	hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
	107	hBg->GetNbinsX(),
	108	hBg->GetXaxis()->GetXmin(),
	109	hBg->GetXaxis()->GetXmax(),
	110	nSec, 0, 2*TMath::Pi());
6c63da78	111	hMultTrVtx = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
	112	hBg->GetNbinsX(),
	113	hBg->GetXaxis()->GetXmin(),
	114	hBg->GetXaxis()->GetXmax(),
	115	nSec, 0, 2*TMath::Pi());
18d4b9ae	116	hHits = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
cc066cb9	117	hBg->GetNbinsX(),
	118	hBg->GetXaxis()->GetXmin(),
	119	hBg->GetXaxis()->GetXmax());
b85ea106	120	hHits->Sumw2();
b85ea106	121	fOutputList->Add(hHits);
cc066cb9	122
7c3e5162	123	hMult->Sumw2();
7c3e5162	124	fOutputList->Add(hMult);
6c63da78	125
	126	hMultTrVtx->Sumw2();
	127	fOutputList->Add(hMultTrVtx);
b85ea106	128
7c3e5162	129	}
b85ea106	130	}
b85ea106	131	}
7c3e5162	132
bb8a464f	133	for(Int_t i = 0; i< nVtxbins; i++) {
	134
	135	hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
	136	Form("primmult_vtxbin%d",i),
5ea1e0a9	137	hBgTmp->GetNbinsX(),
	138	hBgTmp->GetXaxis()->GetXmin(),
	139	hBgTmp->GetXaxis()->GetXmax());
bb8a464f	140	hPrimVertexBin->Sumw2();
bb8a464f	141	fOutputList->Add(hPrimVertexBin);
04f1ff3d	142	//SPD part
	143	TH2F* hSPDMult = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i),
	144	hBgTmp->GetNbinsX(),
	145	hBgTmp->GetXaxis()->GetXmin(),
	146	hBgTmp->GetXaxis()->GetXmax(),
	147	20, 0, 2*TMath::Pi());
	148	hSPDMult->Sumw2();
	149	fOutputList->Add(hSPDMult);
	150	TH2F* hSPDMultTrVtx = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i),
	151	hBgTmp->GetNbinsX(),
	152	hBgTmp->GetXaxis()->GetXmin(),
	153	hBgTmp->GetXaxis()->GetXmax(),
	154	20, 0, 2*TMath::Pi());
	155	hSPDMultTrVtx->Sumw2();
	156	fOutputList->Add(hSPDMultTrVtx);
bb8a464f	157	}
bb8a464f	158
7c3e5162	159	fNevents.SetBins(nVtxbins,0,nVtxbins);
7c3e5162	160	fNevents.SetName("nEvents");
bb8a464f	161	fNMCevents.SetBins(nVtxbins,0,nVtxbins);
bb8a464f	162	fNMCevents.SetName("nMCEvents");
7c3e5162	163	fOutputList->Add(&fNevents);
bb8a464f	164	fOutputList->Add(&fNMCevents);
7c3e5162	165
	166	}
	167	//_____________________________________________________________________
	168	void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t /option*/)
	169	{
	170	if(fStandalone) {
	171	fInputList = (TList*)GetInputData(0);
7c3e5162	172	}
	173	}
	174	//_____________________________________________________________________
	175	void AliFMDAnalysisTaskDndeta::Exec(Option_t /option*/)
	176	{
b64db9b1	177	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
9f55be54	178	fVertexString = (TObjString*)fInputList->At(0);
9f55be54	179
7c3e5162	180	Int_t vtxbin = fVertexString->GetString().Atoi();
9f55be54	181
7c3e5162	182	fNevents.Fill(vtxbin);
9f55be54	183
7c3e5162	184	for(UShort_t det=1;det<=3;det++) {
7c3e5162	185	Int_t nRings = (det==1 ? 1 : 2);
	186	for (UShort_t ir = 0; ir < nRings; ir++) {
	187	Char_t ringChar = (ir == 0 ? 'I' : 'O');
9f55be54	188
b64db9b1	189	TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
6c63da78	190	TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
7c3e5162	191
7c3e5162	192	TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
6c63da78	193	TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
1282ce49	194
7c3e5162	195	hMultTotal->Add(hMultInput);
6c63da78	196	hMultTotalTrVtx->Add(hMultInputTrVtx);
9f55be54	197
7c3e5162	198	}
7c3e5162	199	}
bb8a464f	200
04f1ff3d	201	//SPD code
	202	TH2F* hMultSPDTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin));
	203	TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin));
	204	TH2F* hMultSPDInput = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
	205	TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
	206	hMultSPDTotal->Add(hMultSPDInput);
	207	hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx);
	208
	209
b64db9b1	210	if(pars->GetProcessPrimary())
bb8a464f	211	ProcessPrimary();
bb8a464f	212
7c3e5162	213	if(fStandalone) {
	214	PostData(0, fOutputList);
	215	}
	216
	217	}
	218	//_____________________________________________________________________
	219	void AliFMDAnalysisTaskDndeta::Terminate(Option_t /option*/) {
	220
7c3e5162	221	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	222
	223	Int_t nVtxbins = pars->GetNvtxBins();
	224
	225	for(UShort_t det=1;det<=3;det++) {
7c3e5162	226	Int_t nRings = (det==1 ? 1 : 2);
7c3e5162	227	for (UShort_t ir = 0; ir < nRings; ir++) {
1282ce49	228	Char_t ringChar = (ir == 0 ? 'I' : 'O');
7c3e5162	229	for(Int_t i =0; i<nVtxbins; i++) {
b64db9b1	230
6c63da78	231	TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
	232	//TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
	233	TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i));
25f47050	234
1282ce49	235	TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
f1eff0e8	236	TH1D* hMultProjTrVtx = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
04f1ff3d	237	//fOutputList->Add(hMultTrVtx);
1282ce49	238	fOutputList->Add(hMultProj);
b64db9b1	239	fOutputList->Add(hMultProjTrVtx);
7c3e5162	240	}
	241	}
	242	}
04f1ff3d	243
	244	for(Int_t i =0; i<nVtxbins; i++) {
	245
	246	TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i));
	247	TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i));
	248
	249	TH1D* hMultProj = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY());
	250	TH1D* hMultProjTrVtx = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY());
	251
	252	fOutputList->Add(hMultProj);
	253	fOutputList->Add(hMultProjTrVtx);
	254
	255	}
	256
5c7816c0	257	std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" events"<<std::endl;
bb8a464f	258	}
	259	//_____________________________________________________________________
	260	void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
cc066cb9	261
b64db9b1	262	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	263	AliMCEvent* mcEvent = eventHandler->MCEvent();
	264	if(!mcEvent)
	265	return;
	266
cc066cb9	267	fLastTrackByStrip.Reset(-1);
cc066cb9	268
bb8a464f	269	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	270
	271	AliMCParticle* particle = 0;
b64db9b1	272	AliStack* stack = mcEvent->Stack();
bb8a464f	273
bb8a464f	274	TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
b64db9b1	275	AliHeader* header = mcEvent->Header();
cc066cb9	276	AliGenEventHeader* genHeader = header->GenEventHeader();
	277
	278	TArrayF vertex;
	279	genHeader->PrimaryVertex(vertex);
	280	if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
	281	return;
	282	Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
	283	Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
	284	Int_t vertexBin = (Int_t)vertexBinDouble;
0a2f2742	285
bb8a464f	286	Bool_t firstTrack = kTRUE;
bb0a45c3	287
	288	// we loop over the primaries only unless we need the hits (diagnostics running slowly)
	289	Int_t nTracks = stack->GetNprimary();
b64db9b1	290	if(pars->GetProcessHits())
bb0a45c3	291	nTracks = stack->GetNtrack();
bb0a45c3	292
bb8a464f	293	for(Int_t i = 0 ;i<nTracks;i++) {
7aad0c47	294	particle = (AliMCParticle*) mcEvent->GetTrack(i);
bb8a464f	295	if(!particle)
bb8a464f	296	continue;
bb0a45c3	297
bb8a464f	298	if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
bb8a464f	299	hPrimary->Fill(particle->Eta());
cc066cb9	300
6289b3e8	301
bb8a464f	302	TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
	303	hPrimVtxBin->Fill(particle->Eta());
	304	if(firstTrack) {
	305	fNMCevents.Fill(vertexBin);
	306	firstTrack = kFALSE;
	307	}
b3546e91	308
bb0a45c3	309	}
f1eff0e8	310	if(pars->GetProcessHits()) {
b85ea106	311
bb0a45c3	312	for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
cc066cb9	313
55cadbf9	314	AliTrackReference* ref = particle->GetTrackReference(j);
	315	UShort_t det,sec,strip;
	316	Char_t ring;
	317	if(ref->DetectorId() != AliTrackReference::kFMD)
	318	continue;
	319	AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
	320	Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
	321	if(particle->Charge() != 0 && i != thisStripTrack ) {
	322	//Double_t x,y,z;
	323
	324	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1TMath::Log(TMath::Tan(0.5theta));
18d4b9ae	325	TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
b85ea106	326
b85ea106	327
55cadbf9	328	hHits->Fill(eta);
b85ea106	329
55cadbf9	330	Float_t nstrips = (ring =='O' ? 256 : 512);
	331
	332	fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
cc066cb9	333
55cadbf9	334	if(strip >0)
	335	fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
	336	if(strip < (nstrips - 1))
	337	fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
	338
	339	}
b85ea106	340
b85ea106	341
b3546e91	342	}
b85ea106	343
b85ea106	344
bb0a45c3	345	}
b85ea106	346
	347
	348
bb0a45c3	349	}
7c3e5162	350	}
	351	//_____________________________________________________________________
	352	//
	353	//
	354	// EOF