// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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 <TROOT.h>
#include <TFile.h>
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
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
gROOT->cd();
f->Close();
}
-/******************************************************************************/
+//=============================================================================
#ifndef ALIDANAL_H
#define ALIDANAL_H
-///////////////////////////////////////////////////////////////////////////////
-// Parent analysis class
-///////////////////////////////////////////////////////////////////////////////
-
#include <TNamed.h>
#include <TObjArray.h>
#include <TDatabasePDG.h>
-/*****************************************************************************/
+//=============================================================================
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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005
+
+//=============================================================================
+// two-particle correlation analyzer
+//=============================================================================
+
+#include <TROOT.h>
+#include <TRandom2.h>
+#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 (j<i), mixed - all
+ // thus, proper rotation is either by 180, or by 170 AND 190, etc.
+
+ double cent = (ev0->Centrality()+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; i<ev0->NParticles(); i++) {
+ if (!ev0->ParticleGood(i,fPid0)) continue;
+ for (int j=0; j<ev1->NParticles(); j++) {
+ if (ev0==ev1 && j<i && fPid0==fPid1 ) continue;
+ if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
+ if (!ev1->ParticleGood(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
+ }
+ }
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
+
+//=============================================================================
+// event global variable analyzer
+//=============================================================================
+
+#include <TROOT.h>
+#include <TH1.h>
+#include <TH2.h>
+#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);
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2008
+
+//=============================================================================
+// pt-fluctuations analyzer
+//=============================================================================
+
+#include <TROOT.h>
+#include <TRandom2.h>
+#include <TVector2.h>
+#include <TMath.h>
+#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; i<ev0->NParticles(); 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; j<ev1->NParticles(); j++) {
+ if (ev0==ev1 && j==i) continue;
+ if (ev0==ev1 && j<i && fPid0==fPid1) continue;
+ if (!ev1->ParticleGood(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);
+ }
+ }
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
+
+//=============================================================================
+// single particle analyzer
+//=============================================================================
+
+#include <TROOT.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TParticlePDG.h>
+#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; i<ev->NParticles(); 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);
+ }
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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
--- /dev/null
+//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<TTree*> (GetInputData(0));
+ if (!tree) {
+ Printf("ERROR: Could not read chain from input slot 0");
+ } else {
+ fEv0->AttachTree(tree);
+ }
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
+ (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");
+}
+/*****************************************************************************/
--- /dev/null
+// 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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
+
+//=============================================================================
+// parent class of all events; analyzers access data via this class
+//=============================================================================
+
+#include <TObject.h>
+#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; i<NParticles(); i++) {
+ if (!ParticleGood(i,0)) continue;
+ double pt = ParticlePt(i);
+ if (pt>2.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;
+}
+//=============================================================================
+
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
+
+#ifndef ALIDEVENT_H
+#define ALIDEVENT_H
+
+#include <TObject.h>
+#include <TMath.h>
+
+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
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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;
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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 <TFile.h>
+#include <TDirectory.h>
+#include <TAxis.h>
+#include <TH2.h>
+#include <TMath.h>
+#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; i<fNdim; i++) ax[i]->Copy(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; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
+ for (int i=fNdim-1; i>0; 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; fNdim<fMaxNdim; fNdim++) {
+ sprintf(axnam,"axis%d",fNdim);
+ ax[fNdim] = (TAxis *) gDirectory->Get(axnam);
+ if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]);
+ else break;
+ }
+ f->Close();
+
+ fMbins[fNdim-1] = 1;
+ for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
+ for (int i=fNdim-1; i>0; 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; i<n; i++) nbins *= ax[i]->GetNbins();
+ 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<fNdim; i++) {
+ if (k[i]<0) return -1;
+ if (k[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; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1;
+ return MulToOne(k);
+}
+//=============================================================================
+void AliDHN::OneToMul(Int_t n, Int_t *k)
+{
+ // Calculate the n-dim indices k[fNdim] from 1-dim index n.
+ // Valid n should be between 0 and GetNbinsX()-1.
+ // Valid results will be between 0 and fNbins[i]-1.
+
+ div_t idi; // integer division structure
+ for (int i=0; i<fNdim; i++) {
+ idi = div(n,fMbins[i]);
+ k[i] = idi.quot; // quotient
+ n = idi.rem; // remainder
+ }
+}
+//=============================================================================
+void AliDHN::Fill(Double_t *xx, Double_t w)
+{
+ // Fill the histogram. The array xx holds the abscissa information, w is the
+ // weigth. The 1-dim histogram is filled using the standard Fill method
+ // (direct access to the arrays was tried and was not faster).
+
+ int nbin = MulToOne(xx);
+ if (nbin != -1) TH1D::Fill(nbin+1,w);
+}
+//=============================================================================
+Int_t AliDHN::Write()
+{
+ // Save the 1-dim histo and the axes in a subdirectory on file. This might
+ // not be the most elegant way but it is very simple and backward and forward
+ // compatible.
+
+ Int_t nbytes = 0;
+ TH1D histo(*this);
+ histo.SetName("histo"); // hadd bug fix
+ TDirectory *dest = gDirectory->mkdir(GetName());
+ if (dest) {
+ dest->cd();
+ nbytes += histo.Write("histo");
+ char axnam[1000];
+ for (int i=0; i<fNdim; i++) {
+ sprintf(axnam,"axis%d",i);
+ nbytes += fAxis[i].Write(axnam);
+ }
+ printf(" %s stored in %s\n",GetName(),dest->GetPath());
+ 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; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i);
+ char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
+ char *eame = Form("%s_error",name);
+ AliDHN *his = new AliDHN(name,fNdim-1,ax); // result histogram
+ AliDHN *eis = new AliDHN(eame,fNdim-1,ax); // temporary storage for errors
+
+ // sum up the content and errors squared
+
+ Int_t *k = new Int_t[fNdim]; // old hist multiindex
+ Int_t *m = new Int_t[fNdim-1]; // new hist multiindex
+ for (int i=0; i<GetNbinsX(); i++) {
+ OneToMul(i,k);
+ if (k[dim]+1<first) continue;
+ if (k[dim]+1>last) continue;
+ int n = 0;
+ for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j];
+ n = his->MulToOne(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; i<his->GetNbinsX(); 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; i<fNbins[dim]; i++) yy[i]=0;
+ for (int i=0; i<fNbins[dim]; i++) ey[i]=0;
+ Int_t *k = new Int_t[fNdim];
+ for (int i=0; i<GetNbinsX(); i++) {
+ OneToMul(i,k);
+ int isgood = 1; // bin within the range?
+ for (int j=0; j<fNdim; j++) {
+ if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
+ if (last) if (last[j]) if (k[j]+1>last[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; i<his->GetNbinsX(); 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; i<fNdim; i++) {
+ kfirst[i] = fAxis[i].FindBin(first[i]);
+ klast[i] = fAxis[i].FindBin(last[i]);
+ }
+ return ProjectOn(nam,dim,kfirst,klast);
+}
+//=============================================================================
+TH2D *AliDHN::ProjectOn(char *nam, Int_t dim0, Int_t dim1, Int_t *first, Int_t *last)
+{
+ // Project on dim1 vs dim0. 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 (dim0<0 || dim0>fNdim-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; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]];
+ for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]];
+ for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0;
+ for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0;
+ Int_t *k = new Int_t[fNdim];
+ for (int i=0; i<GetNbinsX(); i++) {
+ OneToMul(i,k);
+ int isgood = 1; // bin within the range?
+ for (int j=0; j<fNdim; j++) {
+ if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
+ if (last) if (last[j]) if (k[j]+1>last[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; i<his->GetNbinsX(); i++)
+ for (int j=0; j<his->GetNbinsY(); 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<fNbins[dim0]; i++) delete [] yy[i];
+ for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i];
+ delete [] yy;
+ delete [] ey;
+ delete [] k;
+ // if (name!=nam) delete [] name;
+
+ return his;
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
+
+#ifndef AliDHN_H
+#define AliDHN_H
+
+#include <TH1.h>
+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
--- /dev/null
+//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 <TROOT.h>
+#include <TTree.h>
+#include <TSystem.h>
+#include <TStopwatch.h>
+#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 && nwish<nmax) nact = nwish;
+ printf("analyzing %d events of %d\n\n",nact,nmax);
+ TStopwatch sto;
+
+ sto.Reset();
+ sto.Start();
+
+ // find groups of of similar events
+ // first allocate nact entries to each array; later reduce
+
+ printf("%4d MB; booking event index arrays\n",Mem());
+ TArrayI *ar[kcen][kphi][kzve];
+ int nar[kcen][kphi][kzve];
+ for (int i=0; i<kcen; i++)
+ for (int j=0; j<kphi; j++)
+ for (int k=0; k<kzve; k++) {
+ ar[i][j][k] = new TArrayI(nact);
+ nar[i][j][k] = 0;
+ }
+ printf("%4d MB; filling event index arrays\n",Mem());
+ for (int i=0; i<nact; i++) {
+ fTree0->GetEvent(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; i<kcen; i++)
+ for (int j=0; j<kphi; j++)
+ for (int k=0; k<kzve; k++) {
+ nall += nar[i][j][k];
+ // suppress scarse bins
+ if (nar[i][j][k]<minn) {
+ nsup += nar[i][j][k];
+ nar[i][j][k]=0;
+ }
+ ar[i][j][k]->Set(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; i<kcen; i++)
+ for (int j=0; j<kphi; j++)
+ for (int k=0; k<kzve; k++) {
+
+ // true pairs
+
+ int nevents = ar[i][j][k]->GetSize(); // number of events in this bin
+ for (int l=0; l<nevents; l++) {
+ if (pr) printf("\revent %6d of %6d %3d%3d%3d",l,nevents,i,j,k);
+ int m = ar[i][j][k]->At(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; imix<mixingFactor; imix++) {
+ int step = mixingStep[imix];
+ for (int phase=0; phase<step; phase++) for (int l=phase; l<nevents; l+=step) {
+ int m0 = ar[i][j][k]->At(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());
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec 2007
+
+#ifndef ALIDLOOP_H
+#define ALIDLOOP_H
+
+#include <TObject.h>
+
+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
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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");
+}
+//=============================================================================
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2002
+
+#ifndef ALIDPAIR_H
+#define ALIDPAIR_H
+
+#include <TObject.h>
+#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
--- /dev/null
+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
--- /dev/null
+#-*- Mode: Makefile -*-
+
+SRCS = $(shell cd UNICOR; ls AliD*.cxx)
+
+HDRS = $(SRCS:.cxx=.h)
+
+CINTAUTOLINK = 1
+
--- /dev/null
+// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 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<nfil+nskip; i++) {
+ ascii_in >> filnam;
+ if (ascii_in.eof()) break;
+ if (i>=nskip) chain->Add(filnam);
+ }
+ }
+ chain->Lookup();
+ chain->GetListOfFiles()->Print();
+ return chain;
+}
+//=============================================================================
--- /dev/null
+{
+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();
+}
--- /dev/null
+{
+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);
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff