From: misko Date: Thu, 16 Oct 2008 13:10:25 +0000 (+0000) Subject: ALICE version of the Universal Correlation analysis package UNICOR. X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=7148817a63b3fbcd8369983d94582430e5c99455;p=u%2Fmrichter%2FAliRoot.git ALICE version of the Universal Correlation analysis package UNICOR. --- diff --git a/UNICOR/AliDAnal.cxx b/UNICOR/AliDAnal.cxx index 8de31f73526..cc79ae1d9ae 100644 --- a/UNICOR/AliDAnal.cxx +++ b/UNICOR/AliDAnal.cxx @@ -1,8 +1,10 @@ // Author: Dariusz Miskowiec 2007 -/////////////////////////////////////////////////////////////////////////////// -// Parent analysis class -/////////////////////////////////////////////////////////////////////////////// +//============================================================================= +// parent class of all analyzers +// keeps the obj array of histograms filled by the daughter +// takes care of storing them on file at the end +//============================================================================= #include #include @@ -10,10 +12,10 @@ ClassImp(AliDAnal) -TDatabasePDG *AliDAnal::fPDG = TDatabasePDG::Instance(); +TDatabasePDG AliDAnal::fgPDG; -/*****************************************************************************/ -AliDAnal::AliDAnal(char *nam) : TNamed(nam,nam), fHistos() +//============================================================================= +AliDAnal::AliDAnal(char *nam) : TNamed(nam,nam), fHistos() { // constructor @@ -23,7 +25,7 @@ AliDAnal::AliDAnal(char *nam) : TNamed(nam,nam), fHistos() printf("%s object named %s created\n",ClassName(),GetName()); } -/*****************************************************************************/ +//============================================================================= void AliDAnal::Save(const char *outfil, const char *mode) { // store histograms on file in a directory named after the object @@ -37,4 +39,4 @@ void AliDAnal::Save(const char *outfil, const char *mode) gROOT->cd(); f->Close(); } -/******************************************************************************/ +//============================================================================= diff --git a/UNICOR/AliDAnal.h b/UNICOR/AliDAnal.h index ee3ac959a07..73068d95a83 100644 --- a/UNICOR/AliDAnal.h +++ b/UNICOR/AliDAnal.h @@ -3,27 +3,23 @@ #ifndef ALIDANAL_H #define ALIDANAL_H -/////////////////////////////////////////////////////////////////////////////// -// Parent analysis class -/////////////////////////////////////////////////////////////////////////////// - #include #include #include -/*****************************************************************************/ +//============================================================================= class AliDAnal : public TNamed { public: - AliDAnal(char *nam); // constructor - virtual ~AliDAnal() {printf("AliDAnal object %s destroyed\n",GetName());} - void Save(const char *outfil, const char *mode="update"); // save histograms on root file + AliDAnal(char *nam); // constructor + virtual ~AliDAnal() {printf("%s object named %s deleted\n",ClassName(),GetName());} + void Save(const char *outfil, const char *mode="update"); // save histograms protected: - static TDatabasePDG *fPDG; // particle database - TObjArray fHistos; // histograms + static TDatabasePDG fgPDG; // particle database + TObjArray fHistos; // histograms ClassDef(AliDAnal,1) }; -/*****************************************************************************/ +//============================================================================= #endif diff --git a/UNICOR/AliDAnalCorrel.cxx b/UNICOR/AliDAnalCorrel.cxx new file mode 100644 index 00000000000..8b3e91c5289 --- /dev/null +++ b/UNICOR/AliDAnalCorrel.cxx @@ -0,0 +1,97 @@ +// Author: Dariusz Miskowiec 2005 + +//============================================================================= +// two-particle correlation analyzer +//============================================================================= + +#include +#include +#include "AliDEvent.h" +#include "AliDHN.h" +#include "AliDAnalCorrel.h" + +ClassImp(AliDAnalCorrel) + +//============================================================================= +AliDAnalCorrel::AliDAnalCorrel(Char_t *nam, Double_t emi, Double_t ema, + Int_t pid0, Int_t pid1): + AliDAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fPa() +{ + // constructor + // emi and ema define the rapidity range for histogram + + TParticlePDG *part0 = AliDAnal::fgPDG.GetParticle(fPid0); + TParticlePDG *part1 = AliDAnal::fgPDG.GetParticle(fPid1); + fMass0 = part0? part0->Mass() : 0; + fMass1 = part1? part1->Mass() : 0; + + double pi = TMath::Pi(); + TAxis *ax[8]; + ax[0] = new TAxis(3,-0.5,2.5);ax[0]->SetTitle("trumixrot"); + ax[1] = new TAxis(5,0,0.5); ax[1]->SetTitle("centrality"); + ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("pair y"); + ax[3] = new TAxis(8,-pi,pi); ax[3]->SetTitle("pair phi"); // wrt event plane + double a0[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0}; + ax[4] = new TAxis(7,a0); ax[4]->SetTitle("(pair pt)/2 (GeV)"); + ax[5] = new TAxis(8,0,pi); ax[5]->SetTitle("q-theta"); + ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi"); + double a1[100]; + for (int i=0;i<20;i++) a1[i]=i*0.005; + for (int i=0;i<45;i++) a1[20+i]=0.1+i*0.02; + ax[7] = new TAxis(64,a1); ax[7]->SetTitle("q (GeV/c)"); + AliDHN *pair = new AliDHN("pair",8,ax); + for (int i=0; i<8; i++) delete ax[i]; + fHistos.Add(pair); + gROOT->cd(); + printf("%s object named %s created\n",ClassName(),GetName()); +} +//============================================================================= +void AliDAnalCorrel::Process(Int_t tmr, AliDEvent *ev0, AliDEvent *ev1, Double_t phirot) +{ + // process pairs from one or two (if mixing) events + // tmr tells which histogram (bins) to fill: tru,mix,rot + + static TRandom2 ran; + AliDHN *pair = (AliDHN*) fHistos.At(0); + + // mixing-and-rotating-proof centrality and reaction plane angle + // (but not rotation-proof for rotation angles much different from 0 and 180) + // true and rotated pairs are within the triangle (jCentrality()+ev1->Centrality())/2.0; + double q0x,q0y,q1x,q1y; + ev0->RP(q0x,q0y); + ev1->RP(q1x,q1y); + double rpphi = atan2(q0y+q1y,q0x+q1x); + + // loop over pairs + + for (int i=0; iNParticles(); i++) { + if (!ev0->ParticleGood(i,fPid0)) continue; + for (int j=0; jNParticles(); j++) { + if (ev0==ev1 && jParticleGood(j,fPid1)) continue; + if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i), + ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot)) continue; + fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i)); + fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot); + if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap(); + fPa.CalcLAB(); + fPa.CalcPairCM(); + if (fPa.QCM()==0) return; // should not be too frequent + double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi); + pair->Fill(tmr, // 0 for tru, 1 for mix, 2 for rot + cent, // centrality + fPa.Rapidity(), // pair rapidity + phi, // pair phi wrt reaction plane + fPa.Pt()/2.0, // half of pair pt + fPa.QCMTheta(), // polar angle of Q + fPa.QCMPhiOut(), // azimuthal angle of Q w.r.t. out + fPa.QCM(), // |p2-p1| in c.m.s. + 1.0); // weigth + } + } +} +//============================================================================= diff --git a/UNICOR/AliDAnalCorrel.h b/UNICOR/AliDAnalCorrel.h new file mode 100644 index 00000000000..b93b6b4136c --- /dev/null +++ b/UNICOR/AliDAnalCorrel.h @@ -0,0 +1,32 @@ +// Author: Dariusz Miskowiec 2005 + +#ifndef ALIDANALCORREL_H +#define ALIDANALCORREL_H +#include "AliDAnal.h" +#include "AliDPair.h" +class TH1D; +class TH2D; +class AliDEvent; + +//============================================================================= +class AliDAnalCorrel : public AliDAnal { + + public: + AliDAnalCorrel(Char_t *nam="correl", + Double_t emi=-1, Double_t ema=1, + Int_t pid0=0, Int_t pid1=0); // constructor + virtual ~AliDAnalCorrel(){} // destructor + // process one (tru) or two (mix) events + void Process(Int_t tmr, AliDEvent *ev0, AliDEvent *ev1, Double_t phirot); + + protected: + Int_t fPid0; // particle species 0 + Int_t fPid1; // particle species 1 + Double_t fMass0; // mass 0 + Double_t fMass1; // mass 1 + AliDPair fPa; // pair buffer for calculations + + ClassDef(AliDAnalCorrel,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/AliDAnalGlobal.cxx b/UNICOR/AliDAnalGlobal.cxx new file mode 100644 index 00000000000..b3780d86244 --- /dev/null +++ b/UNICOR/AliDAnalGlobal.cxx @@ -0,0 +1,53 @@ +// Author: Dariusz Miskowiec 2007 + +//============================================================================= +// event global variable analyzer +//============================================================================= + +#include +#include +#include +#include "AliDEvent.h" +#include "AliDAnalGlobal.h" + +ClassImp(AliDAnalGlobal) + +//============================================================================= +AliDAnalGlobal::AliDAnalGlobal(Char_t *nam) : AliDAnal(nam) +{ + // constructor + + TH1D *mult = new TH1D("mult","mult",5000,-0.5,4999.5); + mult->SetXTitle("multiplicity"); + TH1D *cent = new TH1D("cent","cent",100,0,1); + cent->SetXTitle("centrality"); + TH2D *dire = new TH2D("dire","dire",100,-40,40,100,-40,40); + dire->SetXTitle("Qx (GeV)"); + dire->SetYTitle("Qy (GeV)"); + TH1D *zver = new TH1D("zver","zver",120,-1.2,1.2); + zver->SetXTitle("normalized z-vertex"); + fHistos.Add(mult); + fHistos.Add(cent); + fHistos.Add(dire); + fHistos.Add(zver); + gROOT->cd(); + printf("%s object named %s created\n",ClassName(),GetName()); +} +//============================================================================= +void AliDAnalGlobal::Process(AliDEvent *ev) +{ + // fill event variable histograms + + TH1D *mult = (TH1D*) fHistos.At(0); + TH1D *cent = (TH1D*) fHistos.At(1); + TH2D *dire = (TH2D*) fHistos.At(2); + TH1D *zver = (TH1D*) fHistos.At(3); + + mult->Fill(ev->NParticles(),1.0); + cent->Fill(ev->Centrality(),1.0); + Double_t qx,qy; + ev->RP(qx,qy); + dire->Fill(qx,qy,1.0); + zver->Fill(ev->Zver(),1.0); +} +//============================================================================= diff --git a/UNICOR/AliDAnalGlobal.h b/UNICOR/AliDAnalGlobal.h new file mode 100644 index 00000000000..42651b42a6f --- /dev/null +++ b/UNICOR/AliDAnalGlobal.h @@ -0,0 +1,20 @@ +// Author: Dariusz Miskowiec 2007 + +#ifndef ALIDANALGLOBAL_H +#define ALIDANALGLOBAL_H + +#include "AliDAnal.h" +class AliDEvent; + +//============================================================================= +class AliDAnalGlobal : public AliDAnal { + + public: + AliDAnalGlobal(Char_t *nam="global"); // constructor + virtual ~AliDAnalGlobal(){} // destructor + void Process(AliDEvent *ev); // fill histograms + + ClassDef(AliDAnalGlobal,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/AliDAnalPtfluc.cxx b/UNICOR/AliDAnalPtfluc.cxx new file mode 100644 index 00000000000..635310ccbf9 --- /dev/null +++ b/UNICOR/AliDAnalPtfluc.cxx @@ -0,0 +1,94 @@ +// Author: Dariusz Miskowiec 2008 + +//============================================================================= +// pt-fluctuations analyzer +//============================================================================= + +#include +#include +#include +#include +#include "AliDEvent.h" +#include "AliDHN.h" +#include "AliDAnalPtfluc.h" + +ClassImp(AliDAnalPtfluc) + +//============================================================================= +AliDAnalPtfluc::AliDAnalPtfluc(Char_t *nam, Int_t pid0, Int_t pid1) : + AliDAnal(nam), fPid0(pid0), fPid1(pid1) +{ + // constructor + + TAxis *ax[5]; + ax[0] = new TAxis(2,-0.5,1.5); ax[0]->SetTitle("trumix"); + ax[1] = new TAxis(9,0,0.9); ax[1]->SetTitle("centrality"); + ax[2] = new TAxis(6,-0.5,5.5); ax[2]->SetTitle("n-pt0-pt1-pt00-pt11-pt01"); + ax[3] = new TAxis(48,-180,180); ax[3]->SetTitle("dphi (deg)"); + ax[4] = new TAxis(40,-2,2); ax[4]->SetTitle("deta"); + AliDHN *pair = new AliDHN("pair",5,ax); + for (int i=0; i<5; i++) delete ax[i]; + fHistos.Add(pair); + gROOT->cd(); + printf("%s object named %s created\n",ClassName(),GetName()); +} +//============================================================================= +void AliDAnalPtfluc::Process(Int_t tmr, AliDEvent *ev0, AliDEvent *ev1) +{ + // process pairs from one or two (if mixing) events + + double ptmin=0.1; // GeV + double ptmax=1.5; // GeV + double etamin=-9; + double etamax=9; + + // mixing-and-rotating-proof centrality + + double cent = (ev0->Centrality()+ev1->Centrality())/2.0; + + // loop over pairs + + AliDHN *pair = (AliDHN*) fHistos.At(0); + static TRandom2 ran; + for (int i=0; iNParticles(); i++) { + if (!ev0->ParticleGood(i,fPid0)) continue; + double eta0 = ev0->ParticleEta(i); + double phi0 = ev0->ParticlePhi(i); + double pt0 = ev0->ParticlePt(i); + if (eta0 < etamin) continue; + if (eta0 > etamax) continue; + if (pt0 < ptmin) continue; + if (pt0 > ptmax) continue; + for (int j=0; jNParticles(); j++) { + if (ev0==ev1 && j==i) continue; + if (ev0==ev1 && jParticleGood(j,fPid1)) continue; + double eta1 = ev1->ParticleEta(j); + double phi1 = ev1->ParticlePhi(j); + double pt1 = ev1->ParticlePt(j); + if (eta1 < etamin) continue; + if (eta1 > etamax) continue; + if (pt1 < ptmin) continue; + if (pt1 > ptmax) continue; + double deta = eta1-eta0; + double dphi = phi1-phi0; + // randomize order + if (ran.Rndm()<0.5) { + double buf = pt0; + pt0 = pt1; + pt1 = buf; + deta = -deta; + dphi = -dphi; + } + dphi = TVector2::Phi_mpi_pi(dphi); + dphi*=TMath::RadToDeg(); + pair->Fill(tmr, cent, 0, dphi, deta, 1); // number of pairs + pair->Fill(tmr, cent, 1, dphi, deta, pt0); + pair->Fill(tmr, cent, 2, dphi, deta, pt1); + pair->Fill(tmr, cent, 3, dphi, deta, pt0*pt0); + pair->Fill(tmr, cent, 4, dphi, deta, pt1*pt1); + pair->Fill(tmr, cent, 5, dphi, deta, pt0*pt1); + } + } +} +//============================================================================= diff --git a/UNICOR/AliDAnalPtfluc.h b/UNICOR/AliDAnalPtfluc.h new file mode 100644 index 00000000000..7fc275cedc0 --- /dev/null +++ b/UNICOR/AliDAnalPtfluc.h @@ -0,0 +1,24 @@ +// Author: Dariusz Miskowiec 2005 + +#ifndef ALIDANALPTFLUC_H +#define ALIDANALPTFLUC_H + +#include "AliDAnal.h" +class AliDEvent; +class AliDHN; + +//============================================================================= +class AliDAnalPtfluc : public AliDAnal { + + public: + AliDAnalPtfluc(Char_t *nam="correl", Int_t pid0=0, Int_t pid1=0); // constructor + virtual ~AliDAnalPtfluc(){} // destructor + void Process(Int_t tmr, AliDEvent *ev0, AliDEvent *ev1); // process event(s) + + protected: + Int_t fPid0; // particle species 0 + Int_t fPid1; // particle species 1 + ClassDef(AliDAnalPtfluc,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/AliDAnalSingle.cxx b/UNICOR/AliDAnalSingle.cxx new file mode 100644 index 00000000000..1dea57a16cd --- /dev/null +++ b/UNICOR/AliDAnalSingle.cxx @@ -0,0 +1,70 @@ +// Author: Dariusz Miskowiec 2007 + +//============================================================================= +// single particle analyzer +//============================================================================= + +#include +#include +#include +#include +#include "AliDHN.h" +#include "AliDEvent.h" +#include "AliDAnalSingle.h" + +ClassImp(AliDAnalSingle) + +//============================================================================= +AliDAnalSingle::AliDAnalSingle(Char_t *nam, Double_t emi, Double_t ema, Int_t pid) : + AliDAnal(nam), fPid(pid), fMass(0.0) +{ + // constructor + // emi and ema define the rapidity range for histograms + + fPid = pid; + TParticlePDG *part = AliDAnal::fgPDG.GetParticle(fPid); + fMass = part? part->Mass() : 0; + + double pi = TMath::Pi(); + TAxis *ax[10]; + ax[0] = new TAxis(30,-1,1); ax[0]->SetTitle("vertex z"); + ax[1] = new TAxis(80,emi,ema); ax[1]->SetTitle("eta"); + ax[2] = new TAxis(90,-pi,pi); ax[2]->SetTitle("phi"); + AliDHN *zep = new AliDHN("zep",3,ax); + for (int i=0; i<3; i++) delete ax[i]; + + ax[0] = new TAxis(20,0,1); ax[0]->SetTitle("centrality"); + ax[1] = new TAxis(80,emi,ema); ax[1]->SetTitle("y"); + ax[2] = new TAxis(80,0,2); ax[2]->SetTitle("pt (GeV)"); + AliDHN *cyp = new AliDHN("cyp",3,ax); + for (int i=0; i<3; i++) delete ax[i]; + + ax[0] = new TAxis(10,emi,ema); ax[0]->SetTitle("eta"); + ax[1] = new TAxis(150,0,3); ax[1]->SetTitle("p (GeV)"); + ax[2] = new TAxis(150,0.5,3.5);ax[2]->SetTitle("sqrt(dedx (mips))"); + AliDHN *epd = new AliDHN("epd",3,ax); + for (int i=0; i<3; i++) delete ax[i]; + + fHistos.Add(zep); + fHistos.Add(cyp); + fHistos.Add(epd); + gROOT->cd(); + printf("%s object named %s created\n",ClassName(),GetName()); +} +//============================================================================= +void AliDAnalSingle::Process(AliDEvent *ev) +{ + // fill single particle histograms + + AliDHN *zep = (AliDHN*) fHistos.At(0); + AliDHN *cyp = (AliDHN*) fHistos.At(1); + AliDHN *epd = (AliDHN*) fHistos.At(2); + for (int i=0; iNParticles(); i++) { + if (!ev->ParticleGood(i,fPid)) continue; + zep->Fill(ev->Zver(),ev->ParticleEta(i),ev->ParticlePhi(i),1.0); + double y = fMass>0? ev->ParticleY(i,fMass) : ev->ParticleEta(i); + cyp->Fill(ev->Centrality(),y,ev->ParticlePt(i),1.0); + epd->Fill(ev->ParticleEta(i),ev->ParticleP(i),TMath::Sqrt(ev->ParticleDedx(i)),1.0); + } +} +//============================================================================= diff --git a/UNICOR/AliDAnalSingle.h b/UNICOR/AliDAnalSingle.h new file mode 100644 index 00000000000..a540bf0391e --- /dev/null +++ b/UNICOR/AliDAnalSingle.h @@ -0,0 +1,27 @@ +// Author: Dariusz Miskowiec 2007 + +#ifndef ALIDANALSINGLE_H +#define ALIDANALSINGLE_H + +#include "AliDAnal.h" +class AliDEvent; +class AliDHN; + +//============================================================================= +class AliDAnalSingle : public AliDAnal { + + public: + AliDAnalSingle(Char_t *nam="single", + Double_t emi=-1, Double_t ema=1, + Int_t pid=0); // constructor + virtual ~AliDAnalSingle(){} // destructor + void Process(AliDEvent *ev); // fill histograms + + protected: + Int_t fPid; // pid; 0 means all + Double_t fMass; // mass (if pid!=0) + + ClassDef(AliDAnalSingle,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/AliDAnalysisTask.cxx b/UNICOR/AliDAnalysisTask.cxx new file mode 100644 index 00000000000..da2cd5ab66b --- /dev/null +++ b/UNICOR/AliDAnalysisTask.cxx @@ -0,0 +1,110 @@ +//Author: Dariusz Miskowiec 2007 + +//============================================================================= +// my analysis task +//============================================================================= +#include "TChain.h" +#include "AliESDInputHandler.h" +#include "AliAnalysisManager.h" +#include "AliDAnalysisTask.h" +#include "AliDAnalGlobal.h" +#include "AliDAnalSingle.h" +#include "AliDAnalCorrel.h" +#include "AliDAnalPtfluc.h" +#include "AliDEventAliceESD.h" + +ClassImp(AliDAnalysisTask) + +/*****************************************************************************/ +AliDAnalysisTask::AliDAnalysisTask() : + AliAnalysisTask("dali", ""), + fESD(0), fEv0(0), fEv1(0), + fDag(0), fPim(0), fPip(0), + fCnn(0), fCpp(0), fPtf(0), + fOutputList(0) +{ + // constructor + + fEv0 = new AliDEventAliceESD(); + DefineInput(0, TChain::Class()); + DefineOutput(0, TList::Class()); +} +/*****************************************************************************/ +void AliDAnalysisTask::CreateOutputObjects() +{ + // executed once on each worker + + fDag = new AliDAnalGlobal("dag"); + fAll = new AliDAnalSingle("all",fEv0->Etamin(),fEv0->Etamax(),0); + fPim = new AliDAnalSingle("pim",fEv0->Etamin(),fEv0->Etamax(),-211); + fPip = new AliDAnalSingle("pip",fEv0->Etamin(),fEv0->Etamax(), 211); + fCnn = new AliDAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211); + fCpp = new AliDAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211); + fPtf = new AliDAnalPtfluc("ptf",0,0); + fOutputList = new TList(); + fOutputList->Add(fDag); + fOutputList->Add(fAll); + fOutputList->Add(fPim); + fOutputList->Add(fPip); + fOutputList->Add(fCnn); + fOutputList->Add(fCpp); + fOutputList->Add(fPtf); +} +/*****************************************************************************/ +void AliDAnalysisTask::ConnectInputData(Option_t *) +{ +// connect ESD or AOD here +// called on each input data change. + + TTree* tree = dynamic_cast (GetInputData(0)); + if (!tree) { + Printf("ERROR: Could not read chain from input slot 0"); + } else { + fEv0->AttachTree(tree); + } + AliESDInputHandler *esdH = dynamic_cast + (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + + if (!esdH) Printf("ERROR: Could not get ESDInputHandler"); + else fESD = esdH->GetEvent(); +} +/*****************************************************************************/ +void AliDAnalysisTask::Exec(Option_t */*option*/) +{ + // process one event + + + if (!fESD) { + Printf("ERROR: fESD not available"); + return; + } + // If AttachTree aka ReadFromTree above cannot be done then this is + // the alternative: shallow copy of the current alice esd event to fEv0 + // memcpy((AliESDEvent*)fEv0, fESD, sizeof(AliESDEvent)); + + if (!fEv0->Good()) return; + fDag->Process(fEv0); + fPim->Process(fEv0); + fAll->Process(fEv0); + fPip->Process(fEv0); + fCnn->Process(0,fEv0,fEv0,0); + fCnn->Process(2,fEv0,fEv0,TMath::DegToRad()*180); + fCpp->Process(0,fEv0,fEv0,0); + fCpp->Process(2,fEv0,fEv0,TMath::DegToRad()*180); + fPtf->Process(0,fEv0,fEv0); + PostData(0, fOutputList); +} +/*****************************************************************************/ +void AliDAnalysisTask::Terminate(Option_t *) +{ + // terminate + printf("terminate\n"); + fDag->Save("unicor-result.root","recreate"); + fAll->Save("unicor-result.root"); + fPim->Save("unicor-result.root"); + fPip->Save("unicor-result.root"); + fCnn->Save("unicor-result.root"); + fCpp->Save("unicor-result.root"); + fPtf->Save("unicor-result.root"); +} +/*****************************************************************************/ diff --git a/UNICOR/AliDAnalysisTask.h b/UNICOR/AliDAnalysisTask.h new file mode 100644 index 00000000000..1adc1d21951 --- /dev/null +++ b/UNICOR/AliDAnalysisTask.h @@ -0,0 +1,48 @@ +// Author: Dariusz Miskowiec 2007 + +#ifndef ALIDANALYSISTASK_H +#define ALIDANALYSISTASK_H + +#include "AliESDEvent.h" +#include "AliAnalysisTask.h" +class AliDEventAliceESD; +class AliDAnalGlobal; +class AliDAnalSingle; +class AliDAnalCorrel; +class AliDAnalPtfluc; + +/*****************************************************************************/ +class AliDAnalysisTask : public AliAnalysisTask { + + public: + AliDAnalysisTask(); // constructor + virtual ~AliDAnalysisTask() {} // destructor + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *option); + virtual void Terminate(Option_t *); + virtual void LocalInit() {} + virtual Bool_t Notify() {return kTRUE;} + virtual Bool_t NotifyBinChange() {return kTRUE;} + virtual void FinishTaskOutput() {} + + protected: + AliESDEvent *fESD; //! ESD event + AliDEventAliceESD *fEv0; //! data/analysis interface + AliDEventAliceESD *fEv1; //! another for event mixing + AliDAnalGlobal *fDag; //! global analysis + AliDAnalSingle *fAll; //! single analysis + AliDAnalSingle *fPim; //! single analysis + AliDAnalSingle *fPip; //! single analysis + AliDAnalCorrel *fCnn; //! correlation analysis pi-pi- + AliDAnalCorrel *fCpp; //! correlation analysis pi+pi+ + AliDAnalPtfluc *fPtf; //! pt-fluctuation analysis + TList *fOutputList; // list of output objects + AliDAnalysisTask(const AliDAnalysisTask&); + AliDAnalysisTask& operator=(const AliDAnalysisTask&); + + ClassDef(AliDAnalysisTask,1) +}; +/*****************************************************************************/ + +#endif diff --git a/UNICOR/AliDEvent.cxx b/UNICOR/AliDEvent.cxx new file mode 100644 index 00000000000..438d1404a5b --- /dev/null +++ b/UNICOR/AliDEvent.cxx @@ -0,0 +1,58 @@ +// Author: Dariusz Miskowiec 2007 + +//============================================================================= +// parent class of all events; analyzers access data via this class +//============================================================================= + +#include +#include "AliDEvent.h" + +ClassImp(AliDEvent) + +//============================================================================= +void AliDEvent::RP(Double_t &qx, Double_t &qy, Int_t harmonic) const +{ + // simplest flow vector + + qx=0; + qy=0; + for (int i=0; i2.0) pt = 2.0; // from 2 GeV flow saturates anyway + qx += pt*cos(harmonic*ParticlePhi(i)); + qy += pt*sin(harmonic*ParticlePhi(i)); + } +} +//============================================================================= +Double_t AliDEvent::ParticleEta(Int_t i) const +{ + // pseudorapidity + + double the = ParticleTheta(i); + if (the<0.0001) return 10; + else if (the>TMath::Pi()-0.0001) return -10; + return -TMath::Log(TMath::Tan(the/2)); +} +//============================================================================= +Double_t AliDEvent::ParticleY(Int_t i, Double_t mass) const +{ + // rapidity + + double pp = ParticleP(i); + double ee = sqrt(fabs(mass*mass + pp*pp)); + double pz = ParticlePz(i); + double yy = log((ee+pz)/(ee-pz))/2; + return yy; +} +//============================================================================= +Double_t AliDEvent::PairDCA(Int_t i0, Int_t i1, + Double_t *x0, Double_t *y0, + Double_t *x1, Double_t *y1) const +{ + // distance of closest approach for pair of tracks; under construction + + return 10; +} +//============================================================================= + diff --git a/UNICOR/AliDEvent.h b/UNICOR/AliDEvent.h new file mode 100644 index 00000000000..d731b473e68 --- /dev/null +++ b/UNICOR/AliDEvent.h @@ -0,0 +1,49 @@ +// Author: Dariusz Miskowiec 2007 + +#ifndef ALIDEVENT_H +#define ALIDEVENT_H + +#include +#include + +class TTree; + +//============================================================================= +class AliDEvent : public TObject { + public: + AliDEvent() : TObject(){printf("%s object named %s created\n",ClassName(),GetName());} + virtual ~AliDEvent() {printf("%s object named %s deleted\n",ClassName(),GetName());} + + // interface part + + virtual void AttachTree(TTree *tr) = 0; + virtual Double_t Etamin() const = 0; // experiment's acceptance + virtual Double_t Etamax() const = 0; + virtual Bool_t Good() const = 0; + virtual Double_t Centrality() = 0; // centrality (0,1); 0 is most central + virtual void RP(Double_t &qx, Double_t &qy) const = 0; + virtual Double_t RPphi() const = 0; + virtual Double_t Zver() const = 0; // z-vertex (-1,1) + virtual Int_t NParticles() const = 0; + + virtual Bool_t ParticleGood(Int_t i, Int_t pidi) const = 0; + virtual Double_t ParticleP(Int_t i) const = 0; + virtual Double_t ParticleTheta(Int_t i) const = 0; + virtual Double_t ParticlePhi(Int_t i) const = 0; + virtual Double_t ParticleDedx(Int_t i) const = 0; + virtual Bool_t PairGood(Double_t p0, Double_t the0, Double_t phi0, + Double_t p1, Double_t the1, Double_t phi1) const = 0; + + // toolkit part + + void RP(Double_t &qx, Double_t &qy, Int_t harmonic) const; + Double_t ParticlePt(Int_t i) const {return ParticleP(i)*sin(ParticleTheta(i));} + Double_t ParticlePz(Int_t i) const {return ParticleP(i)*cos(ParticleTheta(i));} + Double_t ParticleEta(Int_t i) const; + Double_t ParticleY(Int_t i, Double_t mass) const; + Double_t PairDCA(Int_t i0, Int_t i1, Double_t *x0, Double_t *y0, Double_t *x1, Double_t *y1) const; + + ClassDef(AliDEvent,0) +}; +#endif +//============================================================================= diff --git a/UNICOR/AliDEventAliceESD.cxx b/UNICOR/AliDEventAliceESD.cxx new file mode 100644 index 00000000000..816905b7b7a --- /dev/null +++ b/UNICOR/AliDEventAliceESD.cxx @@ -0,0 +1,83 @@ +// Author: Dariusz Miskowiec 2008 + +//============================================================================= +// AliDEvent-AliESD interface // +//============================================================================= + +#include "AliESDEvent.h" +#include "AliDEventAliceESD.h" + +ClassImp(AliDEventAliceESD) + +//============================================================================= +AliDEventAliceESD::AliDEventAliceESD() : AliDEvent(), AliESDEvent() +{ + // constructor + + printf("%s object created\n",AliDEvent::ClassName()); +} +//============================================================================= +AliDEventAliceESD::~AliDEventAliceESD() +{ + // destructor + +} +//============================================================================= +Bool_t AliDEventAliceESD::Good() const +{ + // event cuts + + if (fabs(Zver())>1) return kFALSE; + if (AliESDEvent::GetPrimaryVertex()->GetZRes()>0.1) return kFALSE; + return kTRUE; +} +//============================================================================= +Bool_t AliDEventAliceESD::ParticleGood(Int_t i, Int_t pidi) const +{ + // track cuts and particle id cuts; pidi=0 means take all tracks + // consider using the standard ESDcut + + // track quality cuts + + AliESDtrack *track = AliESDEvent::GetTrack(i); + if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0; // TPC refit + if (track->GetTPCNcls() < 120) return 0; // number of TPC clusters + const AliExternalTrackParam *tp = track->GetTPCInnerParam(); + if (!tp) return 0; + + Float_t r,z; + track->GetImpactParameters(r,z); + // if (fabs(z)>0.2) return 0; // impact parameter in z + // if (fabs(r)>0.1) return 0; // impact parameter in xy + + // pid + + if (pidi==0) return 1; + if (!track->IsOn(AliESDtrack::kTPCpid)) return 0; + Double_t p[AliPID::kSPECIES]; + track->GetESDpid(p); + Int_t q = tp->Charge(); + if (pidi == -211) return p[AliPID::kPion]>0.5 && q==-1; + else if (pidi == 211) return p[AliPID::kPion]>0.5 && q==1; + else return 0; +} +//============================================================================= +Bool_t AliDEventAliceESD::PairGood(Double_t p0, Double_t the0, Double_t phi0, + Double_t p1, Double_t the1, Double_t phi1) const { + + // two-track separation cut + + double r = 85; // TPC entrance radius in cm + double x0 = r*sin(the0)*cos(phi0); + double x1 = r*sin(the1)*cos(phi0); + double y0 = r*sin(the0)*sin(phi0); + double y1 = r*sin(the1)*sin(phi1); + double z0 = r*cos(the0); + double z1 = r*cos(the1); + double dx = x1-x0; + double dy = y1-y0; + double dz = z1-z0; + double dist2 = dx*dx+dy*dy+dz*dz; + return dist2>2*2; +} +//============================================================================= diff --git a/UNICOR/AliDEventAliceESD.h b/UNICOR/AliDEventAliceESD.h new file mode 100644 index 00000000000..ce3c83ff35d --- /dev/null +++ b/UNICOR/AliDEventAliceESD.h @@ -0,0 +1,38 @@ +// Author: Dariusz Miskowiec 2008 + +#ifndef ALIDEVENTALICEESD_H +#define ALIDEVENTALICEESD_H + +#include "TVector2.h" +#include "AliDEvent.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" + +//============================================================================= +class AliDEventAliceESD : public AliDEvent, public AliESDEvent { + + public: + AliDEventAliceESD(); + virtual ~AliDEventAliceESD(); + Double_t Etamin() const {return -0.75;} + Double_t Etamax() const {return 0.75;} + void AttachTree(TTree *tr) {ReadFromTree(tr);} + Bool_t Good() const; + Double_t Centrality() {return 0.9999*exp(-NParticles()/20.0);} // OK for pp + void RP(Double_t &qx, Double_t &qy) const {AliDEvent::RP(qx,qy,2);} + Double_t RPphi() const {Double_t qx,qy; RP(qx,qy); return TMath::ATan2(qy,qx);} + Double_t Zver() const {return AliESDEvent::GetPrimaryVertex()->GetZv()/10.0;} + Int_t NParticles() const {return AliESDEvent::GetNumberOfTracks();} + + Bool_t ParticleGood(Int_t i, Int_t pidi=0) const; + Double_t ParticleP(Int_t i) const {return AliESDEvent::GetTrack(i)->GetTPCInnerParam()->P();} + Double_t ParticleTheta(Int_t i) const {return AliESDEvent::GetTrack(i)->GetTPCInnerParam()->Theta();} + Double_t ParticlePhi(Int_t i) const {return TVector2::Phi_mpi_pi(AliESDEvent::GetTrack(i)->GetTPCInnerParam()->Phi());} + Double_t ParticleDedx(Int_t i) const {return AliESDEvent::GetTrack(i)->GetTPCsignal()/47.0;} + Bool_t PairGood(Double_t p0, Double_t the0, Double_t phi0, + Double_t p1, Double_t the1, Double_t phi1) const; + // alternative: GetTPCInnerParam, GetConstrainedParam + ClassDef(AliDEventAliceESD,0) +}; +#endif +//============================================================================= diff --git a/UNICOR/AliDHN.cxx b/UNICOR/AliDHN.cxx new file mode 100644 index 00000000000..5bd0b04acbc --- /dev/null +++ b/UNICOR/AliDHN.cxx @@ -0,0 +1,366 @@ +// Author: Dariusz Miskowiec 2007 + +//============================================================================= +// multidimensional histogram +// Physically, the data are kept in one single one-dimensional histogram. +// Projecting on 1, 2, and n-1 dimensions is implemented. +// The histogram can be saved on root file as the one-dimensional data +// histogram and the axes, thus eternal forward compatibility is ensured. +//============================================================================= + +#include +#include +#include +#include +#include +#include "AliDHN.h" + +ClassImp(AliDHN) + +//============================================================================= +AliDHN::AliDHN(Char_t *nam, Int_t ndim, TAxis **ax) + : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5), + fNdim(ndim) +{ + // Constructor for building from scratch. + + // Above, we just have managed to create a single one-dimensional histogram + // with the number of bins equal to the product of the numbers of bins in all + // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5. + + for (int i=0; iCopy(fAxis[i]); + + // for speed reasons, number of bins of each axis is stored on fNbins + // and the running product of the latter on fMbins + // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1} + + fMbins[fNdim-1] = 1; + for (int i=0; i0; i--) fMbins[i-1] = fMbins[i]*fNbins[i]; + printf("%d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX()); +} +//============================================================================= +AliDHN::AliDHN(Char_t *filnam, Char_t *nam) + : TH1D(*((TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"))), + fNdim(0) +{ + // Constructor for reading from file. + + TFile *f = TFile::Open(filnam,"read"); + f->cd(nam); + TAxis *ax[fMaxNdim]; + char axnam[1000]; + for (fNdim=0; fNdimGet(axnam); + if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]); + else break; + } + f->Close(); + + fMbins[fNdim-1] = 1; + for (int i=0; i0; i--) fMbins[i-1] = fMbins[i]*fNbins[i]; + + if (GetNbinsX()!=Albins(fNdim,ax)) { + printf("number of bins of histo %d differs from product of nbins of axes %d\n", + GetNbinsX(),Albins(fNdim,ax)); + printf("bombing\n"); + exit(-1); + } + printf("%d-dimensional histogram %s with %d bins read from file %s\n", + fNdim,nam,GetNbinsX(),filnam); +} +//============================================================================= +Int_t AliDHN::Albins(Int_t n, TAxis **ax) +{ + // Calculate product of nbins of ax[0]...ax[n-1] + // (= total number of bins of the multidimensional histogram to be made). + + Int_t nbins = 1; + // while (n--) nbins *= ax[n]->GetNbins(); + for (int i=0; iGetNbins(); + return nbins; +} +//============================================================================= +Int_t AliDHN::MulToOne(Int_t *k) const +{ + // Calculate the 1-dim index n from n-dim indices k[fNdim]. + // Valid k[i] should be between 0 and fNbins[i]-1. + // Valid result will be between 0 and GetNbinsX()-1. + // Return -1 if under- or overflow in any dimension. + + Int_t n = 0; + for (int i=0; i=fNbins[i]) return -1; + n += fMbins[i]*k[i]; + } + return n; +} +//============================================================================= +Int_t AliDHN::MulToOne(Double_t *x) +{ + // Calculate the 1-dim index n from n-dim vector x, representing the + // abscissa of the n-dim histogram. The result will be between 0 and + // GetNbinsX()-1. + + Int_t k[fMaxNdim]; + for (int i=0; imkdir(GetName()); + if (dest) { + dest->cd(); + nbytes += histo.Write("histo"); + char axnam[1000]; + for (int i=0; iGetPath()); + dest->cd(".."); + } + return nbytes; +} +//============================================================================= +AliDHN *AliDHN::ProjectAlong(char *nam, Int_t dim, Int_t first, Int_t last) +{ + // Reduce dimension dim by summing up its bins between first and last. + // Use root convention: bin=1 is the first bin, bin=nbins is the last. + // Return the resulting fNdim-1 dimensional histogram. + + if (dim<0 || dim>fNdim-1) return 0; + if (first<0) first = 1; + if (last<0) last = fNbins[dim]; + + // create new (reduced) histogram + + TAxis *ax[fMaxNdim-1]; + int n=0; + for (int i=0; ilast) continue; + int n = 0; + for (int j=0; jMulToOne(m); + his->AddBinContent(n+1,GetBinContent(i+1)); + eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1)); + } + + // combine content and errors in one histogram + + for (int i=0; iGetNbinsX(); i++) { + his->SetBinError(i+1,TMath::Sqrt(eis->GetBinContent(i+1))); + } + + his->SetLineColor(this->GetLineColor()); + his->SetFillColor(this->GetFillColor()); + his->SetMarkerColor(this->GetMarkerColor()); + his->SetMarkerStyle(this->GetMarkerStyle()); + + // some cleanup + + delete eis; + return his; +} +//============================================================================= +TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Int_t *first, Int_t *last) +{ + // Project on dimension dim. Use only bins between first[i] and last[i]. + // Use root convention: bin=1 is the first bin, bin=nbins is the last. + // first[i]=0 means from the first bin + // last[i]=0 means till the last bin + + if (dim<0 || dim>fNdim-1) return 0; + + // calculate the projection; lowest index 0 + + double *yy = new double[fNbins[dim]]; // value + double *ey = new double[fNbins[dim]]; // error + for (int i=0; ilast[j]) {isgood = 0; break;} + } + if (isgood) { + yy[k[dim]]+=GetBinContent(i+1); + ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1); + } + } + + // make the projection histogram + // use name nam if specified; otherwise generate one + + TH1D *his; + char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim); + if (fAxis[dim].IsVariableBinSize()) + his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray()); + else + his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax()); + his->SetXTitle(fAxis[dim].GetTitle()); + // or his->GetXaxis()->ImportAttributes(ax); + his->Sumw2(); + his->SetLineColor(this->GetLineColor()); + his->SetFillColor(this->GetFillColor()); + his->SetMarkerColor(this->GetMarkerColor()); + his->SetMarkerStyle(this->GetMarkerStyle()); + for (int i=0; iGetNbinsX(); i++) { + his->SetBinContent(i+1,yy[i]); + his->SetBinError(i+1,TMath::Sqrt(ey[i])); + } + + // some cleanup + + delete [] yy; + delete [] ey; + delete [] k; + // if (name!=nam) delete [] name; + + return his; +} +//============================================================================= +TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Double_t *first, Double_t *last) +{ + // Project on dimension dim. Use only bins between first[i] and last[i]. + + if (dim<0 || dim>fNdim-1) return 0; + Int_t kfirst[fMaxNdim]; + Int_t klast[fMaxNdim]; + for (int i=0; ifNdim-1) return 0; + if (dim1<0 || dim1>fNdim-1) return 0; + + // calculate the projection + + double **yy = new double*[fNbins[dim0]]; // value + double **ey = new double*[fNbins[dim0]]; // error + for (int i=0; ilast[j]) {isgood = 0; break;} + } + if (isgood) { + yy[k[dim0]][k[dim1]]+=GetBinContent(i+1); + ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1); + } + } + + // make the projection histogram + // use name nam if specified; otherwise generate one + + TH2D *his=0; + char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0); + if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) + his = new TH2D(name,name, + fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), + fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); + else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize()) + his = new TH2D(name,name, + fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), + fNbins[dim1],fAxis[dim1].GetXbins()->GetArray()); + else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) + his = new TH2D(name,name, + fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(), + fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); + else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize()) + his = new TH2D(name,name, + fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(), + fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax()); + his->SetXTitle(fAxis[dim0].GetTitle()); + his->SetYTitle(fAxis[dim1].GetTitle()); + // or his->GetXaxis()->ImportAttributes(ax0); + // or his->GetYaxis()->ImportAttributes(ax1); + his->Sumw2(); + his->SetLineColor(this->GetLineColor()); + his->SetFillColor(this->GetFillColor()); + his->SetMarkerColor(this->GetMarkerColor()); + his->SetMarkerStyle(this->GetMarkerStyle()); + for (int i=0; iGetNbinsX(); i++) + for (int j=0; jGetNbinsY(); j++) { + his->SetBinContent(i+1,j+1,yy[i][j]); + his->SetBinError(i+1,j+1,TMath::Sqrt(ey[i][j])); + } + + // some cleanup + + for (int i=0; i 2007 + +#ifndef AliDHN_H +#define AliDHN_H + +#include +class TH2D; +class TAxis; + +const Int_t fMaxNdim=10; // maximum number of dimensions + +//============================================================================= +class AliDHN : public TH1D { + + public: + AliDHN() : TH1D(), fNdim(0) {printf("AliDHN object created\n");} + AliDHN(Char_t *nam, Int_t ndim, TAxis **ax); // constructor from scratch + AliDHN(Char_t *filename, Char_t *name); // constructor from file + virtual ~AliDHN() {printf("AliDHN object %s deleted\n",GetName());} + Int_t GetNdim() const {return fNdim;} + TAxis *GetAxis(Int_t i) {return &fAxis[i];} + void Fill(Double_t *xx, Double_t y=1); // fill histo + void Fill(Double_t x0=0, Double_t x1=0, + Double_t x2=0, Double_t x3=0, + Double_t x4=0, Double_t x5=0, + Double_t x6=0, Double_t x7=0, + Double_t x8=0, Double_t x9=0, + Double_t x10=0) { // fill histo; fNdim-th arg is weight + Double_t xx[fMaxNdim+1] = {x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10}; + Fill(xx,xx[fNdim]);} + Int_t Write(); // save histo and axis on file + Int_t Write(const char *, Int_t, Int_t) {return Write();} + // project along (integrate over) one axis + AliDHN *ProjectAlong(char *nam, Int_t dim, Int_t first=-1, Int_t last=-1); + // project on 1-dim histogram + TH1D *ProjectOn(char *nam, Int_t dim, Int_t *first=0, Int_t *last=0); + // project on 1-dim histogram + TH1D *ProjectOn(char *nam, Int_t dim, Double_t *first, Double_t *last); + // project on 2-dim histogram + TH2D *ProjectOn(char *nam, Int_t dim0, Int_t dim1, Int_t *first=0, Int_t *last=0); + + protected: + + Int_t fNdim; // number of dimensions + TAxis fAxis[fMaxNdim]; // axes + Int_t fNbins[fMaxNdim]; // {fAxis[0]->GetNbins(),fAxis[1]->... + Int_t fMbins[fMaxNdim]; // {...[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1} + + static Int_t Albins(Int_t n, TAxis **ax); // product of nbins of ax[0]...ax[n-1] + Int_t MulToOne(Int_t *k) const; // calc 1-dim index from n-dim indices + Int_t MulToOne(Double_t *x); // calc 1-dim index from n-dim vector + void OneToMul(Int_t n, Int_t *k); // calc n-dim indices from 1-dim index + + ClassDef(AliDHN,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/AliDLoop.cxx b/UNICOR/AliDLoop.cxx new file mode 100644 index 00000000000..f9c00600c5e --- /dev/null +++ b/UNICOR/AliDLoop.cxx @@ -0,0 +1,231 @@ +//Author: Dariusz Miskowiec 2007 + +//============================================================================= +// simple event loop manager +// First, events are classified according to centrality, reaction plane angle, +// and z-vertex. Then, single particles and true pairs are processed. Finally, +// event mixing is done within event classes. For this, the tree and the event +// passed as arguments are cloned. +// The class is, at present, deliberately left in a macro-like design, with all +// the essential things happening inside the function AliDLoop::Run. This is +// subject to future polishing if the class proves useful. +//============================================================================= + +#include +#include +#include +#include +#include "AliDLoop.h" +#include "AliDAnalGlobal.h" +#include "AliDAnalSingle.h" +#include "AliDAnalCorrel.h" +#include "AliDAnalPtfluc.h" +#include "AliDEvent.h" + +ClassImp(AliDLoop) + +//============================================================================= +AliDLoop::AliDLoop(TTree *tr, AliDEvent *ev0, char *outfil): + TObject(), + fTree0(tr), fTree1(0x0), + fEv0(ev0), fEv1(0x0), + fOutputFilename(outfil) +{ + // constructor + // Clone the tree such that during the event mixing the two events get + // populated from two different trees. This costs some memory (37 MB for + // ceres) but allows to avoid switching one tree between two different + // events, especially difficult for alice because alice trees do not like + // to be attached more than once. + + fTree1 = (TTree*) tr->Clone(); + fEv1 = (AliDEvent*) ev0->Clone(); + fEv0->AttachTree(fTree0); + fEv1->AttachTree(fTree1); +} +//============================================================================= +AliDLoop::AliDLoop(const AliDLoop &loop) : + TObject(loop), fTree0(loop.fTree0), fTree1(loop.fTree1), fEv0(loop.fEv0), + fEv1(loop.fEv1), fOutputFilename(loop.fOutputFilename) +{ + //copy constructor, shallow copy, just to make the compiler happy +} +//============================================================================= +AliDLoop &AliDLoop::operator=(const AliDLoop &source) +{ + // substitution operator, shallow copy, just to make the compiler happy + + TObject::operator=(source); + fTree0 = source.fTree0; + fTree1 = source.fTree1; + fEv0 = source.fEv0; + fEv1 = source.fEv1; + fOutputFilename = source.fOutputFilename; + return *this; +} +//============================================================================= +Int_t AliDLoop::Mem() const +{ + // return memory in MB occupied by this process + + ProcInfo_t info; + gSystem->GetProcInfo(&info); + Int_t memmb = info.fMemVirtual/1024; + return memmb; +} +//============================================================================= +void AliDLoop::Run(int nwish) +{ + // process nwish events from fTree0; nwish=-1 (default) means process all + + // control parameters + + + // const int kcen = 10; // number of centrality bins for mixing + // const int kphi = 8; // number of phi bins for mixing + // const int kzve = 13; // number of z-vertex bins for mixing + const int kcen = 1; // number of centrality bins for mixing + const int kphi = 1; // number of phi bins for mixing + const int kzve = 1; // number of z-vertex bins for mixing + int mixingFactor = 1; // number of wished event pairs/number of events + int mixingStep[10] = {2,3,5,7,11,13,17,19,23,29}; + if (mixingFactor>=10) {printf("mixing factor too high\n"); exit(-1);} + int minn = 20*mixingFactor; // minimum number of events per bin + int pr = !gROOT->IsBatch(); // more printing in interactive mode + + // initialize analysis + + Double_t etami = fEv0->Etamin(); + Double_t etama = fEv0->Etamax(); + AliDAnalGlobal *dag = new AliDAnalGlobal("dag"); + AliDAnalSingle *all = new AliDAnalSingle("all",etami,etama,0); + AliDAnalSingle *pim = new AliDAnalSingle("pim",etami,etama,-211); + AliDAnalSingle *pip = new AliDAnalSingle("pip",etami,etama, 211); + AliDAnalCorrel *cnn = new AliDAnalCorrel("cnn",etami,etama,-211,-211); + AliDAnalCorrel *cpp = new AliDAnalCorrel("cpp",etami,etama, 211, 211); + AliDAnalPtfluc *ptf = new AliDAnalPtfluc("ptf",-211,-211); + + // initialize number of events etc + + printf("\ndetermining number of events in chain..."); + int nmax = fTree0->GetEntries(); + int nact = nmax; + if (nwish>-1 && nwishGetEvent(i); + if (!fEv0->Good()) continue; + int icen = (int) (kcen*fEv0->Centrality()); + int iphi = (int) (kphi*(fEv0->RPphi()/TMath::Pi()+1)/2); + int izve = (int) (kzve*(fEv0->Zver()+1.0)/2.0); + if (pr) printf("\revent %5d %2d%2d%2d",i,icen,iphi,izve); + if (icen<0 || icen>=kcen) continue; + if (iphi<0 || iphi>=kphi) continue; + if (izve<0 || izve>=kzve) continue; + ar[icen][iphi][izve]->AddAt(i,nar[icen][iphi][izve]); + nar[icen][iphi][izve]++; + } + if (pr) printf("\n"); + + printf("%4d MB; shrinking event index arrays\n",Mem()); + int nall = 0; + int nsup = 0; + for (int i=0; iSet(nar[i][j][k]); + } + printf("%4d MB; %d events sorted into bins\n",Mem(),nall); + printf("%d events suppressed because bin occupancy was below %d\n",nsup,minn); + printf("sorting events took %.1f s CPU %.1f s real %.1f ms real per event\n\n", + sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall); + if ((nall-=nsup) < 1) return; + + // actual loop + + printf("%4d MB; starting the main loop\n",Mem()); + sto.Reset(); + sto.Start(); + for (int i=0; iGetSize(); // number of events in this bin + for (int l=0; lAt(l); + fTree0->GetEvent(m); + dag->Process(fEv0); + all->Process(fEv0); + pim->Process(fEv0); + pip->Process(fEv0); + cnn->Process(0,fEv0,fEv0,0); + cnn->Process(2,fEv0,fEv0,TMath::DegToRad()*180); + cpp->Process(0,fEv0,fEv0,0); + cpp->Process(2,fEv0,fEv0,TMath::DegToRad()*180); + ptf->Process(0,fEv0,fEv0); + } + if (pr && nevents) printf("\n"); + + // mixed pairs + // loop somewhat optimized to reduce jumping back and forth + + for (int imix=0; imixAt(l); + int next = (l+step)%nevents; + int m1 = ar[i][j][k]->At(next); + if (pr) printf("\rmixing %5d and %5d",m0,m1); + fTree0->GetEvent(m0); + fTree1->GetEvent(m1); + //printf(" event0 %f event1 %f\n",fEv0->ParticleP(0),fEv1->ParticleP(0)); + cnn->Process(1,fEv0,fEv1,0); + cpp->Process(1,fEv0,fEv1,0); + ptf->Process(1,fEv0,fEv1); + } + if (pr && nevents) printf("\r"); + } + } + + sto.Stop(); + printf("event loop took %.1f s CPU %.1f s real %.1f ms real per event\n\n", + sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall); + + // save analysis results + + dag->Save(fOutputFilename.Data(),"recreate"); + all->Save(fOutputFilename.Data()); + pim->Save(fOutputFilename.Data()); + pip->Save(fOutputFilename.Data()); + cnn->Save(fOutputFilename.Data()); + cpp->Save(fOutputFilename.Data()); + ptf->Save(fOutputFilename.Data()); +} +//============================================================================= diff --git a/UNICOR/AliDLoop.h b/UNICOR/AliDLoop.h new file mode 100644 index 00000000000..95993c0da63 --- /dev/null +++ b/UNICOR/AliDLoop.h @@ -0,0 +1,34 @@ +// Author: Dariusz Miskowiec 2007 + +#ifndef ALIDLOOP_H +#define ALIDLOOP_H + +#include + +class TTree; +class AliDEvent; + +//============================================================================= +class AliDLoop : public TObject { + + public: + // constructor + AliDLoop(TTree *tr, AliDEvent *ev0, char *outfil="result.root"); + AliDLoop(const AliDLoop &loop); // copy constructor + AliDLoop &operator=(const AliDLoop &loop); // substitution operator + virtual ~AliDLoop() {} // destructor + void Run(int n=-1); // process n events; n=-1 means all + + protected: + TTree *fTree0; //! input event tree + TTree *fTree1; //! clone of fTree0 for event mixing + AliDEvent *fEv0; //! event + AliDEvent *fEv1; //! clone of fEv0 for event mixing + TString fOutputFilename; //! output filename + Int_t Mem() const; // virtual memory in MB + + ClassDef(AliDLoop,1) +}; +//============================================================================= + +#endif diff --git a/UNICOR/AliDPair.cxx b/UNICOR/AliDPair.cxx new file mode 100644 index 00000000000..e160d225453 --- /dev/null +++ b/UNICOR/AliDPair.cxx @@ -0,0 +1,21 @@ +// Author: Dariusz Miskowiec 2002 + +//============================================================================= +// particle (track) pair +// Allows to calculate the kinematic pair variables typically used in +// two-particle correlation analyses. +//============================================================================= + +#include "AliDPair.h" + +ClassImp(AliDPair) + +//============================================================================= +AliDPair::AliDPair() : p0(), p1(), p(), q(), beta(), betat(), betaz(), ubeta(), + ubetat(), ubetaz(), CMp(), CMq(), buf() +{ + // constructor + + printf("AliDPair object created\n"); +} +//============================================================================= diff --git a/UNICOR/AliDPair.h b/UNICOR/AliDPair.h new file mode 100644 index 00000000000..c6e44d561a1 --- /dev/null +++ b/UNICOR/AliDPair.h @@ -0,0 +1,60 @@ +// Author: Dariusz Miskowiec 2002 + +#ifndef ALIDPAIR_H +#define ALIDPAIR_H + +#include +#include "TLorentzVector.h" + +//============================================================================= +class AliDPair : public TObject { + + protected: + TLorentzVector p0; // LAB four-momentum of track 0 + TLorentzVector p1; // LAB four-momentum of track 1 + TLorentzVector p; // LAB total four-momentum + TLorentzVector q; // LAB four-momentum difference + TVector3 beta; // LAB pair velocity + TVector3 betat; // LAB pair velocity transverse + TVector3 betaz; // LAB pair velocity along z + TVector3 ubeta; // LAB pair velocity direction + TVector3 ubetat; // LAB pair velocity transverse direction + TVector3 ubetaz; // LAB pair velocity along z direction (hm) + TLorentzVector CMp; // CM total four-momentum + TLorentzVector CMq; // CM four-momentum difference + TLorentzVector buf; // dummy buffer for swapping + + public: + AliDPair(); // constructor + virtual ~AliDPair() {printf("AliDPair object deleted\n");} + void Set0(double m,double p,double theta,double phi) {p0.SetXYZM(0,0,p,m); p0.SetTheta(theta); p0.SetPhi(phi);} + void Set1(double m,double p,double theta,double phi) {p1.SetXYZM(0,0,p,m); p1.SetTheta(theta); p1.SetPhi(phi);} + void SetMXYZ0(double m,double px,double py,double pz) {p0.SetXYZM(px,py,pz,m);} + void SetMXYZ1(double m,double px,double py,double pz) {p1.SetXYZM(px,py,pz,m);} + void CalcLAB() {p=p0+p1; q=p1-p0; beta=p.BoostVector(); + betaz.SetXYZ(0,0,beta.Z()); betat=beta; betat.SetZ(0); + ubeta=beta.Unit(); ubetat=betat.Unit(); ubetaz=betaz.Unit();} + double Rapidity() {return p.Rapidity();} + double Pt() {return p.Pt();} + double Phi() {return p.Phi();} + double DTheta() {return p1.Theta()-p0.Theta();} + double DPhi() {return TVector2::Phi_mpi_pi(p1.Phi()-p0.Phi());} + void CalcPairCM() {CMp=p; CMp.Boost(-beta); CMq=q; CMq.Boost(-beta);} + void CalcLcmsCM() {CMp=p; CMp.Boost(-betaz); CMq=q; CMq.Boost(-betaz);} + void Swap() {buf=p0; p0=p1; p1=buf; q=-q; CMq=-CMq;} + double Minv() {return p.M();} + double Qinv2() {return -q.M2();} + double QCM() {return CMq.Vect().Mag();} + double QCMpar() {return CMq.Vect()*ubeta;} + double QCMper() {return sqrt(fabs(QCM()*QCM()-QCMpar()*QCMpar()));} + double QCMout() {return +CMq.Vect().X()*ubetat.X()+CMq.Vect().Y()*ubetat.Y();} + double QCMside() {return -CMq.Vect().X()*ubetat.Y()+CMq.Vect().Y()*ubetat.X();} + double QCMlong() {return CMq.Vect().Z();} + double QCMTheta() {return CMq.Theta();} + double QCMPhi() {return CMq.Phi();} + double QCMPhiOut() {return TMath::ATan2(QCMside(),QCMout());} // phi w.r.t. out + + ClassDef(AliDPair,1) +}; +//============================================================================= +#endif diff --git a/UNICOR/README b/UNICOR/README new file mode 100644 index 00000000000..bb74b957c25 --- /dev/null +++ b/UNICOR/README @@ -0,0 +1,29 @@ +The UNICOR package performs a universal two-particle analysis. The +two-particle analysis classes like AliDAnalCorrel or AliDAnalPtfluc +are complemented by standard event variable and single particle +analyses AliDAnalGlobal and AliDAnalSingle. The package is universal +in the sense that the same analysis runs on data from ALICE, CERES, +CBM... The files in this directory have been automatically converted +to the aliroot standard (Ali prefix was added to all class and file +names) and thus should not be edited. + +Below are two examples of how to run the analysis on ESD events +(after editing the input data path in run-alone.C and run-on-train.C). + +1) root run-alone.C + +A standalone (non-train) analysis of ESD events, including event +mixing. AliDLoop class is doing the event loop and mixing. + +2) root run-on-train.C + +No event mixing is done. + +In both cases an (identical) output file unicor-result.root should +be produced, containing the (multidimensional) histograms coming out +of the analysis. + +In the second case in addition a file kuku.root is produced, containing +the analysis objects. This is a meaningless by-product. + +Dariusz Miskowiec 2008 diff --git a/UNICOR/libUNICOR.pkg b/UNICOR/libUNICOR.pkg new file mode 100644 index 00000000000..27652d45c64 --- /dev/null +++ b/UNICOR/libUNICOR.pkg @@ -0,0 +1,8 @@ +#-*- Mode: Makefile -*- + +SRCS = $(shell cd UNICOR; ls AliD*.cxx) + +HDRS = $(SRCS:.cxx=.h) + +CINTAUTOLINK = 1 + diff --git a/UNICOR/makechain.C b/UNICOR/makechain.C new file mode 100644 index 00000000000..5917ac75e02 --- /dev/null +++ b/UNICOR/makechain.C @@ -0,0 +1,38 @@ +// Author: Dariusz Miskowiec 2008 + +//============================================================================= +TChain *makechainexp(char *exp, char *path, int nfil=1000000, int nskip=0) { + char treename[1000]; + if (strcmp(exp,"ceres3c2")==0) sprintf(treename,"T"); + if (strcmp(exp,"aliceesd")==0) sprintf(treename,"esdTree"); + if (strcmp(exp,"cbm")==0) sprintf(treename,"cbmsim"); + return makechain(treename, path, nfil, nskip); +} +//============================================================================= +TChain *makechain(char *name, char *path, int nfil=1000000, int nskip=0) { + + // make chain of root files + // if path ends with "root" then add all these files; + // otherwise interprete path as the list of files and add nfil files + // after skipping nskip + + printf("path=%s\n",path); + + TChain *chain = new TChain(name); + TString str(path); + if (str.EndsWith("root")) chain->Add(path); + else { + fstream ascii_in; + ascii_in.open(path, ios::in); + char filnam[1000]; + for (int i=0; i> filnam; + if (ascii_in.eof()) break; + if (i>=nskip) chain->Add(filnam); + } + } + chain->Lookup(); + chain->GetListOfFiles()->Print(); + return chain; +} +//============================================================================= diff --git a/UNICOR/run-alone.C b/UNICOR/run-alone.C new file mode 100644 index 00000000000..fb61b4fd8cd --- /dev/null +++ b/UNICOR/run-alone.C @@ -0,0 +1,16 @@ +{ +gROOT->LoadMacro("makechain.C"); +tr = makechain("esdTree","/u/sma/data/mc/v4-13-Rev-01/fulltrd_pdc_0.txt",10); + +gSystem->Load("libEG.so"); +gSystem->Load("libTree.so"); +gSystem->Load("libVMC.so"); +gSystem->Load("libSTEERBase.so"); +gSystem->Load("libESD.so"); +gSystem->Load("libANALYSIS"); +gSystem->Load("libUNICOR.so"); + +d0=new AliDEventAliceESD(); +lo = new AliDLoop(tr,d0,"unicor-result.root"); +lo->Run(); +} diff --git a/UNICOR/run-on-train.C b/UNICOR/run-on-train.C new file mode 100644 index 00000000000..6f45fd8dcc2 --- /dev/null +++ b/UNICOR/run-on-train.C @@ -0,0 +1,26 @@ +{ +gROOT->LoadMacro("makechain.C"); +tr = makechain("esdTree","/u/sma/data/mc/v4-13-Rev-01/fulltrd_pdc_0.txt",10); + +gSystem->Load("libEG.so"); +gSystem->Load("libPhysics.so"); +gSystem->Load("libTree.so"); +gSystem->Load("libVMC.so"); +gSystem->Load("libSTEERBase.so"); +gSystem->Load("libESD.so"); +gSystem->Load("libANALYSIS"); +gSystem->Load("libUNICOR.so"); + +AliAnalysisManager *mgr = new AliAnalysisManager("TestManager"); +AliVEventHandler* esdH = new AliESDInputHandler; +mgr->SetInputEventHandler(esdH); +AliDAnalysisTask *mytask = new AliDAnalysisTask(); +mgr->AddTask(mytask); +AliAnalysisDataContainer *cinput = mgr->CreateContainer("cinput",TChain::Class(),AliAnalysisManager::kInputContainer); +AliAnalysisDataContainer *coutpt = mgr->CreateContainer("clist1", TList::Class(),AliAnalysisManager::kOutputContainer,"kuku.root"); +mgr->ConnectInput (mytask,0,cinput); +mgr->ConnectOutput(mytask,0,coutpt); +mgr->InitAnalysis(); +mgr->PrintStatus(); +mgr->StartAnalysis("local",tr); +}