From e6f3f2fe2f5e41936ccbfc0d3dd013ddbf26dc8c Mon Sep 17 00:00:00 2001 From: loizides Date: Tue, 2 Nov 2010 14:45:30 +0000 Subject: [PATCH] Updates from Alberica to run centrality tasks. --- PWG2/EVCHAR/AliAnalysisTaskCentrality.cxx | 527 ++++++++++++++++++++++ PWG2/EVCHAR/AliAnalysisTaskCentrality.h | 130 ++++++ PWG2/EVCHAR/AliCentralityBy1D.cxx | 111 +++++ PWG2/EVCHAR/AliCentralityBy1D.h | 43 ++ PWG2/EVCHAR/AliCentralityByFunction.cxx | 210 +++++++++ PWG2/EVCHAR/AliCentralityByFunction.h | 51 +++ PWG2/EVCHAR/AliCentralityGlauberFit.cxx | 288 ++++++++++++ PWG2/EVCHAR/AliCentralityGlauberFit.h | 68 +++ PWG2/EVCHAR/runAliCentralityBy1D.C | 24 + PWG2/EVCHAR/runAliCentralityByFunction.C | 44 ++ PWG2/EVCHAR/runAliCentralityGlauberFit.C | 25 + PWG2/EVCHAR/runCentrality.C | 25 +- 12 files changed, 1533 insertions(+), 13 deletions(-) create mode 100644 PWG2/EVCHAR/AliAnalysisTaskCentrality.cxx create mode 100644 PWG2/EVCHAR/AliAnalysisTaskCentrality.h create mode 100644 PWG2/EVCHAR/AliCentralityBy1D.cxx create mode 100644 PWG2/EVCHAR/AliCentralityBy1D.h create mode 100644 PWG2/EVCHAR/AliCentralityByFunction.cxx create mode 100644 PWG2/EVCHAR/AliCentralityByFunction.h create mode 100644 PWG2/EVCHAR/AliCentralityGlauberFit.cxx create mode 100644 PWG2/EVCHAR/AliCentralityGlauberFit.h create mode 100644 PWG2/EVCHAR/runAliCentralityBy1D.C create mode 100644 PWG2/EVCHAR/runAliCentralityByFunction.C create mode 100644 PWG2/EVCHAR/runAliCentralityGlauberFit.C diff --git a/PWG2/EVCHAR/AliAnalysisTaskCentrality.cxx b/PWG2/EVCHAR/AliAnalysisTaskCentrality.cxx new file mode 100644 index 00000000000..bbe7318ac9b --- /dev/null +++ b/PWG2/EVCHAR/AliAnalysisTaskCentrality.cxx @@ -0,0 +1,527 @@ +/************************************************************************** + * Copyright(c) 1998-2008, 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. * + **************************************************************************/ + +///////////////////////////////////////////////////////////// +// // +// Class to analyze centrality measurements // +// // +///////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliAnalysisManager.h" +#include "AliVEvent.h" +#include "AliESD.h" +#include "AliESDEvent.h" +#include "AliESDHeader.h" +#include "AliESDInputHandler.h" +#include "AliESDZDC.h" +#include "AliESDFMD.h" +#include "AliESDVZERO.h" +#include "AliMultiplicity.h" +#include "AliAODHandler.h" +#include "AliAODEvent.h" +#include "AliAODVertex.h" +#include "AliAODMCHeader.h" +#include "AliMCEvent.h" +#include "AliMCEventHandler.h" +#include "AliMCParticle.h" +#include "AliStack.h" +#include "AliHeader.h" +#include "AliAODMCParticle.h" +#include "AliAnalysisTaskSE.h" +#include "AliGenEventHeader.h" +#include "AliGenHijingEventHeader.h" +#include "AliPhysicsSelectionTask.h" +#include "AliPhysicsSelection.h" +#include "AliBackgroundSelection.h" +#include "AliAnalysisTaskCentrality.h" + +ClassImp(AliAnalysisTaskCentrality) + + +//________________________________________________________________________ +AliAnalysisTaskCentrality::AliAnalysisTaskCentrality(): + AliAnalysisTaskSE(), + fDebug(0), + fAnalysisInput("ESD"), + fIsMCInput(kFALSE), + fOutput(0x0), + hEzdc(0x0), + hEzem(0x0), + hNtracks(0x0), + hNtracklets(0x0), + hNclusters0(0x0), + hmultV0(0x0), + hmultFMD(0x0), + hEzemvsEzdc(0x0), + hNtracksvsEzdc(0x0), + hNtrackletsvsEzdc(0x0), + hNclusters0vsEzdc(0x0), + hmultV0vsEzdc(0x0), + hmultFMDvsEzdc(0x0), + hNtracksvsEzem(0x0), + hNtrackletsvsEzem(0x0), + hNclusters0vsEzem(0x0), + hmultV0vsEzem(0x0), + hmultFMDvsEzem(0x0), + hNtracksvsmultV0(0x0), + hNtrackletsvsmultV0(0x0), + hNclusters0vsmultV0(0x0), + hNtracksvsmultFMD(0x0), + hNtrackletsvsmultFMD(0x0), + hNclusters0vsmultFMD(0x0), + hmultV0vsmultFMD(0x0) +{ + // Default constructor +} + +//________________________________________________________________________ +AliAnalysisTaskCentrality::AliAnalysisTaskCentrality(const char *name): + AliAnalysisTaskSE(name), + fDebug(0), + fAnalysisInput("ESD"), + fIsMCInput(kFALSE), + fOutput(0x0), + hEzdc(0x0), + hEzem(0x0), + hNtracks(0x0), + hNtracklets(0x0), + hNclusters0(0x0), + hmultV0(0x0), + hmultFMD(0x0), + hEzemvsEzdc(0x0), + hNtracksvsEzdc(0x0), + hNtrackletsvsEzdc(0x0), + hNclusters0vsEzdc(0x0), + hmultV0vsEzdc(0x0), + hmultFMDvsEzdc(0x0), + hNtracksvsEzem(0x0), + hNtrackletsvsEzem(0x0), + hNclusters0vsEzem(0x0), + hmultV0vsEzem(0x0), + hmultFMDvsEzem(0x0), + hNtracksvsmultV0(0x0), + hNtrackletsvsmultV0(0x0), + hNclusters0vsmultV0(0x0), + hNtracksvsmultFMD(0x0), + hNtrackletsvsmultFMD(0x0), + hNclusters0vsmultFMD(0x0), + hmultV0vsmultFMD(0x0) +{ + // Default constructor + + // Output slot #1 writes into a TList container + DefineOutput(1, TList::Class()); + +} + +//________________________________________________________________________ +AliAnalysisTaskCentrality& AliAnalysisTaskCentrality::operator=(const AliAnalysisTaskCentrality& c) +{ + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c); + } + return *this; +} + +//________________________________________________________________________ +AliAnalysisTaskCentrality::AliAnalysisTaskCentrality(const AliAnalysisTaskCentrality& ana): + AliAnalysisTaskSE(ana), + fDebug(ana.fDebug), + fAnalysisInput(ana.fDebug), + fIsMCInput(ana.fIsMCInput), + fOutput(ana.fOutput), + hEzdc(ana.hEzdc), + hEzem(ana.hEzem), + hNtracks(ana.hNtracks), + hNtracklets(ana.hNtracklets), + hNclusters0(ana.hNclusters0), + hmultV0(ana.hmultV0), + hmultFMD(ana.hmultFMD), + hEzemvsEzdc(ana.hEzemvsEzdc), + hNtracksvsEzdc(ana.hNtracksvsEzdc), + hNtrackletsvsEzdc(ana.hNtrackletsvsEzdc), + hNclusters0vsEzdc(ana.hNclusters0vsEzdc), + hmultV0vsEzdc(ana.hmultV0vsEzdc), + hmultFMDvsEzdc(ana.hmultFMDvsEzdc), + hNtracksvsEzem(ana.hNtracksvsEzem), + hNtrackletsvsEzem(ana.hNtrackletsvsEzem), + hNclusters0vsEzem(ana.hNclusters0vsEzem), + hmultV0vsEzem(ana.hmultV0vsEzem), + hmultFMDvsEzem(ana.hmultFMDvsEzem), + hNtracksvsmultV0(ana.hNtracksvsmultV0), + hNtrackletsvsmultV0(ana.hNtrackletsvsmultV0), + hNclusters0vsmultV0(ana.hNclusters0vsmultV0), + hNtracksvsmultFMD(ana.hNtracksvsmultFMD), + hNtrackletsvsmultFMD(ana.hNtrackletsvsmultFMD), + hNclusters0vsmultFMD(ana.hNclusters0vsmultFMD), + hmultV0vsmultFMD(ana.hmultV0vsmultFMD) +{ + // + // Copy Constructor + // +} + +//________________________________________________________________________ + AliAnalysisTaskCentrality::~AliAnalysisTaskCentrality() + { + // Destructor + if(fOutput){ + delete fOutput; fOutput=0; + } + } + +//________________________________________________________________________ +void AliAnalysisTaskCentrality::UserCreateOutputObjects() +{ + + // Create the output containers + if(fDebug>1) printf("AnalysisTaskZDCpp::UserCreateOutputObjects() \n"); + + // Several histograms are more conveniently managed in a TList + fOutput = new TList(); + fOutput->SetOwner(); + fOutput->SetName("OutputHistos"); + + hEzdc = new TH1F("hEzdc","hEzdc",500,0,150); + hEzem = new TH1F("hEzem","hEzem",500,0,5); + hNtracks = new TH1F("hNtracks","hNtracks",500,0,17000); + hNtracklets = new TH1F("hNtracklets","hNtracklets",500,0,10000); + hNclusters0 = new TH1F("hNclusters0","hNclusters0",500,0,15000); + hmultV0 = new TH1F("hmultV0","hmultV0",500,0,30000); + hmultFMD = new TH1F("hmultFMD","hmultFMD",500,0,24000); + + hEzemvsEzdc = new TProfile("hEzemvsEzdc","hEzemvsEzdc",500,0,5,""); + hNtracksvsEzdc = new TProfile("hNtracksvsEzdc","hNtracksvsEzdc",500,0,17000,""); + hNtrackletsvsEzdc = new TProfile("hNtrackletsvsEzdc","hNtrackletsvsEzdc",500,0,10000,""); + hNclusters0vsEzdc = new TProfile("hNclusters0vsEzdc","hNclusters0vsEzdc",500,0,15000,""); + hmultV0vsEzdc = new TProfile("hmultV0vsEzdc","hmultV0vsEzdc",500,0,30000,""); + hmultFMDvsEzdc = new TProfile("hmultFMDvsEzdc","hmultFMDvsEzdc",500,0,24000,""); + hNtracksvsEzem = new TProfile("hNtracksvsEzem","hNtracksvsEzem",500,0,17000,""); + hNtrackletsvsEzem = new TProfile("hNtrackletsvsEzem","hNtrackletsvsEzem",500,0,10000,""); + hNclusters0vsEzem = new TProfile("hNclusters0vsEzem","hNclusters0vsEzem",500,0,15000,""); + hmultV0vsEzem = new TProfile("hmultV0vsEzem","hmultV0vsEzem",500,0,30000,""); + hmultFMDvsEzem = new TProfile("hmultFMDvsEzem","hmultFMDvsEzem",500,0,24000,""); + hNtracksvsmultV0 = new TProfile("hNtracksvsmultV0","hNtracksvsmultV0",500,0,17000,""); + hNtrackletsvsmultV0 = new TProfile("hNtrackletsvsmultV0","hNtrackletsvsmultV0",500,0,10000,""); + hNclusters0vsmultV0 = new TProfile("hNclusters0vsmultV0","hNclusters0vsmultV0",500,0,15000,""); + hNtracksvsmultFMD = new TProfile("hNtracksvsmultFMD","hNtracksvsmultFMD",500,0,17000,""); + hNtrackletsvsmultFMD= new TProfile("hNtrackletsvsmultFMD","hNtrackletsvsmultFMD",500,0,10000,""); + hNclusters0vsmultFMD= new TProfile("hNclusters0vsmultFMD","hNclusters0vsmultFMD",500,0,15000,""); + hmultV0vsmultFMD = new TProfile("hmultV0vsmultFMD","hmultV0vsmultFMD",500,0,30000,""); + + hEzdc ->GetXaxis()->SetTitle("E_{ZDC}[TeV]"); + hEzem ->GetXaxis()->SetTitle("E_{ZEM}[TeV]"); + hNtracks ->GetXaxis()->SetTitle("N_{tracks}"); + hNtracklets ->GetXaxis()->SetTitle("N_{tracklets}"); + hNclusters0 ->GetXaxis()->SetTitle("N_{clusters0}"); + hmultV0 ->GetXaxis()->SetTitle("V0 mult"); + hmultFMD ->GetXaxis()->SetTitle("FMD mult"); + + hEzemvsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hNtracksvsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hNtrackletsvsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hNclusters0vsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hmultV0vsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hmultFMDvsEzdc ->GetYaxis()->SetTitle("E_{ZDC}[TeV]"); + hNtracksvsEzem ->GetYaxis()->SetTitle("E_{ZEM}[TeV]"); + hNtrackletsvsEzem ->GetYaxis()->SetTitle("E_{ZEM}[TeV]"); + hNclusters0vsEzem ->GetYaxis()->SetTitle("E_{ZEM}[TeV]"); + hmultV0vsEzem ->GetYaxis()->SetTitle("E_{ZEM}[TeV]"); + hmultFMDvsEzem ->GetYaxis()->SetTitle("E_{ZEM}[TeV]"); + hNtracksvsmultV0 ->GetYaxis()->SetTitle("V0 mult"); + hNtrackletsvsmultV0 ->GetYaxis()->SetTitle("V0 mult"); + hNclusters0vsmultV0 ->GetYaxis()->SetTitle("V0 mult"); + hNtracksvsmultFMD ->GetYaxis()->SetTitle("FMD mult"); + hNtrackletsvsmultFMD->GetYaxis()->SetTitle("FMD mult"); + hNclusters0vsmultFMD->GetYaxis()->SetTitle("FMD mult"); + hmultV0vsmultFMD ->GetYaxis()->SetTitle("FMD mult"); + + hEzemvsEzdc ->GetXaxis()->SetTitle("E_{ZEM}[TeV]"); + hNtracksvsEzdc ->GetXaxis()->SetTitle("N_{tracks}"); + hNtrackletsvsEzdc ->GetXaxis()->SetTitle("N_{tracklets}"); + hNclusters0vsEzdc ->GetXaxis()->SetTitle("N_{clusters0}"); + hmultV0vsEzdc ->GetXaxis()->SetTitle("V0 mult"); + hmultFMDvsEzdc ->GetXaxis()->SetTitle("FMD mult"); + hNtracksvsEzem ->GetXaxis()->SetTitle("N_{tracks}"); + hNtrackletsvsEzem ->GetXaxis()->SetTitle("N_{tracklets}"); + hNclusters0vsEzem ->GetXaxis()->SetTitle("N_{clusters0}"); + hmultV0vsEzem ->GetXaxis()->SetTitle("V0 mult"); + hmultFMDvsEzem ->GetXaxis()->SetTitle("FMD mult"); + hNtracksvsmultV0 ->GetXaxis()->SetTitle("N_{tracks}"); + hNtrackletsvsmultV0 ->GetXaxis()->SetTitle("N_{tracklets}"); + hNclusters0vsmultV0 ->GetXaxis()->SetTitle("N_{clusters0}"); + hNtracksvsmultFMD ->GetXaxis()->SetTitle("N_{tracks}"); + hNtrackletsvsmultFMD->GetXaxis()->SetTitle("N_{tracklets}"); + hNclusters0vsmultFMD->GetXaxis()->SetTitle("N_{clusters}"); + hmultV0vsmultFMD ->GetXaxis()->SetTitle("V0 mult"); + + fOutput->Add(hEzdc); + fOutput->Add(hEzem); + fOutput->Add(hNtracks); + fOutput->Add(hNtracklets); + fOutput->Add(hNclusters0); + fOutput->Add(hmultV0); + fOutput->Add(hmultFMD); + + fOutput->Add(hEzemvsEzdc); + fOutput->Add(hNtracksvsEzdc); + fOutput->Add(hNtrackletsvsEzdc); + fOutput->Add(hNclusters0vsEzdc); + fOutput->Add(hmultV0vsEzdc); + fOutput->Add(hmultFMDvsEzdc); + fOutput->Add(hNtracksvsEzem); + fOutput->Add(hNtrackletsvsEzem); + fOutput->Add(hNclusters0vsEzem); + fOutput->Add(hmultV0vsEzem); + fOutput->Add(hmultFMDvsEzem); + fOutput->Add(hNtracksvsmultV0); + fOutput->Add(hNtrackletsvsmultV0); + fOutput->Add(hNclusters0vsmultV0); + fOutput->Add(hNtracksvsmultFMD); + fOutput->Add(hNtrackletsvsmultFMD); + fOutput->Add(hNclusters0vsmultFMD); + fOutput->Add(hmultV0vsmultFMD); + + PostData(1, fOutput); + +} + +//________________________________________________________________________ +void AliAnalysisTaskCentrality::UserExec(Option_t */*option*/) +{ + // Execute analysis for current event: + if(fDebug>1) printf(" **** AliAnalysisTaskCentrality::UserExec() \n"); + + if(fAnalysisInput.CompareTo("ESD")==0){ + + AliVEvent* event = InputEvent(); + AliESDEvent* esd = dynamic_cast(event); + + fNev++; + + fNTracks = event->GetNumberOfTracks(); + fNPmdTracks = esd->GetNumberOfPmdTracks(); + + AliESDVZERO* esdV0 = esd->GetVZEROData(); + fMultV0A=esdV0->GetMTotV0A(); + fMultV0C=esdV0->GetMTotV0C(); + + if(fIsMCInput){ + + AliMCEvent* mcEvent = MCEvent(); + if (!mcEvent) { + printf(" Could not retrieve MC event!!!\n"); + return; + } + + fNmyTracks_gen = 0; + AliStack *stack = 0x0; // needed for MC studies + stack = MCEvent()->Stack(); + for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) { + //get properties of mc particle + AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack); + // Primaries only + if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue; + //charged tracks only + if (mcP->Particle()->GetPDG()->Charge() == 0) continue; + //same cuts as on ESDtracks +// if(TMath::Abs(mcP->Eta())>0.9)continue; +// if(mcP->Pt()<0.2)continue; +// if(mcP->Pt()>200)continue; + + fNmyTracks_gen ++; + } + + AliGenEventHeader* genHeader = mcEvent->GenEventHeader(); + if(!genHeader){ + printf(" Event generator header not available!!!\n"); + return; + } + + if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){ + fbMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter(); + Int_t specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn(); + Int_t specProtonProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp(); + Int_t specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn(); + Int_t specProtonTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp(); + fNpartTargMC = 208.-(specNeutronTarg+specProtonTarg); + fNpartProjMC = 208.-(specNeutronProj+specProtonProj); + fNNColl = ((AliGenHijingEventHeader*) genHeader)->NN(); + fNNwColl = ((AliGenHijingEventHeader*) genHeader)->NNw(); + fNwNColl = ((AliGenHijingEventHeader*) genHeader)->NwN(); + fNwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw(); + } + + } + + fBeamEnergy = esd->GetBeamEnergy(); + + // ***** Trigger selection + TString triggerClass = esd->GetFiredTriggerClasses(); + sprintf(fTrigClass,"%s",triggerClass.Data()); + + const AliESDVertex *vertex = esd->GetPrimaryVertexSPD(); + fxVertex = vertex->GetX(); + fyVertex = vertex->GetY(); + fzVertex = vertex->GetZ(); + if(vertex->IsFromVertexer3D()) fVertexer3d = kTRUE; + else fVertexer3d = kFALSE; + Double_t vertex3[3]; + vertex->GetXYZ(vertex3); + + const AliMultiplicity *mult = esd->GetMultiplicity(); + fNTracklets = mult->GetNumberOfTracklets(); + + for(Int_t ilay=0; ilay<6; ilay++){ + fNClusters[ilay] = mult->GetNumberOfITSClusters(ilay); + } + fNSingleClusters = mult->GetNumberOfSingleClusters(); + + for(Int_t ilay=0; ilay<2; ilay++){ + fNChips[ilay] = mult->GetNumberOfFiredChips(ilay); + } + + + AliESDFMD *fmd = esd->GetFMDData(); + Float_t totalMultA = 0; + Float_t totalMultC = 0; + const Float_t fFMDLowCut = 0.4; + + 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'); + 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 FMDmult = fmd->Multiplicity(det,ring,sec,strip); + if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue; + + Float_t nParticles=0; + + if(FMDmult > fFMDLowCut) { + nParticles = 1.; + } + + if (det<3) totalMultA = totalMultA + nParticles; + else totalMultC = totalMultC + nParticles; + + } + } + } + } + fMultFMDA = totalMultA; + fMultFMDC = totalMultC; + + + AliESDZDC *esdZDC = esd->GetESDZDC(); + + fESDFlag = esdZDC->GetESDQuality(); + + fZNCEnergy = (Float_t) (esdZDC->GetZDCN1Energy()); + fZPCEnergy = (Float_t) (esdZDC->GetZDCP1Energy()); + fZNAEnergy = (Float_t) (esdZDC->GetZDCN2Energy()); + fZPAEnergy = (Float_t) (esdZDC->GetZDCP2Energy()); + fZEM1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0)); + fZEM2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1)); + + fbZDC = esdZDC->GetImpactParameter(); + fNpartZDC = esdZDC->GetZDCParticipants(); + fbZDCA = esdZDC->GetImpactParamSideA(); + fNpartZDCA = esdZDC->GetZDCPartSideA(); + fbZDCC = esdZDC->GetImpactParamSideC(); + fNpartZDCC = esdZDC->GetZDCPartSideC(); + + const Double_t * towZNC = esdZDC->GetZN1TowerEnergy(); + const Double_t * towZPC = esdZDC->GetZP1TowerEnergy(); + const Double_t * towZNA = esdZDC->GetZN2TowerEnergy(); + const Double_t * towZPA = esdZDC->GetZP2TowerEnergy(); + + for(Int_t it=0; it<5; it++){ + fZNCtower[it] = (Float_t) (towZNC[it]); + fZPCtower[it] = (Float_t) (towZPC[it]); + fZNAtower[it] = (Float_t) (towZNA[it]); + fZPAtower[it] = (Float_t) (towZPA[it]); + } + + Double_t xyZNC[2]={-99.,-99.}, xyZNA[2]={-99.,-99.}; + esdZDC->GetZNCentroidInPbPb(fBeamEnergy, xyZNC, xyZNA); + for(Int_t it=0; it<2; it++){ + fCentrZNC[it] = xyZNC[it]; + fCentrZNA[it] = xyZNA[it]; + } + + // filling histos + + hEzdc ->Fill((fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hEzem ->Fill(fZEM1Energy+fZEM2Energy); + hNtracks ->Fill(fNTracks); + hNtracklets ->Fill(fNTracklets); + hNclusters0 ->Fill(fNClusters[0]); + hmultV0 ->Fill(fMultV0A+fMultV0C); + hmultFMD ->Fill(fMultFMDA+fMultFMDC); + + hEzemvsEzdc ->Fill(fZEM1Energy+fZEM2Energy, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hNtracksvsEzdc ->Fill(fNTracks, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hNtrackletsvsEzdc ->Fill(fNTracklets, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hNclusters0vsEzdc ->Fill(fNClusters[0], (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hmultV0vsEzdc ->Fill(fMultV0A+fMultV0C, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hmultFMDvsEzdc ->Fill(fMultFMDA+fMultFMDC, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.); + hNtracksvsEzem ->Fill(fNTracks, fZEM1Energy+fZEM2Energy); + hNtrackletsvsEzem ->Fill(fNTracklets, fZEM1Energy+fZEM2Energy); + hNclusters0vsEzem ->Fill(fNClusters[0], fZEM1Energy+fZEM2Energy); + hmultV0vsEzem ->Fill(fMultV0A+fMultV0C, fZEM1Energy+fZEM2Energy); + hmultFMDvsEzem ->Fill(fMultFMDA+fMultFMDC, fZEM1Energy+fZEM2Energy); + hNtracksvsmultV0 ->Fill(fNTracks,fMultV0A+fMultV0C); + hNtrackletsvsmultV0 ->Fill(fNTracklets,fMultV0A+fMultV0C); + hNclusters0vsmultV0 ->Fill(fNClusters[0],fMultV0A+fMultV0C); + hNtracksvsmultFMD ->Fill(fNTracks,fMultFMDA+fMultFMDC); + hNtrackletsvsmultFMD->Fill(fNTracklets,fMultFMDA+fMultFMDC); + hNclusters0vsmultFMD->Fill(fNClusters[0],fMultFMDA+fMultFMDC); + hmultV0vsmultFMD ->Fill(fMultV0A+fMultV0C,fMultFMDA+fMultFMDC); + + } + else if(fAnalysisInput.CompareTo("AOD")==0){ + //AliAODEvent *aod = dynamic_cast (InputEvent()); + // to be implemented + printf(" AOD analysis not yet implemented!!!\n\n"); + return; + } + PostData(1, fOutput); +} + + + +//________________________________________________________________________ +void AliAnalysisTaskCentrality::Terminate(Option_t */*option*/) +{ + // Terminate analysis +} + + diff --git a/PWG2/EVCHAR/AliAnalysisTaskCentrality.h b/PWG2/EVCHAR/AliAnalysisTaskCentrality.h new file mode 100644 index 00000000000..2b2d935e57b --- /dev/null +++ b/PWG2/EVCHAR/AliAnalysisTaskCentrality.h @@ -0,0 +1,130 @@ +#ifndef ALIANALYSISTASKCENTRALITY_H +#define ALIANALYSISTASKCENTRALITY_H + +/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//***************************************************** +// Class AliAnalysisTaskCentrality +// author: Alberica Toia +//***************************************************** + +#include "AliAnalysisTaskSE.h" + +class TROOT; +class TSystem; +class TList; +class TFile; +class TH2F; + +class AliAnalysisTaskCentrality : public AliAnalysisTaskSE { + + public: + + AliAnalysisTaskCentrality(); + AliAnalysisTaskCentrality(const char *name); + AliAnalysisTaskCentrality& operator= (const AliAnalysisTaskCentrality& ana); + AliAnalysisTaskCentrality(const AliAnalysisTaskCentrality& c); + virtual ~AliAnalysisTaskCentrality(); + + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + + virtual void SetDebugLevel(Int_t level) {fDebug = level;} + void SetInput(const char* input) {fAnalysisInput = input;} + void SetMCInput() {fIsMCInput = kTRUE;} + + + private: + + Int_t fDebug; // Debug flag + TString fAnalysisInput; // "ESD", "AOD" + Bool_t fIsMCInput; // true when input is MC + // + TList *fOutput; //! list send on output slot 0 + TH1F *hEzdc; //! + TH1F *hEzem; //! + TH1F *hNtracks; //! + TH1F *hNtracklets; //! + TH1F *hNclusters0; //! + TH1F *hmultV0; //! + TH1F *hmultFMD; //! + + TProfile *hEzemvsEzdc; //! + TProfile *hNtracksvsEzdc; //! + TProfile *hNtrackletsvsEzdc; //! + TProfile *hNclusters0vsEzdc; //! + TProfile *hmultV0vsEzdc; //! + TProfile *hmultFMDvsEzdc; //! + TProfile *hNtracksvsEzem; //! + TProfile *hNtrackletsvsEzem; //! + TProfile *hNclusters0vsEzem; //! + TProfile *hmultV0vsEzem; //! + TProfile *hmultFMDvsEzem; //! + TProfile *hNtracksvsmultV0; //! + TProfile *hNtrackletsvsmultV0; //! + TProfile *hNclusters0vsmultV0; //! + TProfile *hNtracksvsmultFMD; //! + TProfile *hNtrackletsvsmultFMD; //! + TProfile *hNclusters0vsmultFMD; //! + TProfile *hmultV0vsmultFMD; //! + + Int_t fNev; // event counter + Float_t fBeamEnergy; // beam energy + Int_t fNmyTracks_gen; // no. generated primary charged tracks + char fTrigClass[100]; // fired trigger classes + // + Double_t fxVertex; // X vertex from ITS + Double_t fyVertex; // Y vertex from ITS + Double_t fzVertex; // Z vertex from ITS + Bool_t fVertexer3d; // Is vertex from 3d vertexer? + // + Double_t fbMC; // impact parameter from MC + Int_t fNpartTargMC; // no. of participants for target nucleus from MC + Int_t fNpartProjMC; // no. of participants for project nucleus from MC + Int_t fNNColl; // Number of N-N collisions + Int_t fNNwColl; // Number of N-Nwounded collisions + Int_t fNwNColl; // Number of Nwounded-N collisons + Int_t fNwNwColl; // Number of Nwounded-Nwounded collisions + // + Int_t fNTracklets; // no. tracklets + Int_t fNSingleClusters; // no. single clusters + Int_t fNClusters[6]; // no. clusters on 6 ITS layers + Int_t fNChips[2]; // no. chips on 2 SPD layers + // + Double_t fbZDC; // impact parameter from ZDC reco + Int_t fNpartZDC; // no. of participants from ZDC reco + Double_t fbZDCA; // impact parameter from ZDC reco side A + Int_t fNpartZDCA; // no. of part. from ZDC reco side A + Double_t fbZDCC; // impact parameter from ZDC reco side C + Int_t fNpartZDCC; // no. of part. from ZDC reco side C + // + UInt_t fESDFlag; // ZDC ESD flags + Float_t fZNCEnergy; // ZNC Energy + Float_t fZPCEnergy; // ZPC Energy + Float_t fZNAEnergy; // ZNA Energy + Float_t fZPAEnergy; // ZPA Energy + Float_t fZEM1Energy; // ZEM1 Energy + Float_t fZEM2Energy; // ZEM2 Energy + Float_t fZNCtower[5]; // ZNC 5 tower signals + Float_t fZPCtower[5]; // ZPC 5 tower signals + Float_t fZNAtower[5]; // ZNA 5 tower signals + Float_t fZPAtower[5]; // ZPA 5 tower signals + Float_t fCentrZNC[2]; // centroid over ZNC + Float_t fCentrZNA[2]; // centroid over ZNA + + Int_t fNTracks; // no. tracks + Int_t fNPmdTracks; // no. PMD tracks + Float_t fMultV0A; // multiplicity from V0 reco side A + Float_t fMultV0C; // multiplicity from V0 reco side C + Float_t fMultFMDA; // multiplicity from FMD on detector A + Float_t fMultFMDC; // multiplicity from FMD on detector C + + ClassDef(AliAnalysisTaskCentrality,1); + +}; + +#endif + diff --git a/PWG2/EVCHAR/AliCentralityBy1D.cxx b/PWG2/EVCHAR/AliCentralityBy1D.cxx new file mode 100644 index 00000000000..07257beb3e9 --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityBy1D.cxx @@ -0,0 +1,111 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 1D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliCentralityBy1D.h" + + +ClassImp(AliCentralityBy1D) + +//______________________________________________________________________________ + + +AliCentralityBy1D::AliCentralityBy1D() { + // standard constructor +} + +AliCentralityBy1D::~AliCentralityBy1D() { + // destructor +} + +void AliCentralityBy1D::AddHisto(TString name) { + histnames.push_back(name); +} + +void AliCentralityBy1D::SetPercentileFile(TString filename) { + outrootfilename = filename; +} + +void AliCentralityBy1D::SetPercentileCrossSection(Float_t xsec) { + percentXsec = xsec; +} + +void AliCentralityBy1D::MakePercentiles(TString infilename) { + TH1D *thist; + + TFile *outrootfile; + + // open inrootfile, outrootfile + inrootfile = new TFile(infilename); + outrootfile = new TFile(outrootfilename,"RECREATE"); + + // loop over all distribution names + + vector::const_iterator hni; + for(hni=histnames.begin(); hni!=histnames.end(); hni++) { + thist = MakePercentHisto(*hni); + SaveHisto(thist,outrootfile); + delete thist; //?? + } + // close inrootfile, outrootfile + inrootfile->Close(); + outrootfile->Close(); + +} + +TH1D * AliCentralityBy1D::MakePercentHisto(TString hdistributionName) { + TH1D *htemp = (TH1D*) (inrootfile->Get(hdistributionName)); + TH1D *hpercent = (TH1D*) htemp->Clone("hpercent"); + hpercent->SetName(hdistributionName.Append("_percentile")); + hpercent->Reset(); + + for (int ibin=1; ibin<=htemp->GetNbinsX(); ibin++) { + + hpercent->SetBinContent(ibin, percentXsec * + htemp->Integral(ibin,htemp->GetNbinsX()) / + htemp->Integral(1,htemp->GetNbinsX())); + + } + + delete htemp; + return hpercent; + +} + +void AliCentralityBy1D::SaveHisto(TH1D *hist, TFile *outrootfile) { + outrootfile->cd(); + hist->Write(); +} + + + diff --git a/PWG2/EVCHAR/AliCentralityBy1D.h b/PWG2/EVCHAR/AliCentralityBy1D.h new file mode 100644 index 00000000000..b49383996df --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityBy1D.h @@ -0,0 +1,43 @@ +#ifndef ALICENTRALITYBY1D_H +#define ALICENTRALITYBY1D_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 1D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include "TObject.h" + +class AliCentralityBy1D : public TObject { + + public: + + AliCentralityBy1D(); + virtual ~AliCentralityBy1D(); + + void SetPercentileFile(TString outrootfilename); + void SetPercentileCrossSection(Float_t percentXsec); + void AddHisto(TString name); + void MakePercentiles(TString infilename); + + private: + + TFile *inrootfile; + + TString outrootfilename; + vector histnames; + Float_t percentXsec; + + TH1D * MakePercentHisto(TString hdistributionName); + void SaveHisto(TH1D *hist, TFile *outrootfile); + + ClassDef(AliCentralityBy1D, 1) +}; +#endif + + diff --git a/PWG2/EVCHAR/AliCentralityByFunction.cxx b/PWG2/EVCHAR/AliCentralityByFunction.cxx new file mode 100644 index 00000000000..fabc87cbe44 --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityByFunction.cxx @@ -0,0 +1,210 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 2D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliCentralityByFunction.h" + + +ClassImp(AliCentralityByFunction) + +//______________________________________________________________________________ + +Double_t fitf_pol2(Double_t* x, Double_t* par) { + Double_t fitValue = + (par[0]+ + par[1]*x[0]+ + par[2]*x[0]*x[0]); + return fitValue; +} +Double_t fitf_pol3(Double_t* x, Double_t* par) { + Double_t fitValue = + (par[0]+ + par[1]*x[0]+ + par[2]*x[0]*x[0]+ + par[3]*x[0]*x[0]*x[0]); + return fitValue; +} +Double_t fitf_pol4(Double_t* x, Double_t* par) { + Double_t fitValue = + (par[0]+ + par[1]*x[0]+ + par[2]*x[0]*x[0]+ + par[3]*x[0]*x[0]*x[0]+ + par[4]*x[0]*x[0]*x[0]*x[0]); + return fitValue; +} +Double_t fitf_pol6(Double_t* x, Double_t* par) { + Double_t fitValue = + (par[0]+ + par[1]*x[0]+ + par[2]*x[0]*x[0]+ + par[3]*x[0]*x[0]*x[0]+ + par[4]*x[0]*x[0]*x[0]*x[0]+ + par[5]*x[0]*x[0]*x[0]*x[0]*x[0]+ + par[6]*x[0]*x[0]*x[0]*x[0]*x[0]*x[0]); + return fitValue; +} + +AliCentralityByFunction::AliCentralityByFunction() { + // standard constructor + fitter["fitf_pol2"] = new TF1("fitf_pol2",fitf_pol2,0,1,3); + fitter["fitf_pol2"]->SetLineColor(kRed); + + fitter["fitf_pol3"] = new TF1("fitf_pol3",fitf_pol3,0,1,4); + fitter["fitf_pol3"]->SetLineColor(kRed); + + fitter["fitf_pol4"] = new TF1("fitf_pol4",fitf_pol4,0,1,5); + fitter["fitf_pol4"]->SetLineColor(kRed); + + fitter["fitf_pol6"] = new TF1("fitf_pol6",fitf_pol6,0,1,6); + fitter["fitf_pol6"]->SetLineColor(kRed); +} + +AliCentralityByFunction::~AliCentralityByFunction() { + // destructor +} + +void AliCentralityByFunction::AddHisto(TString name) { + histnames.push_back(name); +} + +void AliCentralityByFunction::SetPercentileFile(TString filename) { + outrootfilename = filename; +} + +void AliCentralityByFunction::SetPercentileCrossSection(Float_t xsec) { + percentXsec = xsec; +} + +void AliCentralityByFunction::SetFitFunction(TString distribution, TString func, Double_t xmin, Double_t xmax) { + fitfunc[distribution] = func; + fitter[fitfunc[distribution]]->SetRange(xmin,xmax); +} + +void AliCentralityByFunction::MakePercentiles(TString infilename) { + TH1D *hpercentile; + + // open inrootfile, outrootfile + inrootfile = new TFile(infilename); + outrootfile = new TFile(outrootfilename,"RECREATE"); + + // loop over all distribution names + vector::const_iterator hni; + for(hni=histnames.begin(); hni!=histnames.end(); hni++) { + hpercentile = FitHisto(*hni); + outrootfile->cd(); + hpercentile->Write(); + } + // close inrootfile, outrootfile + inrootfile->Close(); + outrootfile->Close(); + +} + + TH1D *AliCentralityByFunction::FitHisto(TString hdistributionName) { + TH2D *hdist = (TH2D*) (inrootfile->Get(hdistributionName)); + TProfile *profile =hdist->ProfileX(); + // fitter[fitfunc[hdistributionName]]->SetRange(0,profile->GetBinCenter(profile->GetNbinsX())); + profile->Fit(fitter[fitfunc[hdistributionName]], "RNM"); + + outrootfile->cd(); + profile->Write(); + fitter[fitfunc[hdistributionName]]->Write(hdistributionName.Append("_fit")); + + return MakePercentHisto(hdist); +} + +TH1D * AliCentralityByFunction::MakePercentHisto(TH2D *histo) { + TH2D *htemp = (TH2D*)histo->Clone("htemp"); + TString hdistributionName = htemp->GetTitle(); + + const int NUM_DIVISION=500; + + TH1D *hpercent = new TH1D("","",NUM_DIVISION,0,htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX())); + hpercent->Reset(); + + double Xpoint, Ypoint, Xcomp, Ycomp, count; + + double Xmax = htemp->GetXaxis()->GetBinCenter(htemp->GetNbinsX()); + double Xmin = 0; + + double delta = (Xmax - Xmin) / NUM_DIVISION; + double slopePerp; // slope of the perpendicular line + double slopeFunction; + + cout << "Start Percentile Histo for distribution " << hdistributionName << " with fit function " << fitfunc[hdistributionName] << endl; + + // move perpendicular line along fitting curve + // and count number of points on the right + for (int i = 0; i <= NUM_DIVISION; i++) { + Xcomp = Xmin + i * delta; + Ycomp = fitter[fitfunc[hdistributionName]]->Eval(Xcomp); + count = 0.0; + slopeFunction = fitter[fitfunc[hdistributionName]]->Derivative(Xcomp,NULL,0.0001); + slopePerp = -1.0 /slopeFunction; + // + // equation of perpendicular line + // (y-Ycomp)/(x-Xcomp) = slopePerp + // + for (int ibiny = 1; ibiny < htemp->GetNbinsY(); ibiny++) { + Ypoint = htemp->GetYaxis()->GetBinCenter(ibiny); + + double XonLine =(Ypoint - Ycomp)/slopePerp + Xcomp; + // double YonLine = slopePerp*(XonLine - Xcomp) +Ycomp; + + for (int ibinx = 1; ibinx < htemp->GetNbinsX(); ibinx++) { + Xpoint = htemp->GetXaxis()->GetBinCenter(ibinx); + // + // (XonLine,Ypoint) lies on the perpendicular + // + if ((Xpoint > XonLine && Xpoint > Xcomp - 0.5)|| + (Xpoint > Xcomp + 0.2)) { + + count += htemp->GetBinContent(ibinx,ibiny); + } + } + } + count = count/htemp->Integral() * 100.0; // change to percentage + hpercent->SetBinContent(i, count); + } + + hpercent->SetName(hdistributionName.Append("_percentile")); + + return hpercent; +} + + + diff --git a/PWG2/EVCHAR/AliCentralityByFunction.h b/PWG2/EVCHAR/AliCentralityByFunction.h new file mode 100644 index 00000000000..e50c58a2c2b --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityByFunction.h @@ -0,0 +1,51 @@ +#ifndef ALICENTRALITYBYFUNCTION_H +#define ALICENTRALITYBYFUNCTION_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 2D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include "TObject.h" + +// forward decl +class TF1; + +class AliCentralityByFunction : public TObject { + + public: + + AliCentralityByFunction(); + virtual ~AliCentralityByFunction(); + + void SetPercentileFile(TString outrootfilename); + void SetPercentileCrossSection(Float_t percentXsec); + // void SetFitFunction(TString distribution, TString func); + void SetFitFunction(TString distribution, TString func, Double_t xmin, Double_t xmax); + void AddHisto(TString name); + void MakePercentiles(TString infilename); + + private: + + TFile *inrootfile; + TString outrootfilename; + TFile *outrootfile; + + vector histnames; + Float_t percentXsec; + mapfitfunc; // mapping from distribution to fit function name + mapfitter; // mapping from fit function name to corresponding TF1 + + TH1D *FitHisto(TString hdistributionName); + TH1D * MakePercentHisto(TH2D *hist); + + ClassDef(AliCentralityByFunction, 1) +}; +#endif + + diff --git a/PWG2/EVCHAR/AliCentralityGlauberFit.cxx b/PWG2/EVCHAR/AliCentralityGlauberFit.cxx new file mode 100644 index 00000000000..5e439ff89c1 --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityGlauberFit.cxx @@ -0,0 +1,288 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 1D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliCentralityGlauberFit.h" + + +ClassImp(AliCentralityGlauberFit) + +//______________________________________________________________________________ + + +AliCentralityGlauberFit::AliCentralityGlauberFit() { + // standard constructor + // glauber file + TFile *f = TFile::Open("hj-unquenched.root"); + glau_ntuple = (TNtuple*) f->Get("gnt"); + glau_ntuple->SetBranchAddress("npart",&npart); + glau_ntuple->SetBranchAddress("ncoll",&ncoll); + glau_ntuple->SetBranchAddress("b",&b); +} + +//-------------------------------------------------------------------------------------------------- + +AliCentralityGlauberFit::~AliCentralityGlauberFit() { + // destructor +} + +//-------------------------------------------------------------------------------------------------- + +void AliCentralityGlauberFit::AddHisto(TString name) { + histnames.push_back(name); +} + +//-------------------------------------------------------------------------------------------------- + +void AliCentralityGlauberFit::SetOutputFile(TString filename) { + outrootfilename = filename; +} + +//-------------------------------------------------------------------------------------------------- + +void AliCentralityGlauberFit::SetGlauberParam( + Int_t fNmu, + Float_t fmulow, + Float_t fmuhigh, + Int_t fNk, + Float_t fklow, + Float_t fkhigh, + Int_t fNalpha, + Float_t falphalow, + Float_t falphahigh) + +{ + Nmu=fNmu; + mulow=fmulow; + muhigh=fmuhigh; + Nk=fNk; + klow=fklow; + khigh=fkhigh; + Nalpha=fNalpha; + alphalow=falphalow; + alphahigh=falphahigh; +} + +//-------------------------------------------------------------------------------------------------- + +Float_t AliCentralityGlauberFit::GetPercentileCrossSection() { + return percentXsec; +} + +//-------------------------------------------------------------------------------------------------- + +void AliCentralityGlauberFit::MakeFits(TString infilename) { + TH1D *hDATA; + TH1D *thistGlau; + + float efflow = 0.60; + float effhigh = 1.00; + + TFile *outrootfile; + + // open inrootfile, outrootfile + inrootfile = new TFile(infilename); + outrootfile = new TFile(outrootfilename,"RECREATE"); + + // loop over all distribution names + vector::const_iterator hni; + for(hni=histnames.begin(); hni!=histnames.end(); hni++) { + hDATA = (TH1D*) (inrootfile->Get(*hni)); + hDATA->Rebin(10); + TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1)); + int chi2min = 9999; + Float_t alpha_min=-1; + Float_t mu_min=-1; + Float_t k_min=-1; + Float_t eff_min=-1; + Float_t alpha, mu, k, eff, chi2; + + for (int nalpha=0;nalphaClone("hGLAU"); + hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU_%d_%d_%d_%d",(int)(100*mu_min),(int)(100*k_min),(int)(100*alpha_min),(int)(100*eff_min)))); + hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.2f_%.2f_%.2f_%.2f",mu_min,k_min,alpha_min,eff_min))); + + heffi = GetTriggerEfficiencyFunction(hDATA, hGLAU); + SaveHisto(hDATA,hGLAU,heffi,outrootfile); + + } + // close inrootfile, outrootfile + inrootfile->Close(); + outrootfile->Close(); + +} + +//-------------------------------------------------------------------------------------------------- + +TH1D * AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, TH1D *hDATA, Bool_t save) { + + TH1D *hSample = NBDhist(mu,k); + TH1D *h1 = new TH1D(Form("fit_%.2f_%.2f_%.2f_%.2f",mu,k,eff,alpha),"",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1)); + TFile *outFile = NULL; + TNtuple *ntuple = NULL; + + if (save) { + outFile = TFile::Open("pippo.root", "RECREATE"); + ntuple = new TNtuple("gnt", "Glauber ntuple", "npart:ncoll:b:ntot"); + } + + Int_t nents = glau_ntuple->GetEntries(); + for (Int_t i=0;iGetEntry(i); + Int_t n = TMath::Nint(TMath::Power(npart,alpha)); + //Int_t n = TMath::Nint(TMath::Power(ncoll,alpha)); + Int_t ntot=0; + for(Int_t j = 0; jGetRandom(); + h1->Fill(ntot); + + if (save) + ntuple->Fill(npart,ncoll,b,ntot); + + } + + if (save) { + ntuple->Write(); + outFile->Close(); + } + + + return h1; + +} + +//-------------------------------------------------------------------------------------------------- + +Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) { + // note that for different values of neff the mc distribution is identical, just re-scaled below + // normalize the two histogram + // scale MC up but leave data alone - that preserves error bars !!!! + + float mcintegral = thistGlau->Integral(); + for (int i=1; i<=thistGlau->GetNbinsX(); i++) { + float scale = (hDATA->Integral()/mcintegral) / ((float) eff); + float temp = scale * thistGlau->GetBinContent(i); + thistGlau->SetBinContent(i,temp); + } + + // calculate the chi2 between MC and real data over some range ???? + int lowchibin = 10; + int highchibin = 80; + float chi2 = 0.0; + for (int i=lowchibin; i<=highchibin; i++) { + if (hDATA->GetBinContent(i) < 1.0) continue; + float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2); + diff = diff / pow(hDATA->GetBinError(i) , 2); + chi2 += diff; + } + chi2 = chi2 / (highchibin - lowchibin + 1); + return chi2; +} + +//-------------------------------------------------------------------------------------------------- + +TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) { + heffi = (TH1D*)hist1->Clone("heffi"); + heffi->Divide(hist2); + return heffi; +} + +//-------------------------------------------------------------------------------------------------- + +Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) { + effi = hist1->Integral()/hist2->Integral(); + return effi; +} + +//-------------------------------------------------------------------------------------------------- + +void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) { + outrootfile->cd(); + hist1->Write(); + hist2->Write(); + heffi->Write(); +} + +//-------------------------------------------------------------------------------------------------- + +Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) +{ + Double_t mudk = mu/k; + Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) * + TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k); + return ret; +} + +//-------------------------------------------------------------------------------------------------- + +TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k) +{ + TH1D *h = new TH1D("htemp","",100,0,100); + h->SetName(Form("nbd_%f_%f",mu,k)); + h->SetDirectory(0); + for (Int_t i=0;i<100;++i) { + Double_t val = NBD(i,mu,k); + h->Fill(i,val); + } + return h; +} + + diff --git a/PWG2/EVCHAR/AliCentralityGlauberFit.h b/PWG2/EVCHAR/AliCentralityGlauberFit.h new file mode 100644 index 00000000000..0e4aeb17dc2 --- /dev/null +++ b/PWG2/EVCHAR/AliCentralityGlauberFit.h @@ -0,0 +1,68 @@ +#ifndef ALICENTRALITYGLAUBERFIT_H +#define ALICENTRALITYGLAUBERFIT_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// class to determine centrality percentiles from 1D distributions // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include "TObject.h" + +class AliCentralityGlauberFit : public TObject { + + public: + + AliCentralityGlauberFit(); + virtual ~AliCentralityGlauberFit(); + + void SetOutputFile(TString outrootfilename); + void AddHisto(TString name); + void MakeFits(TString infilename); + Float_t GetPercentileCrossSection(); + void SetGlauberParam(Int_t fNmu, Float_t fmulow, Float_t fmuhigh, Int_t fNk, Float_t fklow, Float_t fkhigh, Int_t fNalpha, Float_t falphalow, Float_t falphahigh); + + private: + + Int_t Nmu; + Float_t mulow; + Float_t muhigh; + Int_t Nk; + Float_t klow; + Float_t khigh; + Int_t Nalpha; + Float_t alphalow; + Float_t alphahigh; + + TNtuple *glau_ntuple; + Float_t npart; + Float_t ncoll; + Float_t b; + + Float_t effi; + TH1D *heffi; + + TFile *inrootfile; + TString outrootfilename; + vector histnames; + Float_t percentXsec; + + TH1D * NormalizeHisto(TString hdistributionName); + TH1D * GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, TH1D *hDATA, Bool_t save=kFALSE); + Float_t CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff); + TH1D * GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2); + Float_t GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2); + + Double_t NBD(Int_t n, Double_t mu, Double_t k); + TH1D *NBDhist(Double_t mu, Double_t k); + + void SaveHisto(TH1D *hist1,TH1D *hist2,TH1D *heffi, TFile *outrootfile); + ClassDef(AliCentralityGlauberFit, 1) +}; +#endif + + diff --git a/PWG2/EVCHAR/runAliCentralityBy1D.C b/PWG2/EVCHAR/runAliCentralityBy1D.C new file mode 100644 index 00000000000..1a5395a78ee --- /dev/null +++ b/PWG2/EVCHAR/runAliCentralityBy1D.C @@ -0,0 +1,24 @@ +{ + //load libraries + //gSystem->SetBuildDir("/tmp"); + gROOT->LoadMacro("AliCentralityBy1D.cxx+"); + + const char *finname ="/home/alberica/analysis/zdc/out5/output_f8TOT.root"; // name input file + const char *foutname="test_AliCentralityBy1D_95.root"; // name output file + + float percentXsec=95.0; + + AliCentralityBy1D *mPM = new AliCentralityBy1D(); + mPM.AddHisto("hNtracklets"); + mPM.AddHisto("hNclusters0"); + mPM.AddHisto("hNtracks"); + mPM.AddHisto("hmultV0"); + mPM.AddHisto("hmultFMD"); + + mPM.SetPercentileCrossSection(percentXsec); + mPM.SetPercentileFile(foutname); + + mPM.MakePercentiles(finname); + + +} diff --git a/PWG2/EVCHAR/runAliCentralityByFunction.C b/PWG2/EVCHAR/runAliCentralityByFunction.C new file mode 100644 index 00000000000..01bf9053a64 --- /dev/null +++ b/PWG2/EVCHAR/runAliCentralityByFunction.C @@ -0,0 +1,44 @@ +{ + //load libraries + //gSystem->SetBuildDir("/tmp"); + gROOT->LoadMacro("AliCentralityByFunction.cxx+"); + + const char *finname ="/home/alberica/analysis/zdc/out5/output_f8TOT_correlations.root"; // name input file + const char *foutname="test_AliCentralityByFunction.root"; // name output file + + float percentXsec=100.0; + + AliCentralityByFunction *mPM = new AliCentralityByFunction(); + + mPM.AddHisto("hmultV0vsmultFMD_all"); + mPM.SetFitFunction("hmultV0vsmultFMD_all","fitf_pol3", 0,30000); + + mPM.AddHisto("hNtrackletsvsmultV0_all"); + mPM.SetFitFunction("hNtrackletsvsmultV0_all","fitf_pol2",0,10000); + + mPM.AddHisto("hEzemvsEzdc_all"); + mPM.SetFitFunction("hEzemvsEzdc_all","fitf_pol6",0,3.1); + + // mPM.AddHisto("hNtracksvsEzdc_all"); + // mPM.AddHisto("hNtrackletsvsEzdc_all"); + // mPM.AddHisto("hNclusters0vsEzdc_all"); + // mPM.AddHisto("hmultV0vsEzdc_all"); + // mPM.AddHisto("hmultFMDvsEzdc_all"); + // mPM.AddHisto("hNtracksvsEzem_all"); + // mPM.AddHisto("hNtrackletsvsEzem_all"); + // mPM.AddHisto("hNclusters0vsEzem_all"); + // mPM.AddHisto("hmultV0vsEzem_all"); + // mPM.AddHisto("hmultFMDvsEzem_all"); + // mPM.AddHisto("hNtracksvsmultV0_all"); + // mPM.AddHisto("hNclusters0vsmultV0_all"); + // mPM.AddHisto("hNtracksvsmultFMD_all"); + // mPM.AddHisto("hNtrackletsvsmultFMD_all"); + // mPM.AddHisto("hNclusters0vsmultFMD_all"); + + mPM.SetPercentileCrossSection(percentXsec); + mPM.SetPercentileFile(foutname); + + mPM.MakePercentiles(finname); + + +} diff --git a/PWG2/EVCHAR/runAliCentralityGlauberFit.C b/PWG2/EVCHAR/runAliCentralityGlauberFit.C new file mode 100644 index 00000000000..f3838e465bc --- /dev/null +++ b/PWG2/EVCHAR/runAliCentralityGlauberFit.C @@ -0,0 +1,25 @@ +{ + //load libraries + //gSystem->SetBuildDir("/tmp"); + gROOT->LoadMacro("AliCentralityGlauberFit.cxx+"); + + const char *finname ="/home/alberica/analysis/zdc/out5/output_f8TOT.root"; // name input file + const char *foutname="test_AliCentralityGlauberFit.root"; // name output file + + float percentXsec=95.0; + + AliCentralityGlauberFit *mPM = new AliCentralityGlauberFit(); + mPM.AddHisto("hmultV0"); + mPM.SetGlauberParam(5,60.0,80.0,2,0.5,1.5,7,0.9,1.2); + // mPM.AddHisto("hNclusters0"); + // mPM.AddHisto("hNtracks"); + // mPM.AddHisto("hmultV0"); + // mPM.AddHisto("hmultFMD"); + + // mPM.SetPercentileCrossSection(percentXsec); + mPM.SetOutputFile(foutname); + + mPM.MakeFits(finname); + + +} diff --git a/PWG2/EVCHAR/runCentrality.C b/PWG2/EVCHAR/runCentrality.C index 887c4abb065..0c2cc4df7dc 100644 --- a/PWG2/EVCHAR/runCentrality.C +++ b/PWG2/EVCHAR/runCentrality.C @@ -1,4 +1,4 @@ -void runCentrality(const char * type = "a") +void runCentrality(const char * type = "a", const char *mode="grid") { // Load common libraries gSystem->Load("libCore.so"); @@ -13,16 +13,15 @@ void runCentrality(const char * type = "a") gSystem->Load("libANALYSISalice"); // Use AliRoot includes to compile our task gROOT->ProcessLine(".include $ALICE_ROOT/include"); - // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG2/FORWARD/analysis")); // form filename TString filenameStr = Form("%s",type); - filenameStr = TString("LHC10f8")+filenameStr; + filenameStr = TString("LHC10g2")+filenameStr; const char * filename = filenameStr.Data(); // Create and configure the alien handler plugin - gROOT->LoadMacro("CreateAlienHandler.C"); - AliAnalysisGrid *alienHandler = CreateAlienHandler(filename); + gROOT->LoadMacro("CreateAlienHandlerAOD.C"); + AliAnalysisGrid *alienHandler = CreateAlienHandlerAOD(filename); if(!alienHandler) return; // Create the analysis manager @@ -32,8 +31,8 @@ void runCentrality(const char * type = "a") mgr->SetGridHandler(alienHandler); // My task - gROOT->LoadMacro("AliAnalysisTaskCentralityTreeMaker.cxx++g"); - AliAnalysisTaskCentralityTreeMaker *task = new AliAnalysisTaskCentralityTreeMaker("CentralityTask"); + gROOT->LoadMacro("AliAnalysisTaskCentrality.cxx++g"); + AliAnalysisTaskCentrality *task = new AliAnalysisTaskCentrality("CentralityTask"); // Writing (or not) output tree //task->SetTreeFilling(writeTree); task->SetMCInput(); @@ -59,11 +58,11 @@ void runCentrality(const char * type = "a") "cenHistos.root"); mgr->ConnectOutput(task, 1, coutput1); - AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("coutput2",TTree::Class(), - AliAnalysisManager::kOutputContainer, - "cenTree.root"); - coutput2->SetSpecialOutput(); - mgr->ConnectOutput(task, 2, coutput2); +// AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("coutput2",TTree::Class(), +// AliAnalysisManager::kOutputContainer, +// "cenTree.root"); +// // coutput2->SetSpecialOutput(); +// mgr->ConnectOutput(task, 2, coutput2); // Enable debug printouts mgr->SetDebugLevel(2); @@ -72,6 +71,6 @@ void runCentrality(const char * type = "a") mgr->PrintStatus(); // Start analysis in grid. - mgr->StartAnalysis("grid"); + mgr->StartAnalysis(mode); }; -- 2.43.0