Q-bins variable size, finer kt binning.
Analysis frame can now be specified via argument list of AnalCorrel.
Two macros for running without and with event mixing provided; both
can run locally or on grid.
Cosmetics.
#include "AliUnicorAnalPtfluc.h"
#include "AliUnicorAnalHighpt.h"
#include "AliUnicorEventAliceESD.h"
+#include "AliESDEvent.h"
#include "AliAnalysisTaskUnicor.h"
ClassImp(AliAnalysisTaskUnicor)
fOutputList->Add(new AliUnicorAnalSingle("all",fEv0->Etamin(),fEv0->Etamax(),0));
fOutputList->Add(new AliUnicorAnalSingle("pim",fEv0->Etamin(),fEv0->Etamax(),-211));
fOutputList->Add(new AliUnicorAnalSingle("pip",fEv0->Etamin(),fEv0->Etamax(), 211));
- fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211));
- fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211));
- fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211));
+ int frame = AliUnicorAnalCorrel::kLCMS;
+ fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211, frame));
+ fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211, frame));
+ fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211, frame));
fOutputList->Add(new AliUnicorAnalPtfluc("ptf",0,0));
fOutputList->Add(new AliUnicorAnalHighpt("hpt",fEv0->Etamin(),fEv0->Etamax(),0,0));
}
// unicor analysis task
//=============================================================================
#include "AliMultiEventInputHandler.h"
-#include "AliESDHeader.h"
+#include "AliESDEvent.h"
#include "AliUnicorAnalGlobal.h"
#include "AliUnicorAnalSingle.h"
#include "AliUnicorAnalCorrel.h"
fOutputList->Add(new AliUnicorAnalSingle("all",fEv0->Etamin(),fEv0->Etamax(),0));
fOutputList->Add(new AliUnicorAnalSingle("pim",fEv0->Etamin(),fEv0->Etamax(),-211));
fOutputList->Add(new AliUnicorAnalSingle("pip",fEv0->Etamin(),fEv0->Etamax(), 211));
- fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211));
- fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211));
- fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211));
+ int frame = AliUnicorAnalCorrel::kLCMS;
+ fOutputList->Add(new AliUnicorAnalCorrel("cnn",fEv0->Etamin(),fEv0->Etamax(),-211,-211, frame));
+ fOutputList->Add(new AliUnicorAnalCorrel("cpp",fEv0->Etamin(),fEv0->Etamax(), 211, 211, frame));
+ fOutputList->Add(new AliUnicorAnalCorrel("cnp",fEv0->Etamin(),fEv0->Etamax(),-211, 211, frame));
fOutputList->Add(new AliUnicorAnalPtfluc("ptf",0,0));
fOutputList->Add(new AliUnicorAnalHighpt("hpt",fEv0->Etamin(),fEv0->Etamax(),0,0));
}
//=============================================================================
void AliAnalysisTaskUnicorME::UserExec(Option_t */*option*/)
{
- // process one event
+ // process one event pair
+
+ static int nic0=0;
+ static int nic1=0;
+ static int nic2=0;
+ static int nic3=0;
+
+ nic0++;
if (fInputHandler->GetBufferSize() < 2) return;
+ nic1++;
AliESDEvent *esd0 = dynamic_cast<AliESDEvent*>(GetEvent(0));
AliESDEvent *esd1 = dynamic_cast<AliESDEvent*>(GetEvent(1));
- if (!esd0) return;
- if (!esd1) return;
- if (!esd0->GetNumberOfTracks()) return;
- if (!esd1->GetNumberOfTracks()) return;
-
- printf("esd0 nr %3d mult %3d esd1 nr %3d mult %3d\n",
- esd0->GetEventNumberInFile(), esd0->GetNumberOfTracks(),
- esd1->GetEventNumberInFile(), esd1->GetNumberOfTracks());
+ if (!esd0) return;
+ nic2++;
+ // if (esd0->GetEventType() != 7) return; // physics event
fEv0->SetESD(esd0);
- fEv1->SetESD(esd1);
-
if (!fEv0->Good()) return;
+ nic3++;
((AliUnicorAnalGlobal *) fOutputList->At(0))->Process(fEv0);
((AliUnicorAnalSingle *) fOutputList->At(1))->Process(fEv0);
((AliUnicorAnalSingle *) fOutputList->At(2))->Process(fEv0);
((AliUnicorAnalPtfluc *) fOutputList->At(7))->Process(0,fEv0,fEv0);
((AliUnicorAnalHighpt *) fOutputList->At(8))->Process(fEv0,fEv0);
+ if (!esd1) return;
+ // if (esd1->GetEventType() != 7) return; // physics event
+ printf("esd0 nr %3d mult %3d esd1 nr %3d mult %3d\n",
+ esd0->GetEventNumberInFile(), esd0->GetNumberOfTracks(),
+ esd1->GetEventNumberInFile(), esd1->GetNumberOfTracks());
+ fEv1->SetESD(esd1);
if (!fEv1->Good()) return;
((AliUnicorAnalCorrel *) fOutputList->At(4))->Process(1,fEv0,fEv1,0);
((AliUnicorAnalCorrel *) fOutputList->At(5))->Process(1,fEv0,fEv1,0);
((AliUnicorAnalPtfluc *) fOutputList->At(7))->Process(1,fEv0,fEv1);
((AliUnicorAnalHighpt *) fOutputList->At(8))->Process(fEv0,fEv1);
+ printf("counts: %6d %6d %6d %6d\n",nic0,nic1,nic2,nic3);
PostData(1, fOutputList);
}
//=============================================================================
TDatabasePDG AliUnicorAnal::fgPDG;
//=============================================================================
-AliUnicorAnal::AliUnicorAnal(char *nam) : TNamed(nam,nam), fHistos() {
+AliUnicorAnal::AliUnicorAnal(const char *nam) : TNamed(nam,nam), fHistos() {
// constructor
class AliUnicorAnal : public TNamed {
public:
- AliUnicorAnal(char *nam="anal"); // constructor
+ AliUnicorAnal(const char *nam="anal"); // constructor
virtual ~AliUnicorAnal() {} // destructor
Long64_t Merge(const TCollection * const list); // sumup histograms
void Save(const char *outfil, const char *mode="update"); // save histograms
ClassImp(AliUnicorAnalCorrel)
//=============================================================================
-AliUnicorAnalCorrel::AliUnicorAnalCorrel(Char_t *nam, Double_t emi, Double_t ema,
- Int_t pid0, Int_t pid1):
- AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fPa()
+AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t ema,
+ Int_t pid0, Int_t pid1, Int_t frame):
+ AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fFrame(frame), fPa()
{
// constructor
// emi and ema define the rapidity range for histogram
//int npt = 7; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
//int npt = 6; double ptbins[]={0,0.1,0.25,0.35,0.55,1.0,2.0}; // like Adam, except last bin split
- int npt = 7; double ptbins[]={0,0.1,0.25,0.40,0.55,0.7,1.0,2.0}; // like Adam in Mar-2010, + first+last
+ //int npt = 7; double ptbins[]={0,0.1,0.25,0.40,0.55,0.7,1.0,2.0}; // like Adam in Mar-2010, + first+last
+ int npt = 8; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0,2.0}; // for Pb
- double qbins[100];
- for (int i=0;i<20;i++) qbins[i]=i*0.005;
- for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
+ double qbins[200] = {0};
+ // for (int i=0;i<20;i++) qbins[i]=i*0.005;
+ // for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
+ for (int i=0;i<30;i++) qbins[i]=i*0.010;
+ for (int i=0;i<35;i++) qbins[30+i]=0.3+i*0.02;
+ for (int i=0;i<20;i++) qbins[65+i]=1.0+i*0.05;
+ for (int i=0;i<11;i++) qbins[85+i]=2.0+i*0.20;
TAxis *ax[8];
ax[0] = new TAxis(3,-0.5,2.5);ax[0]->SetTitle("trumixrot");
ax[4] = new TAxis(npt,ptbins);ax[4]->SetTitle("kt (GeV/c)"); // pair pt/2
ax[5] = new TAxis(8,0,pi); ax[5]->SetTitle("q-theta");
ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
- //ax[7] = new TAxis(64,qbins); ax[7]->SetTitle("q (GeV/c)");
- ax[7] = new TAxis(150,0,3.0); ax[7]->SetTitle("q (GeV/c)");
+ ax[7] = new TAxis(95,qbins); ax[7]->SetTitle("q (GeV/c)");
+ // ax[7] = new TAxis(700,0,3.5); ax[7]->SetTitle("q (GeV/c)");
AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
for (int i=0; i<8; i++) delete ax[i];
fHistos.Add(pair);
// 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;
+ double q0x=0,q0y=0,q1x=0,q1y=0;
ev0->RP(q0x,q0y);
ev1->RP(q1x,q1y);
double rpphi = atan2(q0y+q1y,q0x+q1x);
ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot)) continue;
twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
fPa.CalcLAB();
- fPa.CalcPairCM();
- //fPa.CalcLcmsCM();
- if (fPa.QCM()==0) return; // should not be too frequent
+ if (fFrame == kPairFrame) fPa.CalcPairCM();
+ if (fFrame == kLCMS) fPa.CalcLcmsCM();
+ if (fPa.QCM()==0) {printf("AliUnicorAnalCorrel: Q=0\n"); return;} // should not be too frequent
double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
double weigth = 1.0;
- // if (tmr==0) weigth = 1.0+0.5*exp(-fPa.QCM()*fPa.QCM()*1*1/0.197/0.197);
+ // if (tmr==0 && fPid0==fPid1) weigth = 1.0+0.6*exp(-pow(fPa.QCM()*5/0.197,2));
+ // if (tmr==0 && fPid0==fPid1) weigth = 1.0+0.6*exp(-fPa.QCM()*2/0.197);
pair->Fill((double) tmr, // 0 for tru, 1 for mix, 2 for rot
cent, // centrality
fPa.Rapidity(), // pair rapidity
class AliUnicorAnalCorrel : public AliUnicorAnal {
public:
- AliUnicorAnalCorrel(Char_t *nam="correl",
- Double_t emi=-1, Double_t ema=1,
- Int_t pid0=0, Int_t pid1=0); // constructor
- virtual ~AliUnicorAnalCorrel(){} // destructor
+ enum AnalysisFrames {kPairFrame, kLCMS, kLAB};
+ AliUnicorAnalCorrel(const char *nam="correl", Double_t emi=-1, Double_t ema=1,
+ Int_t pid0=0, Int_t pid1=0, Int_t frame=0); // constructor
+ virtual ~AliUnicorAnalCorrel(){} // destructor
// process one (tru) or two (mix) events
void Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot);
Int_t fPid1; // particle species 1
Double_t fMass0; // mass 0
Double_t fMass1; // mass 1
+ Int_t fFrame; // analysis frame
AliUnicorPair fPa; // pair buffer for calculations
ClassDef(AliUnicorAnalCorrel,1)
ClassImp(AliUnicorAnalGlobal)
//=============================================================================
-AliUnicorAnalGlobal::AliUnicorAnalGlobal(Char_t *nam) : AliUnicorAnal(nam)
+AliUnicorAnalGlobal::AliUnicorAnalGlobal(const char *nam) : AliUnicorAnal(nam)
{
// constructor
mult->Fill(ev->NGoodParticles(),1.0);
cent->Fill(ev->Centrality(),1.0);
cemu->Fill(ev->Centrality(),ev->NGoodParticles());
- Double_t qx,qy;
+ Double_t qx=0,qy=0;
ev->RP(qx,qy);
dire->Fill(qx,qy,1.0);
zver->Fill(ev->Zver(),1.0);
class AliUnicorAnalGlobal : public AliUnicorAnal {
public:
- AliUnicorAnalGlobal(Char_t *nam="global"); // constructor
- virtual ~AliUnicorAnalGlobal(){} // destructor
- void Process(AliUnicorEvent *ev) const; // fill histograms
+ AliUnicorAnalGlobal(const char *nam="global"); // constructor
+ virtual ~AliUnicorAnalGlobal(){} // destructor
+ void Process(AliUnicorEvent *ev) const; // fill histograms
ClassDef(AliUnicorAnalGlobal,1)
};
ClassImp(AliUnicorAnalHighpt)
//=============================================================================
-AliUnicorAnalHighpt::AliUnicorAnalHighpt(Char_t *nam, Double_t emi, Double_t ema, Int_t pid0,
+AliUnicorAnalHighpt::AliUnicorAnalHighpt(const char *nam, Double_t emi, Double_t ema, Int_t pid0,
Int_t pid1): AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1)
{
// constructor
// mixing-proof centrality and reaction plane angle
double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
- double q0x,q0y,q1x,q1y;
+ double q0x=0,q0y=0,q1x=0,q1y=0;
ev0->RP(q0x,q0y);
ev1->RP(q1x,q1y);
double rpphi = atan2(q0y+q1y,q0x+q1x);
class AliUnicorAnalHighpt : public AliUnicorAnal {
public:
- AliUnicorAnalHighpt(Char_t *nam="highpt",
+ AliUnicorAnalHighpt(const char *nam="highpt",
Double_t emi=-1, Double_t ema=1,
Int_t pid0=0, Int_t pid1=0); // constructor
virtual ~AliUnicorAnalHighpt(){} // destructor
ClassImp(AliUnicorAnalPtfluc)
//=============================================================================
-AliUnicorAnalPtfluc::AliUnicorAnalPtfluc(Char_t *nam, Int_t pid0, Int_t pid1) :
+AliUnicorAnalPtfluc::AliUnicorAnalPtfluc(const char *nam, Int_t pid0, Int_t pid1) :
AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1)
{
// constructor
class AliUnicorAnalPtfluc : public AliUnicorAnal {
public:
- AliUnicorAnalPtfluc(Char_t *nam="correl", Int_t pid0=0, Int_t pid1=0); // constructor
+ AliUnicorAnalPtfluc(const char *nam="correl", Int_t pid0=0, Int_t pid1=0); // constructor
virtual ~AliUnicorAnalPtfluc(){} // destructor
void Process(Int_t tmr, AliUnicorEvent * const ev0, AliUnicorEvent * const ev1); // process event(s)
ClassImp(AliUnicorAnalSingle)
//=============================================================================
-AliUnicorAnalSingle::AliUnicorAnalSingle(Char_t *nam, Double_t emi, Double_t ema, Int_t pid) :
+AliUnicorAnalSingle::AliUnicorAnalSingle(const char *nam, Double_t emi, Double_t ema, Int_t pid) :
AliUnicorAnal(nam), fPid(pid), fMass(0.0)
{
// constructor
class AliUnicorAnalSingle : public AliUnicorAnal {
public:
- AliUnicorAnalSingle(Char_t *nam="single",
+ AliUnicorAnalSingle(const char *nam="single",
Double_t emi=-1, Double_t ema=1,
Int_t pid=0); // constructor
virtual ~AliUnicorAnalSingle(){} // destructor
virtual Double_t RPphi() const = 0;
virtual Double_t Zver() const = 0; // z-vertex (-1,1)
virtual Int_t NParticles() const = 0; // number of tracks
- virtual Int_t NGoodParticles() const = 0; // number of good particles
- virtual Bool_t ParticleGood(Int_t i, Int_t pidi) const = 0;
+ virtual Bool_t ParticleGood(Int_t i, Int_t pidi=0) 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;
// toolkit part
+ Int_t NGoodParticles() const {int n=0; for (int i=0; i<NParticles(); i++) if (ParticleGood(i)) n++; return n;}
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));}
ClassImp(AliUnicorEventAliceESD)
//=============================================================================
- AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fESD(esd), fPhysicsSelection(0)
+ AliUnicorEventAliceESD::AliUnicorEventAliceESD(AliESDEvent *esd) : AliUnicorEvent(), fESD(esd)//, fPhysicsSelection(0)
{
// constructor
// printf("%s object created\n",ClassName());
if (!fESD) fESD = new AliESDEvent();
- fPhysicsSelection = new AliPhysicsSelection();
+ // fPhysicsSelection = new AliPhysicsSelection();
}
//=============================================================================
AliUnicorEventAliceESD::~AliUnicorEventAliceESD()
const AliESDVertex *vtx = fESD->GetPrimaryVertex();
if (!vtx->GetStatus()) return kFALSE;
if (fabs(Zver())>1) return kFALSE;
+
+ /*
+ int hard=0;
+ for (int i=0; i<NParticles(); i++) {
+ if (!ParticleGood(i)) continue;
+ if (ParticlePt(i)<1.0) continue;
+ hard = 1;
+ break;
+ }
+ return hard;
+ */
+
return kTRUE;
}
//=============================================================================
AliESDtrack *track = fESD->GetTrack(i);
if (!track->IsOn(AliESDtrack::kTPCrefit)) return 0; // TPC refit
if (!track->IsOn(AliESDtrack::kITSrefit)) return 0; // ITS refit
- if (track->GetTPCNcls() < 90) return 0; // number of TPC clusters
+ if (track->GetTPCNcls() < 70) return 0; // number of TPC clusters
if (track->GetKinkIndex(0) > 0) return 0; // no kink daughters
const AliExternalTrackParam *tp = GetTrackParam(i);
if (!tp) return 0; // track param
// double phi = ParticlePhi(i);
// if (eta>0 && phi>-8*pi9 && phi<-7*pi9) return 0; // A10
// if (eta>0 && phi> 1*pi9 && phi< 2*pi9) return 0; // A01
+ // if (tp->Pt()<0.2) return 0; // lower pt cutoff
Float_t r,z;
track->GetImpactParametersTPC(r,z);
if (pidi == -211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==-1;
else if (pidi == 211) return p[AliPID::kPion]+p[AliPID::kMuon]>0.5 && q==1;
-
else if (pidi == -321) return p[AliPID::kKaon]>0.5 && q==-1;
else if (pidi == 321) return p[AliPID::kKaon]>0.5 && q==1;
else if (pidi == -2212) return p[AliPID::kProton]>0.5 && q==-1;
else if (pidi == 2212) return p[AliPID::kProton]>0.5 && q==1;
+ else if (pidi == -11) return p[AliPID::kElectron]>0.5 && q==1;
+ else if (pidi == 11) return p[AliPID::kElectron]>0.5 && q==-1;
else return 0;
}
//=============================================================================
-Bool_t AliUnicorEventAliceESD::PairGood(Double_t /*p0*/, Double_t the0, Double_t phi0,
- Double_t /*p1*/, Double_t the1, Double_t phi1) const {
+Bool_t AliUnicorEventAliceESD::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 dthe = the1-the0;
+ // double dphi = TVector2::Phi_mpi_pi(phi1-phi0);
+ // double dpt = p1*sin(the1) - p0*sin(the0);
+ // return (fabs(dthe)>0.010 || fabs(dphi)>0.060);
+ // return (dpt*dphi<0);
return 1;
- double dthe = the1-the0;
- double dphi = TVector2::Phi_mpi_pi(phi1-phi0);
- // double dpt = p1*sin(the1) - p0*sin(the0);
- return (fabs(dthe)>0.010 || fabs(dphi)>0.060);
- //return (dpt*dphi<0);
}
//=============================================================================
#include <cmath>
#include "TVector2.h"
#include "AliESDEvent.h"
-#include "AliPhysicsSelection.h"
+//#include "AliPhysicsSelection.h"
#include "AliExternalTrackParam.h"
#include "AliUnicorEvent.h"
-class AliPhysicsSelection;
+//class AliPhysicsSelection;
//=============================================================================
class AliUnicorEventAliceESD : public AliUnicorEvent {
public:
AliUnicorEventAliceESD(AliESDEvent *esd=0);
- AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fESD(ev.fESD), fPhysicsSelection(ev.fPhysicsSelection) {}
+ //AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fESD(ev.fESD), fPhysicsSelection(ev.fPhysicsSelection){}
+ AliUnicorEventAliceESD(const AliUnicorEventAliceESD &ev): AliUnicorEvent(ev), fESD(ev.fESD) {}
virtual ~AliUnicorEventAliceESD();
AliUnicorEventAliceESD &operator=(const AliUnicorEventAliceESD &source) {fESD = source.fESD; return *this;}
Double_t Etamin() const {return -0.75;}
Double_t Etamax() const {return 0.75;}
void AttachTree(TTree *tr) {fESD->ReadFromTree(tr);}
Bool_t Good() const;
- Double_t Centrality() const {return 0.9999*exp(-NGoodParticles()/7.0);} // OK for pp
+ Double_t Centrality() const {return 0.9999*exp(-NGoodParticles()/2000.0);} // 7 for pp 900 GeV, 12 for pp 7 TeV, 20 extr, 2000 PbPb
void RP(Double_t &qx, Double_t &qy) const {AliUnicorEvent::RP(qx,qy,2);}
Double_t RPphi() const {Double_t qx,qy; RP(qx,qy); return atan2(qy,qx);}
Double_t Zver() const {return fESD->GetPrimaryVertex()->GetZv()/10.0;}
Int_t NParticles() const {return fESD->GetNumberOfTracks();}
- Int_t NGoodParticles() const {int n=0; for (int i=0; i<NParticles(); i++) if (ParticleGood(i)) n++; return n;}
// Int_t NGoodParticles() const {int n=0; for (int i=0; i<fESD->GetMultiplicity()->GetNumberOfTracklets(); i++) if (fabs(fESD->GetMultiplicity()->GetEta(i))<0.8) n++; return n;}
Bool_t ParticleGood(Int_t i, Int_t pidi=0) const;
void SetESD(AliESDEvent * const esd) {fESD = esd;}
AliESDEvent *GetESD() const {return fESD;}
//const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i);}
- const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetConstrainedParam();}
+ //const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetConstrainedParam();}
//const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetInnerParam();} // not at vtx!
- //const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetTPCInnerParam();}
+ const AliExternalTrackParam *GetTrackParam(Int_t i) const {return fESD->GetTrack(i)->GetTPCInnerParam();}
protected:
AliESDEvent *fESD; //! pointer to the actual source of data
- AliPhysicsSelection *fPhysicsSelection; //! interaction event filter
+ // AliPhysicsSelection *fPhysicsSelection; //! interaction event filter
ClassDef(AliUnicorEventAliceESD,0)
};
ClassImp(AliUnicorHN)
//=============================================================================
-AliUnicorHN::AliUnicorHN(Char_t *nam, Int_t ndim, TAxis **ax)
+AliUnicorHN::AliUnicorHN(const char *nam, Int_t ndim, TAxis **ax)
: TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5),
- fNdim(ndim)
+ fNdim(ndim)
{
// Constructor for building from scratch.
// 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=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
+ for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
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());
}
//=============================================================================
-AliUnicorHN::AliUnicorHN(Char_t *filnam, Char_t *nam) :
+AliUnicorHN::AliUnicorHN(const char *filnam, const char *nam) :
TH1D(*(TFile::Open(filnam,"read")?
TFile::Open(filnam,"read")->GetDirectory(nam)?
(TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"):new TH1D():new TH1D()
}
f->Close();
- fMbins[fNdim-1] = 1;
- for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
+ for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
+ for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
if (GetNbinsX()!=Albins(fNdim,ax)) {
// abscissa of the n-dim histogram. The result will be between 0 and
// GetNbinsX()-1.
- Int_t k[fgkMaxNdim];
+ Int_t k[fgkMaxNdim] = {0};
for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1;
return MulToOne(k);
}
return nbytes;
}
//=============================================================================
-AliUnicorHN *AliUnicorHN::ProjectAlong(char *nam, Int_t dim, Int_t first, Int_t last)
+AliUnicorHN *AliUnicorHN::ProjectAlong(const 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.
// create new (reduced) histogram
- TAxis *ax[fgkMaxNdim-1];
+ TAxis *ax[fgkMaxNdim-1] = {0};
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);
+ const char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
+ const char *eame = Form("%s_error",name);
AliUnicorHN *his = new AliUnicorHN(name,fNdim-1,ax); // result histogram
AliUnicorHN *eis = new AliUnicorHN(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
+ int k[fgkMaxNdim] = {0}; // old hist multiindex
+ int m[fgkMaxNdim] = {0}; // new hist multiindex
for (int i=0; i<GetNbinsX(); i++) {
OneToMul(i,k);
if (k[dim]+1<first) continue;
return his;
}
//=============================================================================
-TH1D *AliUnicorHN::ProjectOn(char *nam, Int_t dim, const Int_t * const first, const Int_t * const last) const
+TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first, const Int_t * const last) const
{
// 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.
// use name nam if specified; otherwise generate one
TH1D *his;
- char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
+ const char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
// if (gDirectory->Get(name)) gDirectory->Get(name)->Delete(); // for some reason leads to troubles
if (fAxis[dim].IsVariableBinSize())
his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
return his;
}
//=============================================================================
-TH1D *AliUnicorHN::ProjectOn(char *nam, Int_t dim, const Double_t * const first, const Double_t * const last)
+TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first, const Double_t * const 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[fgkMaxNdim];
- Int_t klast[fgkMaxNdim];
+ Int_t kfirst[fgkMaxNdim] = {0};
+ Int_t klast[fgkMaxNdim] = {0};
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 *AliUnicorHN::ProjectOn(char *nam, Int_t dim0, Int_t dim1, const Int_t * const first, const Int_t * const last) const
+TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * const first, const Int_t * const last) const
{
// 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.
// use name nam if specified; otherwise generate one
TH2D *his=0;
- char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0);
+ const 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(),
class AliUnicorHN : public TH1D {
public:
- AliUnicorHN() : TH1D(), fNdim(0) {} // default contructor
- AliUnicorHN(Char_t *nam, Int_t ndim, TAxis **ax); // constructor from scratch
- AliUnicorHN(Char_t *filename, Char_t *name); // constructor from file
- virtual ~AliUnicorHN() {} // destructor
+ AliUnicorHN(const char *nam="muhi", Int_t ndim=0, TAxis **ax=0); // constructor from scratch
+ AliUnicorHN(const char *filename, const char *name); // constructor from file
+ virtual ~AliUnicorHN() {} // destructor
Int_t GetNdim() const {return fNdim;}
TAxis *GetAxis(Int_t i) const {return (TAxis*) &fAxis[i];}
Int_t Fill(Double_t *xx, Double_t y=1); // fill histo
Int_t Fill(Double_t) {return -1;} // insufficient number of arguments
- Int_t Fill(Double_t x0, Double_t w) {return Fill(&x0,w);} // 1-dim histo fill
+ Int_t Fill(Double_t x0, Double_t w) {Double_t x[1]={x0}; return Fill(x,w);} // 1-dim histo fill
Int_t Fill(Double_t x0, Double_t x1, ...);// 2 or more dim histo fill
Int_t Fill(const char*, Double_t) {return -1;} // overload TH1
Int_t Write(const char *, Int_t, Int_t) const {return Write();}
// project along (integrate over) one axis
- AliUnicorHN *ProjectAlong(char *nam, Int_t dim, Int_t first=0, Int_t last=0);
+ AliUnicorHN *ProjectAlong(const 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, const Int_t * const first=0, const Int_t * const last=0) const;
+ TH1D *ProjectOn(const char *nam, Int_t dim, const Int_t * const first=0, const Int_t * const last=0) const;
// project on 1-dim histogram
- TH1D *ProjectOn(char *nam, Int_t dim, const Double_t * const first, const Double_t * const last);
+ TH1D *ProjectOn(const char *nam, Int_t dim, const Double_t * const first, const Double_t * const last);
// project on 2-dim histogram
- TH2D *ProjectOn(char *nam, Int_t dim0, Int_t dim1, const Int_t * const first=0, const Int_t * const last=0) const;
+ TH2D *ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t * const first=0, const Int_t * const last=0) const;
void OneToMul(Int_t n, Int_t *k) const; // calc n-dim indices from 1-dim index
protected:
void Set1(double m,double p,double theta,double phi) {fP1.SetXYZM(0,0,p,m); fP1.SetTheta(theta); fP1.SetPhi(phi);}
void SetMXYZ0(double m,double px,double py,double pz) {fP0.SetXYZM(px,py,pz,m);}
void SetMXYZ1(double m,double px,double py,double pz) {fP1.SetXYZM(px,py,pz,m);}
- void CalcLAB() {fP=fP0+fP1; fQ=fP1-fP0; fBeta=fP.BoostVector();
+ void CalcLAB() {fCMp=fP=fP0+fP1; fCMq=fQ=fP1-fP0; fBeta=fP.BoostVector();
fBetaz.SetXYZ(0,0,fBeta.Z()); fBetat=fBeta; fBetat.SetZ(0);
fUbeta=fBeta.Unit(); fUbetat=fBetat.Unit(); fUbetaz=fBetaz.Unit();}
double Rapidity() const {return fP.Rapidity();}
double Phi() const {return fP.Phi();}
double DTheta() const {return fP1.Theta()-fP0.Theta();}
double DPhi() const {return TVector2::Phi_mpi_pi(fP1.Phi()-fP0.Phi());}
- void CalcPairCM() {fCMp=fP; fCMp.Boost(-fBeta); fCMq=fQ; fCMq.Boost(-fBeta);}
- void CalcLcmsCM() {fCMp=fP; fCMp.Boost(-fBetaz); fCMq=fQ; fCMq.Boost(-fBetaz);}
+ void CalcPairCM() {fCMp.Boost(-fBeta); fCMq.Boost(-fBeta);}
+ void CalcLcmsCM() {fCMp.Boost(-fBetaz); fCMq.Boost(-fBetaz);}
void Swap() {fBuf=fP0; fP0=fP1; fP1=fBuf; fQ=-fQ; fCMq=-fCMq;}
double Minv() const {return fP.M();}
double Qinv2() const {return -fQ.M2();}
-{
-gSystem->Load("libPhysics.so");
-gSystem->Load("libEG.so");
-gSystem->Load("libTree.so");
-gSystem->Load("libVMC.so");
-gSystem->Load("libSTEERBase.so");
-gSystem->Load("libESD.so");
-gSystem->Load("libAOD.so");
-gSystem->Load("libANALYSIS");
-gSystem->Load("libANALYSISalice");
-gSystem->Load("libPWG2unicor");
-
-gROOT->LoadMacro("makechain.C");
-tr = makechain("esdTree","filelist.txt");
-
-AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
-AliVEventHandler* esdH = new AliESDInputHandler;
-mgr->SetInputEventHandler(esdH);
+//=============================================================================
+void runAsTask(int grid=0) {
+ gSystem->Load("libPhysics.so");
+ gSystem->Load("libEG.so");
+ gSystem->Load("libTree.so");
+ gSystem->Load("libVMC.so");
+ gSystem->Load("libSTEERBase.so");
+ gSystem->Load("libESD.so");
+ gSystem->Load("libAOD.so");
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+ gSystem->Load("libPWG2unicor");
+
+ if (grid) {
+ if (!TGrid::Connect("alien://")) return;
+ TChain *tr = CreateChainFromTags("wn.xml", "ESD");
+ } else {
+ gROOT->LoadMacro("makechain.C");
+ tr = makechain("esdTree","filelist.txt");
+ }
-gROOT->LoadMacro("AddTaskUnicor.C");
+ AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+ AliVEventHandler* esdH = new AliESDInputHandler;
+ mgr->SetInputEventHandler(esdH);
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/UNICOR/AddTaskUnicor.C");
+ AliAnalysisTaskUnicor *mytask = AddTaskUnicor();
-AliAnalysisTaskUnicor *mytask = AddTaskUnicor();
+ mgr->InitAnalysis();
+ mgr->PrintStatus();
+ mgr->StartAnalysis("local",tr);
-mgr->InitAnalysis();
-mgr->PrintStatus();
-mgr->StartAnalysis("local",tr);
-
-TFile::Open("AnalysisResults.root","read");
-gDirectory->Cd("PWG2UNICOR");
-TList *list = (TList *) gDirectory->Get("unilis");
-char *outfil = "unicor-result-as-anal.root";
-for (int i=0; i<list->GetEntries(); i++) {
- AliUnicorAnal *an = (AliUnicorAnal *) list->At(i);
- if (i==0) an->Save(outfil,"recreate");
- else an->Save(outfil);
- delete an;
+ TFile::Open("AnalysisResults.root","read");
+ gDirectory->Cd("PWG2UNICOR");
+ TList *list = (TList *) gDirectory->Get("unilis");
+ char *outfil = "unicor-result-as-anal.root";
+ for (int i=0; i<list->GetEntries(); i++) {
+ AliUnicorAnal *an = (AliUnicorAnal *) list->At(i);
+ if (i==0) an->Save(outfil,"recreate");
+ else an->Save(outfil);
+ delete an;
+ }
}
+//=============================================================================
+TChain* CreateChainFromTags(const char *xmlfile, const char *type="ESD")
+{
+ // Create a chain using tags from the xml file.
+ TAlienCollection* coll = TAlienCollection::Open(xmlfile);
+ if (!coll) {
+ ::Error("CreateChainFromTags", "Cannot create an AliEn collection from %s", xmlfile);
+ return NULL;
+ }
+ TGridResult* tagResult = coll->GetGridResult("",kFALSE,kFALSE);
+ AliTagAnalysis *tagAna = new AliTagAnalysis(type);
+ tagAna->ChainGridTags(tagResult);
+
+ AliRunTagCuts *runCuts = new AliRunTagCuts();
+ AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
+ AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
+ AliEventTagCuts *evCuts = new AliEventTagCuts();
+ // Check if the cuts configuration file was provided
+ if (!gSystem->AccessPathName("ConfigureCuts.C")) {
+ gROOT->LoadMacro("ConfigureCuts.C");
+ ConfigureCuts(runCuts, lhcCuts, detCuts, evCuts);
+ }
+ TChain *chain = tagAna->QueryTags(runCuts, lhcCuts, detCuts, evCuts);
+ if (!chain || !chain->GetNtrees()) return NULL;
+ chain->ls();
+ return chain;
}
+//=============================================================================