#include "AliGeomManager.h"
#include "AliInputEventHandler.h"
#include "AliLog.h"
+#include "AliMCEvent.h"
+#include "AliMCParticle.h"
#include "AliMagF.h"
#include "AliMultiplicity.h"
+#include "AliStack.h"
#include "AliTrackerBase.h"
ClassImp(AliAnalysisTaskEMCALPi0PbPb)
fMarkCells(),
fMinL0Time(-1),
fMaxL0Time(1024),
- fIsGeoMatsSet(0),
- fNEvs(0),
+ fMcMode(0),
fGeom(0),
fReco(0),
+ fIsGeoMatsSet(0),
+ fNEvs(0),
fOutput(0),
fTrClassNamesArr(0),
fEsdEv(0),
cout << " fMinEcc: " << fMinEcc << endl;
cout << " fGeoName: \"" << fGeoName << "\"" << endl;
cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
+ cout << " fIsoDist: " << fIsoDist << endl;
cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
cout << " fTrCuts: " << fTrCuts << endl;
cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
cout << " fMarkCells: " << fMarkCells << endl;
cout << " fMinL0Time: " << fMinL0Time << endl;
cout << " fMaxL0Time: " << fMaxL0Time << endl;
-
- fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
- fReco = new AliEMCALRecoUtils();
+ cout << " fMcMode: " << fMcMode << endl;
+ cout << " fGeom: " << fGeom << endl;
+ cout << " fReco: " << fReco << endl;
+
+ if (!fGeom)
+ fGeom = new AliEMCALGeoUtils(fGeoName,"EMCAL");
+ else {
+ if (fGeom->GetMatrixForSuperModule(0))
+ fIsGeoMatsSet = kTRUE;
+ }
+ if (!fReco)
+ fReco = new AliEMCALRecoUtils();
fTrClassNamesArr = fTrClassNames.Tokenize(" ");
fOutput = new TList();
fOutput->SetOwner();
am->LoadBranch("header");
offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
}
- if (offtrigger & AliVEvent::kFastOnly) {
+ if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
AliWarning(Form("EMCAL not in fast only partition"));
return;
}
recalibrated = kTRUE;
}
}
- if (1 && AODEvent() && !fClusName.IsNull()) {
- TList *l = AODEvent()->GetList();
- TClonesArray *clus = 0;
+ if (1 && !fClusName.IsNull()) {
+ TList *l = 0;
+ if (AODEvent())
+ l = AODEvent()->GetList();
+ else if (fAodEv)
+ l = fAodEv->GetList();
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
- fAodClusters = clus;
+ fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
}
}
if (!clusattached)
am->LoadBranch("CaloClusters");
TList *l = fEsdEv->GetList();
- TClonesArray *clus = 0;
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
- fEsdClusters = clus;
+ fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
}
}
if (1) {
if (!clusattached)
am->LoadBranch("caloClusters");
TList *l = fAodEv->GetList();
- TClonesArray *clus = 0;
if (l) {
- clus = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
- fAodClusters = clus;
+ fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
}
}
if (1) {
if (fDoAfterburner)
ClusterAfterburner();
+ CalcMcInfo();
CalcCaloTriggers();
CalcPrimTracks();
CalcTracks();
FillPionHists();
FillOtherHists();
}
+ FillMcHists();
FillNtuple();
if (fTrainMode) {
{
// Calculate triggers
+ if (fAodEv)
+ return; // information not available in AOD
+
fTriggers->Clear();
AliVCaloCells *cells = fEsdCells;
cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
- Double_t trigpen = 0;
- Double_t trignen = 0;
- for(Int_t j=0; j<cl->fN; ++j) {
- Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
- Short_t pos = -1;
- std::map<Short_t,Short_t>::iterator it = map.find(cid);
- if (it!=map.end())
- pos = it->second;
- if (pos<0)
- continue;
- if (fAmpInTrigger[pos]>0)
- trigpen += cells->GetAmplitude(pos);
- else if (fAmpInTrigger[pos]<0)
- trignen += cells->GetAmplitude(pos);
- }
- if (trigpen>0) {
- cl->fIsTrigM = 1;
- cl->fTrigE = trigpen;
- }
- if (trignen>0) {
- cl->fIsTrigM = 1;
- cl->fTrigMaskE = trignen;
+ if (fAmpInTrigger) { // fill trigger info if present
+ Double_t trigpen = 0;
+ Double_t trignen = 0;
+ for(Int_t j=0; j<cl->fN; ++j) {
+ Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
+ Short_t pos = -1;
+ std::map<Short_t,Short_t>::iterator it = map.find(cid);
+ if (it!=map.end())
+ pos = it->second;
+ if (pos<0)
+ continue;
+ if (fAmpInTrigger[pos]>0)
+ trigpen += cells->GetAmplitude(pos);
+ else if (fAmpInTrigger[pos]<0)
+ trignen += cells->GetAmplitude(pos);
+ }
+ if (trigpen>0) {
+ cl->fIsTrigM = 1;
+ cl->fTrigE = trigpen;
+ }
+ if (trignen>0) {
+ cl->fIsTrigM = 1;
+ cl->fTrigMaskE = trignen;
+ }
}
cl->fIsShared = IsShared(clus);
}
Double_t trkPos[3] = {0};
tParam->GetXYZ(trkPos); //Get the extrapolated global position
- tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2)+TMath::Power(clsPos[1]-trkPos[1],2)+TMath::Power(clsPos[2]-trkPos[2],2) );
+ tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2) +
+ TMath::Power(clsPos[1]-trkPos[1],2) +
+ TMath::Power(clsPos[2]-trkPos[2],2) );
tmpZ = clsPos[2]-trkPos[2];
delete tParam;
} else {
}
}
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
+{
+ // Get Mc truth particle information.
+
+ if (!fMcMode)
+ return;
+
+// todo: treat mc differently for AOD and ESD. Check with Amorsch when we get MC handler also for AOD case.
+
+ AliMCEvent *mcEvent = MCEvent();
+ if (!mcEvent)
+ return;
+
+ const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
+ if (!evtVtx)
+ return;
+
+ mcEvent->PreReadAll();
+ AliStack *stack = mcEvent->Stack();
+
+ Int_t nTracks = mcEvent->GetNumberOfPrimaries();
+ for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
+ AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
+ if (!mcP)
+ continue;
+
+ // pion or eta meson
+ TParticle *tP = mcP->Particle();
+ if(tP->GetPdgCode() == 111) {
+ } else if(tP->GetPdgCode() == 221) {
+ } else
+ continue;
+
+ // primary particle
+ Double_t dR = TMath::Sqrt((tP->Vx()-evtVtx->GetX())*(tP->Vx()-evtVtx->GetX())+(tP->Vy()-evtVtx->GetY())*(tP->Vy()-evtVtx->GetY()));
+ if(dR > 1.0)
+ continue;
+
+ // do not account dalitz decays
+ if(tP->GetNDaughters()!=2)
+ continue;
+
+ // kinematic cuts
+ Double_t pt = tP->Pt() ;
+ if (pt<0.5) {
+ continue;
+ }
+ Double_t rap = tP->Y();
+ if (TMath::Abs(rap)>1){
+ continue;
+ }
+
+ // get daughters
+ TParticle *tD1 = 0;
+ TClonesArray *tD1refs = 0;
+ mcEvent->GetParticleAndTR(tP->GetFirstDaughter(),tD1,tD1refs);
+
+ TParticle *tD2 = 0;
+ TClonesArray *tD2refs = 0;
+ mcEvent->GetParticleAndTR(tP->GetLastDaughter(),tD2,tD2refs);
+
+ TParticle *tD11 = 0;
+ TClonesArray *tD11refs = 0;
+ TParticle *tD12 = 0;
+ TClonesArray *tD12refs = 0;
+ Bool_t d1Dec = 0;
+ if (tD1->GetNDaughters()>2) {
+ AliWarning(Form("Decay photon has more than two daughters: %d", tD1->GetNDaughters()));
+ return;
+ }
+ if (tD1->GetNDaughters()>=1)
+ mcEvent->GetParticleAndTR(tD1->GetFirstDaughter(),tD11,tD11refs);
+ if (tD1->GetNDaughters()==2) {
+ mcEvent->GetParticleAndTR(tD1->GetLastDaughter(),tD12,tD12refs);
+ d1Dec = 1;
+ }
+
+ TParticle *tD21 = 0;
+ TClonesArray *tD21refs = 0;
+ TParticle *tD22 = 0;
+ TClonesArray *tD22refs = 0;
+ Bool_t d2Dec = 0;
+ if (tD2->GetNDaughters()>2) {
+ AliWarning(Form("Decay photon has more than two daughters: %d", tD2->GetNDaughters()));
+ return;
+ }
+ if (tD2->GetNDaughters()>=1)
+ mcEvent->GetParticleAndTR(tD2->GetFirstDaughter(),tD21,tD21refs);
+ if (tD2->GetNDaughters()==2) {
+ mcEvent->GetParticleAndTR(tD2->GetLastDaughter(),tD22,tD22refs);
+ d2Dec = 1;
+ }
+
+ Bool_t pDec = d1Dec && d2Dec;
+ if (!pDec)
+ continue;
+
+ if (tD11) {
+ cout << " Dau 11 with " << tD11->GetNDaughters() << endl;
+ tD11->Print();
+ }
+
+ if (0) {
+ tP->Print();
+ tD1->Print();
+ tD2->Print();
+ if (tD11) {
+ cout << " Dau 11 with " << tD11->GetNDaughters() << endl;
+ tD11->Print();
+ if (tD11refs) {
+ Int_t nrefs = tD11refs->GetEntries();
+ for (Int_t j=0;j<nrefs;++j) {
+ AliTrackReference *tr = (AliTrackReference*)tD11refs->At(j);
+ tr->SetUserId(tr->DetectorId());
+ tr->Print();
+ }
+ }
+ }
+ if (tD12){
+ cout << " Dau 12 with " << tD12->GetNDaughters() << endl;
+ tD12->Print();
+ if (tD12refs) {
+ Int_t nrefs = tD12refs->GetEntries();
+ for (Int_t j=0;j<nrefs;++j) {
+ AliTrackReference *tr = (AliTrackReference*)tD12refs->At(j);
+ tr->SetUserId(tr->DetectorId());
+ tr->Print();
+ }
+ }
+ }
+ if (tD21) {
+ cout << " Dau 21 with " << tD21->GetNDaughters() << endl;
+ tD21->Print();
+ if (tD21refs) {
+ Int_t nrefs = tD21refs->GetEntries();
+ for (Int_t j=0;j<nrefs;++j) {
+ AliTrackReference *tr = (AliTrackReference*)tD21refs->At(j);
+ tr->SetUserId(tr->DetectorId());
+ tr->Print();
+ }
+ }
+ }
+ if (tD22) {
+ cout << " Dau 22 with " << tD22->GetNDaughters() << endl;
+ tD22->Print();
+ if (tD22refs) {
+ Int_t nrefs = tD22refs->GetEntries();
+ for (Int_t j=0;j<nrefs;++j) {
+ AliTrackReference *tr = (AliTrackReference*)tD22refs->At(j);
+ tr->SetUserId(tr->DetectorId());
+ tr->Print();
+ }
+ }
+ }
+ }
+#if 0
+ TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
+ TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
+ //Number of pi0s decayed into acceptance
+ Bool_t inAcc1 = (TMath::Abs(gamma1->Eta())<0.9) ;
+ Bool_t inAcc2 = (TMath::Abs(gamma2->Eta())<0.9) ;
+ Int_t mod1,mod2 ;
+ Double_t x=0.,z=0. ;
+ Bool_t hitPHOS1 = fPHOSgeom->ImpactOnEmc(gamma1, mod1, z,x) ;
+ Bool_t hitPHOS2 = fPHOSgeom->ImpactOnEmc(gamma2, mod2, z,x) ;
+ Bool_t hitEMCAL1= fEMCALgeom->Impact(gamma1) ;
+ Bool_t hitEMCAL2= fEMCALgeom->Impact(gamma2) ;
+
+
+ TParticle *mcP2 = 0;
+ TClonesArray *refs = 0;
+ mcEvent->GetParticleAndTR(iTrack,mcP2,refs);
+ if (!refs)
+ continue;
+
+ Int_t nrefs = refs->GetEntries();
+ if (nrefs == 0)
+ continue;
+
+ AliTrackReference *trEmcal = 0;
+ for (Int_t j=0;j<nrefs;++j) {
+ AliTrackReference *tr = (AliTrackReference*)refs->At(j);
+ Int_t did = tr->DetectorId();
+ if (did!=AliTrackReference::kEMCAL)
+ continue;
+ trEmcal = tr;
+ break;
+ }
+
+ if (!trEmcal)
+ continue;
+
+ Double_t distA = TMath::Sqrt((trEmcal->X()-tP->Vx())*(trEmcal->X()-tP->Vx()) +
+ (trEmcal->Y()-tP->Vy())*(trEmcal->Y()-tP->Vy()) +
+ (trEmcal->Z()-tP->Vz())*(trEmcal->Z()-tP->Vz()));
+ cout << "dist label to particle " << distA << " --------- " << endl;
+ trEmcal->Print();
+ tP->Print();
+
+ Int_t moi = tP->GetFirstMother();
+ if (moi>=0) {
+ AliMCParticle *mcMo = static_cast<AliMCParticle*>(mcEvent->GetTrack(moi));
+ TParticle *tMo = mcMo->Particle();
+ Double_t dist3d = TMath::Sqrt((tMo->Vx()-tP->Vx())*(tMo->Vx()-tP->Vx()) +
+ (tMo->Vy()-tP->Vy())*(tMo->Vy()-tP->Vy()) +
+ (tMo->Vz()-tP->Vz())*(tMo->Vz()-tP->Vz()));
+ cout << "dist 3d " << dist3d << endl;
+ tMo->Print();
+ while (1) {
+ Int_t newm = tMo->GetFirstMother();
+ if (newm>=0) {
+ AliMCParticle *newmo = static_cast<AliMCParticle*>(mcEvent->GetTrack(newm));
+ tMo = newmo->Particle();
+ tMo->Print();
+ } else {
+ break;
+ }
+ }
+ cout << " ------- " << endl;
+ }
+#endif
+ }
+}
+
//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
{
}
}
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
+{
+ // Fill additional MC information histograms.
+
+ if (!fMcMode)
+ return;
+
+ // check if aod or esd mc mode and the fill histos
+}
+
//________________________________________________________________________
void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
{
return 0;
}
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(AliMCParticle *p) const
+{
+ // Print daugher information. TODO
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
+{
+ // Print track ref array. TODO
+}
+
class AliESDTrack;
class AliESDVertex;
class AliESDtrackCuts;
+class AliMCParticle;
class AliStaHeader;
class AliStaVertex;
void SetDoTrackMatWithGeom(Bool_t b) { fDoTrMatGeom = b; }
void SetFillNtuple(Bool_t b) { fDoNtuple = b; }
void SetGeoName(const char *n) { fGeoName = n; }
+ void SetGeoUtils(AliEMCALGeoUtils *geo) { fGeom = geo; }
void SetIsoDist(Double_t d) { fIsoDist = d; }
+ void SetL0TimeRange(Int_t l, Int_t h) { fMinL0Time=l; fMaxL0Time=h; }
void SetMarkCells(const char *n) { fMarkCells = n; }
+ void SetMcMode(Bool_t b) { fMcMode = b; }
void SetMinClusEnergy(Double_t e) { fMinE = e; }
void SetMinEcc(Double_t ecc) { fMinEcc = ecc; }
void SetMinErat(Double_t erat) { fMinErat = erat; }
void SetMinNClustersPerTrack(Double_t m) { fMinNClusPerTr = m; }
void SetNminCells(Int_t n) { fNminCells = n; }
void SetPrimTrackCuts(AliESDtrackCuts *c) { fPrimTrCuts = c; }
+ void SetRecoUtils(AliEMCALRecoUtils *reco) { fReco = reco; }
void SetTrClassNames(const char *n) { fTrClassNames = n; }
void SetTrackCuts(AliESDtrackCuts *c) { fTrCuts = c; }
void SetTrainMode(Bool_t b) { fTrainMode = b; }
void SetUseQualFlag(Bool_t b) { fUseQualFlag = b; }
void SetVertexRange(Double_t z1, Double_t z2) { fVtxZMin=z1; fVtxZMax=z2; }
- void SetL0TimeRange(Int_t l, Int_t h) { fMinL0Time=l; fMaxL0Time=h; }
protected:
- void CalcCaloTriggers();
+ virtual void CalcCaloTriggers();
virtual void CalcClusterProps();
virtual void CalcPrimTracks();
+ virtual void CalcMcInfo();
virtual void CalcTracks();
virtual void ClusterAfterburner();
virtual void FillCellHists();
virtual void FillNtuple();
virtual void FillOtherHists();
virtual void FillPionHists();
+ virtual void FillMcHists();
void FillVertex(AliStaVertex *v, const AliESDVertex *esdv);
void FillVertex(AliStaVertex *v, const AliAODVertex *aodv);
Double_t GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2) const;
Double_t GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius=0.2, Double_t pt=0.) const;
Double_t GetTrigEnergy(const AliVCluster *c) const;
Bool_t IsShared(const AliVCluster *c) const;
+ void PrintDaughters(AliMCParticle *p) const;
+ void PrintTrackRefs(AliMCParticle *p) const;
// input members
TString fCentVar; // variable for centrality determination
TString fMarkCells; // list of mark cells to monitor
Int_t fMinL0Time; // minimum accepted time for trigger
Int_t fMaxL0Time; // maximum accepted time for trigger
+ Bool_t fMcMode; // monte carlo mode
+ AliEMCALGeoUtils *fGeom; // geometry utils
+ AliEMCALRecoUtils *fReco; // reco utils
+
// derived members (ie with ! after //)
Bool_t fIsGeoMatsSet; //!indicate that geo matrices are set
ULong64_t fNEvs; //!accepted events
- AliEMCALGeoUtils *fGeom; //!geometry utils
- AliEMCALRecoUtils *fReco; //!geometry utils
TList *fOutput; //!container of output histograms
TObjArray *fTrClassNamesArr; //!array of trig class names
AliESDEvent *fEsdEv; //!pointer to input esd event
TH2 *fHPionMggAsym; //!histo for pion mass vs. asym
TH2 *fHPionMggDgg; //!histo for pion mass vs. opening angle
TH1 *fHPionInvMasses[21]; //!histos for invariant mass plots
+ // histograms for MC
private:
AliAnalysisTaskEMCALPi0PbPb(const AliAnalysisTaskEMCALPi0PbPb&); // not implemented
AliAnalysisTaskEMCALPi0PbPb &operator=(const AliAnalysisTaskEMCALPi0PbPb&); // not implemented
- ClassDef(AliAnalysisTaskEMCALPi0PbPb, 7) // Analysis task for neutral pions in Pb+Pb
+ ClassDef(AliAnalysisTaskEMCALPi0PbPb, 9) // Analysis task for neutral pions in Pb+Pb
};
#endif
public:
Double32_t fE; //[0,0,16] energy
- Double32_t fR; //[0,0,16] radius
+ Double32_t fR; //[0,0,16] radius (cylinder)
Double32_t fEta; //[0,0,16] eta
Double32_t fPhi; //[0,0,16] phi
UChar_t fN; // number of cells
ClassDef(AliStaCluster,4) // Cluster class
};
+class AliStaTrackRef : public TObject
+{
+ public:
+ AliStaTrackRef() : TObject(), fR(0), fEta(0), fPhi(0), fId(-1), fMo(-1) {}
+
+ public:
+ Double32_t fR; //[0,0,16] r (cylinder)
+ Double32_t fEta; //[0,0,16] eta
+ Double32_t fPhi; //[0,0,16] phi
+ Short_t fId; // id
+ Short_t fMo; // mother index
+
+ ClassDef(AliStaTrackRef,1) // Track reference class
+};
+
class AliStaTrigger : public TObject
{
public:
ClassDef(AliStaTrigger,1) // Trigger class
};
+
+
+#if 0
+class AliStaPart : public TObject
+{
+ public:
+ AliStaTrigger() : TObject(), fE(0), fEta(0), fPhi(0), fAmp(0), fMinTime(0), fMaxTime(0) {}
+
+ public:
+ Double32_t fE; //[0,0,16] pt
+ Double32_t fEta; //[0,0,16] eta
+ Double32_t fPhi; //[0,0,16] phi
+ Int_t fId; //[0,0,16] id
+
+ ClassDef(AliStaTrigger,1) // Trigger class
+};
+#endif
+
#endif