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

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/

/* $Id:$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  Here some comments on what this task does                                //
 
#include <TChain.h>
#include <TF1.h>
#include <TFile.h>
#include <TInterpreter.h>
#include <TList.h>
#include <TMath.h>
#include <TROOT.h>
#include <TSystem.h>
#include <iostream>

#include "AliFMDAnalysisTaskCollector.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliMultiplicity.h"
#include "AliESDVertex.h"
#include "AliFMDAnaParameters.h"
//#include "AliFMDGeometry.h"


ClassImp(AliFMDAnalysisTaskCollector)

//____________________________________________________________________
Double_t  AliFMDAnalysisTaskCollector::TripleLandau(Double_t *x, Double_t *par) {
  
  Double_t energy        = x[0];
  Double_t constant = par[0];
  Double_t mpv      = par[1];
  Double_t sigma    = par[2];
  Double_t alpha    = par[3];
  Double_t beta     = par[4];
 
  Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
			 alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
			 beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
  
  return f;
}
//____________________________________________________________________

AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
: fDebug(0),
  fOutputList(0),
  fArray(0),
  fZvtxDist(0),
  fEvents(0)
{
  // Default constructor
  
 
}
//____________________________________________________________________
AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
    AliAnalysisTaskSE(name),
    fDebug(0),
    fOutputList(0),
    fArray(0),
    fZvtxDist(0),
    fEvents(0)
{
  // Default constructor
  
  DefineOutput(1, TList::Class());
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
{
  // Create the output container
  printf("AnalysisTaskFMD::CreateOutPutData() \n");
  
  fOutputList = new TList();//(TList*)GetOutputData(0);
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fArray     = new TObjArray();
  fArray->SetName("FMD");
  fArray->SetOwner();
  TH1F* hEdist = 0;
  TH1F* hEmptyEdist = 0;
  TH1F* hRingEdist = 0;
  for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
    TObjArray* etaArray = new TObjArray();
    fArray->AddAtAndExpand(etaArray,nEta);
    for(Int_t det =1; det<=3;det++)
      {
	TObjArray* detArray = new TObjArray();
	detArray->SetName(Form("FMD%d",det));
	etaArray->AddAtAndExpand(detArray,det);
	Int_t nRings = (det==1 ? 1 : 2);
	for(Int_t ring = 0;ring<nRings;ring++)
	  {
	    Char_t ringChar = (ring == 0 ? 'I' : 'O');
	    hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
	    hEdist->SetXTitle("#Delta E / E_{MIP}");
	    fOutputList->Add(hEdist);
	    detArray->AddAtAndExpand(hEdist,ring);
	    
	  } 
      }
    
  }
  
  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');
	  hRingEdist = new TH1F(Form("ringFMD%d%c",det,ringChar),Form("ringFMD%d%c",det,ringChar),200,0,6);
	  hRingEdist->SetXTitle("#Delta E / E_{MIP}");
	  fOutputList->Add(hRingEdist);
	  hEmptyEdist = new TH1F(Form("emptyFMD%d%c",det,ringChar),Form("emptyFMD%d%c",det,ringChar),200,0,6);
	  hEmptyEdist->SetXTitle("#Delta E / E_{MIP}");
	  fOutputList->Add(hEmptyEdist);
	  
	}
    }
  fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
  fZvtxDist->SetXTitle("z vertex");
  fOutputList->Add(fZvtxDist);
}

//____________________________________________________________________
void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
{
  
  AliESDEvent* esd = (AliESDEvent*)InputEvent();
  AliESD* old = esd->GetAliESDOld();
  if (old) {
    esd->CopyFromOldESD();
  }
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  TString triggers = esd->GetFiredTriggerClasses();
  //if(!triggers.IsNull()) return;
  //Bool_t trigger = pars->IsEventTriggered(esd);
  
  Bool_t physics = pars->IsEventTriggered(esd);
  Bool_t empty   = pars->IsEventTriggered(esd,AliFMDAnaParameters::kEMPTY);
  //std::cout<<physics<<"   "<<empty<<std::endl;
  //if(!trigger)
  //  physics = kFALSE;
  Double_t vertex[3];
  
  Bool_t vtxStatus =   pars->GetVertex(esd,vertex);
  if(!vtxStatus)
    physics = kFALSE;
  
  if(physics) {
    fZvtxDist->Fill(vertex[2]);
    if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
      physics = kFALSE;
  }
  
  AliESDFMD* fmd = esd->GetFMDData();
  if (!fmd) return;
  if(physics)
    fEvents++;
  else if(empty)
    fEmptyEvents++;
  
  if(!physics && !empty)
    return;
  TH1F* Edist = 0;
  TH1F* emptyEdist = 0;
  TH1F* ringEdist = 0;
  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   ring = (ir == 0 ? 'I' : 'O');
      emptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ring));
      ringEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ring));
      UShort_t nsec = (ir == 0 ? 20  : 40);
      UShort_t nstr = (ir == 0 ? 512 : 256);
      
      for(UShort_t sec =0; sec < nsec;  sec++)  {
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  
	  
	  Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	  if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
	  
	  Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
	  
	  Int_t nEta  = pars->GetEtaBin(eta);
	  //	  std::cout<<det<<"  "<<ring<<"   "<<sec<<"    "<<strip<<"   "<<vertex[2]<<"   "<<eta<<"   "<<nEta<<std::endl;
	  if(physics) {
	    Edist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ring,nEta));	  
	    Edist->Fill(mult);
	    ringEdist->Fill(mult);
	  }
	  else if(empty) {
	    emptyEdist->Fill(mult);
	    
	  }
	  else {
	    AliWarning("Something is wrong - wrong trigger");
	    continue;
	  }
	 
	  
	}
      }
    }
    
    PostData(1, fOutputList); 
    
  }
}
//____________________________________________________________________
  
void AliFMDAnalysisTaskCollector::Terminate(Option_t */*option*/) {
  std::cout<<"Analysed "<<fEvents<<" events and "<<fEmptyEvents<<" empty"<<std::endl;
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 
    
    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');
	TH1F* hRingEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ringChar));
	TH1F* hEmptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ringChar));
	if(fEmptyEvents)
	  hEmptyEdist->Scale(1./(Float_t)fEmptyEvents);
	
	if(fEvents)
	  hRingEdist->Scale(1./(Float_t)fEvents);
	for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
	  TH1F* hEdist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
	  if(fEvents)
	    hEdist->Scale(1./(Float_t)fEvents);
	  
	
      }
    }
  }
  
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {

  //speciesOption:
  //0: p+p Landau fit
  //1: Pb+Pb triple landau convolution fit
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
  
  TFile fin(filename,"UPDATE");
  
  TList* list = (TList*)fin.Get("energyDist");
  
  AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
  
  EnergyDist->SetNetaBins(pars->GetNetaBins());
  EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
    
  TF1* fitFunc = 0;
  
  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');
      TH1F* hEmptyEdist = (TH1F*)list->FindObject(Form("emptyFMD%d%c",det,ringChar));
      TH1F* hRingEdist = (TH1F*)list->FindObject(Form("ringFMD%d%c",det,ringChar));
      fitFunc = FitEnergyDistribution(hEmptyEdist, speciesOption) ;
      if(fitFunc)
	fitFunc->Write(Form("emptyFMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
      fitFunc = FitEnergyDistribution(hRingEdist, speciesOption) ;
      if(fitFunc)
	fitFunc->Write(Form("FMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
      
      
      EnergyDist->SetEmptyEnergyDistribution(det,ringChar,hEmptyEdist);
      EnergyDist->SetRingEnergyDistribution(det,ringChar,hRingEdist);
      for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
	TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
	
	fitFunc = FitEnergyDistribution(hEdist, speciesOption) ;
	if(fitFunc)
	  fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
	EnergyDist->SetEnergyDistribution(det,ringChar,nEta,hEdist);
	
      }
      
    }
    
  }
  
  fin.Write();
  fin.Close();
  
  if(store) {
    TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
    EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
    fcalib.Close();
  }


}
//____________________________________________________________________
TF1* AliFMDAnalysisTaskCollector::FitEnergyDistribution(TH1F* hEnergy, Int_t speciesOption) {
  
  TF1* fitFunc = 0;
  if(hEnergy->GetEntries() != 0) {
	  
    hEnergy->GetXaxis()->SetRangeUser(0.2,hEnergy->GetXaxis()->GetXmax());
    
    if(speciesOption == 0)
      fitFunc =  new TF1("FMDfitFunc","landau",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
    if(speciesOption == 1) {
      fitFunc = new TF1("FMDfitFunc",TripleLandau,hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,5,5);
      fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
      fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
      fitFunc->SetParLimits(1,0.6,1.2);
      fitFunc->SetParLimits(3,0,1);
      fitFunc->SetParLimits(4,0,1);
      
    }
    hEnergy->Fit(fitFunc,"","",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
    
  }
  return fitFunc;
  
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
96ee66b9	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	/* $Id:$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Here some comments on what this task does //
c48a797f	21
c48a797f	22	#include <TChain.h>
96ee66b9	23	#include <TF1.h>
c48a797f	24	#include <TFile.h>
96ee66b9	25	#include <TInterpreter.h>
c48a797f	26	#include <TList.h>
9f55be54	27	#include <TMath.h>
96ee66b9	28	#include <TROOT.h>
96ee66b9	29	#include <TSystem.h>
c48a797f	30	#include <iostream>
	31
	32	#include "AliFMDAnalysisTaskCollector.h"
	33	#include "AliAnalysisManager.h"
	34	#include "AliESDFMD.h"
	35	#include "AliESDEvent.h"
	36	#include "AliAODEvent.h"
	37	#include "AliAODHandler.h"
	38	#include "AliMCEventHandler.h"
	39	#include "AliStack.h"
884cadc9	40	#include "AliMultiplicity.h"
c48a797f	41	#include "AliESDVertex.h"
c48a797f	42	#include "AliFMDAnaParameters.h"
78f6f750	43	//#include "AliFMDGeometry.h"
9f55be54	44
9f55be54	45
c48a797f	46	ClassImp(AliFMDAnalysisTaskCollector)
c48a797f	47
9f55be54	48	//____________________________________________________________________
	49	Double_t AliFMDAnalysisTaskCollector::TripleLandau(Double_t x, Double_t par) {
	50
	51	Double_t energy = x[0];
	52	Double_t constant = par[0];
	53	Double_t mpv = par[1];
	54	Double_t sigma = par[2];
	55	Double_t alpha = par[3];
	56	Double_t beta = par[4];
	57
	58	Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
	59	alphaTMath::Landau(energy,2mpv+2sigmaTMath::Log(2),2*sigma,kTRUE)+
	60	betaTMath::Landau(energy,3mpv+3sigmaTMath::Log(3),3*sigma,kTRUE) );
	61
	62	return f;
	63	}
	64	//____________________________________________________________________
c48a797f	65
	66	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
	67	: fDebug(0),
c48a797f	68	fOutputList(0),
c48a797f	69	fArray(0),
884cadc9	70	fZvtxDist(0),
aa303f0c	71	fEvents(0)
c48a797f	72	{
c48a797f	73	// Default constructor
9f55be54	74
9f55be54	75
c48a797f	76	}
	77	//____________________________________________________________________
	78	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
9f55be54	79	AliAnalysisTaskSE(name),
c48a797f	80	fDebug(0),
c48a797f	81	fOutputList(0),
c48a797f	82	fArray(0),
884cadc9	83	fZvtxDist(0),
aa303f0c	84	fEvents(0)
c48a797f	85	{
c48a797f	86	// Default constructor
9f55be54	87
9f55be54	88	DefineOutput(1, TList::Class());
c48a797f	89	}
c48a797f	90	//____________________________________________________________________
9f55be54	91	void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
c48a797f	92	{
	93	// Create the output container
	94	printf("AnalysisTaskFMD::CreateOutPutData() \n");
	95
	96	fOutputList = new TList();//(TList*)GetOutputData(0);
6289b3e8	97	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c48a797f	98
	99	fArray = new TObjArray();
	100	fArray->SetName("FMD");
	101	fArray->SetOwner();
	102	TH1F* hEdist = 0;
884cadc9	103	TH1F* hEmptyEdist = 0;
	104	TH1F* hRingEdist = 0;
	105	for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
6289b3e8	106	TObjArray* etaArray = new TObjArray();
	107	fArray->AddAtAndExpand(etaArray,nEta);
	108	for(Int_t det =1; det<=3;det++)
	109	{
	110	TObjArray* detArray = new TObjArray();
	111	detArray->SetName(Form("FMD%d",det));
	112	etaArray->AddAtAndExpand(detArray,det);
	113	Int_t nRings = (det==1 ? 1 : 2);
	114	for(Int_t ring = 0;ring<nRings;ring++)
	115	{
	116	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	117	hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
	118	hEdist->SetXTitle("#Delta E / E_{MIP}");
	119	fOutputList->Add(hEdist);
	120	detArray->AddAtAndExpand(hEdist,ring);
884cadc9	121
6289b3e8	122	}
	123	}
	124
	125	}
c48a797f	126
884cadc9	127	for(Int_t det =1; det<=3;det++)
	128	{
	129	Int_t nRings = (det==1 ? 1 : 2);
	130	for(Int_t ring = 0;ring<nRings;ring++)
	131	{
	132	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	133	hRingEdist = new TH1F(Form("ringFMD%d%c",det,ringChar),Form("ringFMD%d%c",det,ringChar),200,0,6);
	134	hRingEdist->SetXTitle("#Delta E / E_{MIP}");
	135	fOutputList->Add(hRingEdist);
	136	hEmptyEdist = new TH1F(Form("emptyFMD%d%c",det,ringChar),Form("emptyFMD%d%c",det,ringChar),200,0,6);
	137	hEmptyEdist->SetXTitle("#Delta E / E_{MIP}");
	138	fOutputList->Add(hEmptyEdist);
	139
	140	}
	141	}
c48a797f	142	fZvtxDist = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
c48a797f	143	fZvtxDist->SetXTitle("z vertex");
c48a797f	144	fOutputList->Add(fZvtxDist);
c48a797f	145	}
c48a797f	146
c48a797f	147	//____________________________________________________________________
9f55be54	148	void AliFMDAnalysisTaskCollector::UserExec(Option_t /option*/)
c48a797f	149	{
c48a797f	150
9f55be54	151	AliESDEvent* esd = (AliESDEvent*)InputEvent();
9f55be54	152	AliESD* old = esd->GetAliESDOld();
c48a797f	153	if (old) {
9f55be54	154	esd->CopyFromOldESD();
c48a797f	155	}
bb8a464f	156	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
884cadc9	157	TString triggers = esd->GetFiredTriggerClasses();
	158	//if(!triggers.IsNull()) return;
	159	//Bool_t trigger = pars->IsEventTriggered(esd);
c48a797f	160
884cadc9	161	Bool_t physics = pars->IsEventTriggered(esd);
aa303f0c	162	Bool_t empty = pars->IsEventTriggered(esd,AliFMDAnaParameters::kEMPTY);
884cadc9	163	//std::cout<<physics<<" "<<empty<<std::endl;
	164	//if(!trigger)
	165	// physics = kFALSE;
c48a797f	166	Double_t vertex[3];
bb8a464f	167
884cadc9	168	Bool_t vtxStatus = pars->GetVertex(esd,vertex);
	169	if(!vtxStatus)
	170	physics = kFALSE;
bb8a464f	171
884cadc9	172	if(physics) {
	173	fZvtxDist->Fill(vertex[2]);
	174	if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
	175	physics = kFALSE;
	176	}
bb8a464f	177
9f55be54	178	AliESDFMD* fmd = esd->GetFMDData();
c48a797f	179	if (!fmd) return;
884cadc9	180	if(physics)
884cadc9	181	fEvents++;
aa303f0c	182	else if(empty)
aa303f0c	183	fEmptyEvents++;
c48a797f	184
aa303f0c	185	if(!physics && !empty)
884cadc9	186	return;
	187	TH1F* Edist = 0;
	188	TH1F* emptyEdist = 0;
	189	TH1F* ringEdist = 0;
c48a797f	190	for(UShort_t det=1;det<=3;det++) {
884cadc9	191
c48a797f	192	Int_t nRings = (det==1 ? 1 : 2);
c48a797f	193	for (UShort_t ir = 0; ir < nRings; ir++) {
c48a797f	194	Char_t ring = (ir == 0 ? 'I' : 'O');
884cadc9	195	emptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ring));
884cadc9	196	ringEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ring));
c48a797f	197	UShort_t nsec = (ir == 0 ? 20 : 40);
c48a797f	198	UShort_t nstr = (ir == 0 ? 512 : 256);
6289b3e8	199
c48a797f	200	for(UShort_t sec =0; sec < nsec; sec++) {
c48a797f	201	for(UShort_t strip = 0; strip < nstr; strip++) {
cc066cb9	202
	203
	204	Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	205	if(mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
cc066cb9	206
78f6f750	207	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
9f55be54	208
884cadc9	209	Int_t nEta = pars->GetEtaBin(eta);
	210	// std::cout<<det<<" "<<ring<<" "<<sec<<" "<<strip<<" "<<vertex[2]<<" "<<eta<<" "<<nEta<<std::endl;
	211	if(physics) {
	212	Edist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ring,nEta));
	213	Edist->Fill(mult);
	214	ringEdist->Fill(mult);
	215	}
aa303f0c	216	else if(empty) {
884cadc9	217	emptyEdist->Fill(mult);
aa303f0c	218
884cadc9	219	}
aa303f0c	220	else {
	221	AliWarning("Something is wrong - wrong trigger");
	222	continue;
884cadc9	223	}
884cadc9	224
c48a797f	225
	226	}
	227	}
	228	}
884cadc9	229
	230	PostData(1, fOutputList);
	231
c48a797f	232	}
884cadc9	233	}
884cadc9	234	//____________________________________________________________________
c48a797f	235
884cadc9	236	void AliFMDAnalysisTaskCollector::Terminate(Option_t /option*/) {
aa303f0c	237	std::cout<<"Analysed "<<fEvents<<" events and "<<fEmptyEvents<<" empty"<<std::endl;
884cadc9	238	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	239
	240
	241	for(Int_t det = 1; det<=3; det++) {
	242	Int_t nRings = (det == 1 ? 1 : 2);
	243	for(Int_t ring = 0;ring<nRings; ring++) {
	244
	245	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	246	TH1F* hRingEdist = (TH1F*)fOutputList->FindObject(Form("ringFMD%d%c",det,ringChar));
	247	TH1F* hEmptyEdist = (TH1F*)fOutputList->FindObject(Form("emptyFMD%d%c",det,ringChar));
aa303f0c	248	if(fEmptyEvents)
aa303f0c	249	hEmptyEdist->Scale(1./(Float_t)fEmptyEvents);
884cadc9	250
	251	if(fEvents)
	252	hRingEdist->Scale(1./(Float_t)fEvents);
	253	for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
	254	TH1F* hEdist = (TH1F*)fOutputList->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
	255	if(fEvents)
	256	hEdist->Scale(1./(Float_t)fEvents);
	257
	258
	259	}
	260	}
	261	}
c48a797f	262
	263	}
	264	//____________________________________________________________________
9f55be54	265	void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption) {
	266
	267	//speciesOption:
	268	//0: p+p Landau fit
	269	//1: Pb+Pb triple landau convolution fit
	270
	271	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	272	pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
	273
	274	TFile fin(filename,"UPDATE");
	275
	276	TList* list = (TList*)fin.Get("energyDist");
	277
	278	AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
	279
	280	EnergyDist->SetNetaBins(pars->GetNetaBins());
	281	EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
	282
884cadc9	283	TF1* fitFunc = 0;
9f55be54	284
884cadc9	285	for(Int_t det = 1; det<=3; det++) {
	286	Int_t nRings = (det == 1 ? 1 : 2);
	287	for(Int_t ring = 0;ring<nRings; ring++) {
	288	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	289	TH1F* hEmptyEdist = (TH1F*)list->FindObject(Form("emptyFMD%d%c",det,ringChar));
	290	TH1F* hRingEdist = (TH1F*)list->FindObject(Form("ringFMD%d%c",det,ringChar));
	291	fitFunc = FitEnergyDistribution(hEmptyEdist, speciesOption) ;
	292	if(fitFunc)
	293	fitFunc->Write(Form("emptyFMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
	294	fitFunc = FitEnergyDistribution(hRingEdist, speciesOption) ;
	295	if(fitFunc)
	296	fitFunc->Write(Form("FMD%d%c_fitfunc",det,ringChar),TObject::kWriteDelete);
	297
	298
	299	EnergyDist->SetEmptyEnergyDistribution(det,ringChar,hEmptyEdist);
	300	EnergyDist->SetRingEnergyDistribution(det,ringChar,hRingEdist);
	301	for(Int_t nEta = 0; nEta < pars->GetNetaBins(); nEta++) {
9f55be54	302	TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
9f55be54	303
884cadc9	304	fitFunc = FitEnergyDistribution(hEdist, speciesOption) ;
884cadc9	305	if(fitFunc)
9f55be54	306	fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
884cadc9	307	EnergyDist->SetEnergyDistribution(det,ringChar,nEta,hEdist);
9f55be54	308
9f55be54	309	}
9f55be54	310
c48a797f	311	}
bb8a464f	312
bb8a464f	313	}
bb8a464f	314
9f55be54	315	fin.Write();
	316	fin.Close();
	317
	318	if(store) {
0b0a4ae5	319	TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
9f55be54	320	EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
	321	fcalib.Close();
	322	}
	323
	324
	325
bb8a464f	326	}
884cadc9	327	//____________________________________________________________________
	328	TF1* AliFMDAnalysisTaskCollector::FitEnergyDistribution(TH1F* hEnergy, Int_t speciesOption) {
	329
	330	TF1* fitFunc = 0;
	331	if(hEnergy->GetEntries() != 0) {
	332
	333	hEnergy->GetXaxis()->SetRangeUser(0.2,hEnergy->GetXaxis()->GetXmax());
	334
	335	if(speciesOption == 0)
	336	fitFunc = new TF1("FMDfitFunc","landau",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
	337	if(speciesOption == 1) {
	338	fitFunc = new TF1("FMDfitFunc",TripleLandau,hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,5,5);
	339	fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
	340	fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
	341	fitFunc->SetParLimits(1,0.6,1.2);
	342	fitFunc->SetParLimits(3,0,1);
	343	fitFunc->SetParLimits(4,0,1);
	344
	345	}
	346	hEnergy->Fit(fitFunc,"","",hEnergy->GetBinCenter(hEnergy->GetMaximumBin())-0.2,3);
	347
	348	}
	349	return fitFunc;
	350
	351	}
c48a797f	352	//____________________________________________________________________
	353	//
	354	// EOF
	355	//