[u/mrichter/AliRoot.git] / PWGLF / 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 "AliFMDAnaCalibBackgroundCorrection.h"
#include "AliFMDAnaCalibEnergyDistribution.h"
#include "AliFMDAnaCalibEventSelectionEfficiency.h"
#include "AliFMDAnaCalibSharingEfficiency.h"

//#include "AliFMDGeometry.h"


ClassImp(AliFMDAnalysisTaskCollector)

//____________________________________________________________________
Double_t  AliFMDAnalysisTaskCollector::TripleLandau(const Double_t *x, Double_t *par) {
  //A convolution of three landaus to fit three MIP peaks

  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),
  fClusters(0),
  fClustersEmpty(0),
  fFirstEvent(kTRUE),
  fParam(0)
{
  // Default constructor
  
 
}
//____________________________________________________________________
AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
    AliAnalysisTaskSE(name),
    //fDebug(0),
    fOutputList(0),
    fArray(0),
    fZvtxDist(0),
    fEvents(0),
    fEmptyEvents(0),
    fClusters(0),
    fClustersEmpty(0),
    fFirstEvent(kTRUE),
    fParam(0)
{
  // Default constructor
  
  DefineOutput(1, TList::Class());
  fParam = AliFMDAnaParameters::Instance();
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
{
  // Create the output container
  printf("AnalysisTaskFMD::CreateOutPutData() \n");
  fFirstEvent = kTRUE;
  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*/)
{
  //Collect data for fitting
  AliESDEvent* esd = (AliESDEvent*)InputEvent();
  AliESD* old = esd->GetAliESDOld();
  if (old) {
    esd->CopyFromOldESD();
  }
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  pars->SetTriggerStatus(esd);
  if (fFirstEvent) { 
    pars->SetParametersFromESD(esd);
    pars->PrintStatus();
    fFirstEvent = kFALSE;
  }

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