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