--- /dev/null
+/**************************************************************************
+ * 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 <TTree.h>
+#include <TList.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TCanvas.h>
+
+#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<AliESDEvent*>(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<AliAODEvent*> (InputEvent());
+ // to be implemented
+ printf(" AOD analysis not yet implemented!!!\n\n");
+ return;
+ }
+ PostData(1, fOutput);
+}
+
+
+
+//________________________________________________________________________
+void AliAnalysisTaskCentrality::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+}
+
+
--- /dev/null
+#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
+
--- /dev/null
+/**************************************************************************
+ * 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 <TNamed.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TROOT.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TGraphErrors.h>
+#include <vector>
+#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<TString>::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();
+}
+
+
+
--- /dev/null
+#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<TString> histnames;
+ Float_t percentXsec;
+
+ TH1D * MakePercentHisto(TString hdistributionName);
+ void SaveHisto(TH1D *hist, TFile *outrootfile);
+
+ ClassDef(AliCentralityBy1D, 1)
+};
+#endif
+
+
--- /dev/null
+/**************************************************************************
+ * 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 <TNamed.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TROOT.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TGraphErrors.h>
+#include <vector>
+#include <TMap.h>
+#include <map>
+#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<TString>::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;
+}
+
+
+
--- /dev/null
+#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<TString> histnames;
+ Float_t percentXsec;
+ map<TString, TString>fitfunc; // mapping from distribution to fit function name
+ map<TString, TF1 *>fitter; // mapping from fit function name to corresponding TF1
+
+ TH1D *FitHisto(TString hdistributionName);
+ TH1D * MakePercentHisto(TH2D *hist);
+
+ ClassDef(AliCentralityByFunction, 1)
+};
+#endif
+
+
--- /dev/null
+/**************************************************************************
+ * 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 <TNamed.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TROOT.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TNtuple.h>
+#include <TStyle.h>
+#include <TGraphErrors.h>
+#include <vector>
+#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<TString>::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;nalpha<Nalpha;nalpha++) {
+ alpha = alphalow + ((float) nalpha) * 0.05;
+ mu=0.0;
+ for (int nmu=0; nmu<Nmu; nmu++) {
+ // mu = (mulow*(1-nalpha*0.1*mulow) ) + ((float) nmu ) * (muhigh - mulow ) / Nmu;
+ mu = mulow + ((float) nmu ) * (muhigh - mulow ) / Nmu;
+
+ for (int nk=0; nk<Nk; nk++) {
+ k = klow + ((float) nk ) * (khigh - klow ) / Nk;
+
+ for (int neff=0; neff<20; neff++) {
+ eff = efflow + ((float) neff) * (effhigh - efflow) / 20.0;
+
+ thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
+ chi2 = CalculateChi2(hDATA,thistGlau,eff);
+ if ( chi2 < chi2min ) {
+ chi2min = chi2;
+ alpha_min=alpha;
+ mu_min=mu;
+ k_min=k;
+ eff_min=eff;
+ }
+ }
+ }
+ }
+ }
+
+ thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
+ hGLAU = (TH1D *) thistGlau->Clone("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;i<nents;++i) {
+ glau_ntuple->GetEntry(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; j<n; ++j)
+ ntot+=hSample->GetRandom();
+ 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;
+}
+
+
--- /dev/null
+#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<TString> 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
+
+
--- /dev/null
+{
+ //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);
+
+
+}
--- /dev/null
+{
+ //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);
+
+
+}
--- /dev/null
+{
+ //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);
+
+
+}
-void runCentrality(const char * type = "a")
+void runCentrality(const char * type = "a", const char *mode="grid")
{
// Load common libraries
gSystem->Load("libCore.so");
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
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();
"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);
mgr->PrintStatus();
// Start analysis in grid.
- mgr->StartAnalysis("grid");
+ mgr->StartAnalysis(mode);
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff