+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-//
-// AliCDMesonUtils
-// for
-// AliAnalysisTaskCDMeson
-//
-// Author:
-// Xianguo Lu <lu@physi.uni-heidelberg.de>
-// continued by
-// Felix Reidt <Felix.Reidt@cern.ch>
-
-#include <TH1.h>
-#include <TH2.h>
-#include <TGeoMatrix.h>
-#include <THnSparse.h>
-#include <TLorentzVector.h>
-#include <TRandom3.h>
-#include <TParticle.h>
-#include <TObjArray.h>
-
-#include "AliCDBManager.h"
-#include "AliCDBEntry.h"
-#include "AliCDBStorage.h"
-#include "AliESDEvent.h"
-#include "AliPIDResponse.h"
-#include "AliVTrack.h"
-#include "AliVParticle.h"
-#include "AliESDtrack.h"
-#include "AliESDtrackCuts.h"
-#include "AliESDFMD.h"
-#include "AliESDVZERO.h"
-#include "AliGeomManager.h"
-#include "AliITSAlignMille2Module.h"
-#include "AliITSsegmentationSPD.h"
-#include "AliMultiplicity.h"
-#include "AliPIDResponse.h"
-#include "AliSPDUtils.h"
-#include "AliTriggerAnalysis.h"
-#include "AliVZEROTriggerData.h"
-#include "AliVZEROCalibData.h"
-
-#include "AliAODTracklets.h"
-#include "AliAODEvent.h"
-
-#include "AliCDMesonBase.h"
-#include "AliCDMesonTracks.h"
-
-#include "AliCDMesonUtils.h"
-
-
-//==============================================================================
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetMassPtCtsOA(Int_t pid,
- const TVector3* momenta[], Double_t & mass,
- Double_t &pt, Double_t &cts, Double_t &oa)
-{
- //
- // Get Mass Pt Cts OA
- //
-
- Double_t tmpp1[3], tmpp2[3];
- momenta[0]->GetXYZ(tmpp1);
- momenta[1]->GetXYZ(tmpp2);
-
- Double_t masstrk = -999;
- if(pid==AliCDMesonBase::kBinPionE || pid==AliCDMesonBase::kBinPion ||
- pid==AliCDMesonBase::kBinSinglePion)
- masstrk = AliPID::ParticleMass(AliPID::kPion);
- else if(pid==AliCDMesonBase::kBinKaonE || pid==AliCDMesonBase::kBinKaon ||
- pid==AliCDMesonBase::kBinSingleKaon)
- masstrk = AliPID::ParticleMass(AliPID::kKaon);
- else if(pid==AliCDMesonBase::kBinProtonE || pid==AliCDMesonBase::kBinProton ||
- pid==AliCDMesonBase::kBinSingleProton)
- masstrk = AliPID::ParticleMass(AliPID::kProton);
- else if(pid==AliCDMesonBase::kBinElectronE ||
- pid==AliCDMesonBase::kBinElectron ||
- pid==AliCDMesonBase::kBinSingleElectron)
- masstrk = AliPID::ParticleMass(AliPID::kElectron);
- else
- masstrk = AliPID::ParticleMass(AliPID::kPion);
-
- const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masstrk, masstrk,
- cts);
- mass = sumv.M();
- pt = sumv.Pt();
-
- oa = GetOA(tmpp1, tmpp2);
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetPWAinfo(Int_t pid, const AliVTrack *trks[],
- Float_t& theta, Float_t& phi, Float_t& mass,
- Float_t momentum[])
-{
- // transforms the coordiante system to the helicity frame and determines
- // invariant mass and 3-vector of the two track system
- //
-
- Double_t tmpp1[3], tmpp2[3]; // 3-vectors of the daughter tracks
- trks[0]->GetPxPyPz(tmpp1);
- trks[1]->GetPxPyPz(tmpp2);
-
- // determine mass of the daugther tracks
- Double_t masstrk = -999;
- if(pid==AliCDMesonBase::kBinPion || pid==AliCDMesonBase::kBinSinglePion)
- masstrk = AliPID::ParticleMass(AliPID::kPion);
- else if(pid==AliCDMesonBase::kBinKaon || pid==AliCDMesonBase::kBinSingleKaon)
- masstrk = AliPID::ParticleMass(AliPID::kKaon);
- else if(pid==AliCDMesonBase::kBinProton ||
- pid==AliCDMesonBase::kBinSingleProton)
- masstrk = AliPID::ParticleMass(AliPID::kProton);
- else if(pid==AliCDMesonBase::kBinElectron ||
- pid==AliCDMesonBase::kBinSingleElectron)
- masstrk = AliPID::ParticleMass(AliPID::kElectron);
- else
- masstrk = AliPID::ParticleMass(AliPID::kPion);
-
- TLorentzVector va, vb; // 4-vectors of the two tracks
- va.SetXYZM(tmpp1[0], tmpp1[1], tmpp1[2], masstrk);
- vb.SetXYZM(tmpp2[0], tmpp2[1], tmpp2[2], masstrk);
-
- TLorentzVector sumv = va+vb; // 4-vector of the pair
- TVector3 sumXZ(sumv.X(), 0, sumv.Z()); // its projection to the xz-plane
-
- // mass and momentum
- mass = sumv.M();
- momentum[0] = sumv.X();
- momentum[1] = sumv.Y();
- momentum[2] = sumv.Z();
-
- // coordinate axis vectors
- const TVector3 x(1.,0,0);
- const TVector3 y(0,1.,0);
- const TVector3 z(0,0,1.);
-
- const Double_t alpha = -z.Angle(sumXZ);
-
- sumXZ.Rotate(alpha, y);
- sumv.Rotate(alpha, y);
- va.Rotate(alpha, y);
- vb.Rotate(alpha, y);
-
- const Double_t beta = z.Angle(sumv.Vect());
- sumv.Rotate(beta, x);
- va.Rotate(beta, x);
- vb.Rotate(beta, x);
-
-
- const TVector3 bv = -sumv.BoostVector();
- va.Boost(bv);
- vb.Boost(bv);
-
- // angles
- theta = va.Theta();
- phi = va.Phi();
-}
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetEventType(const AliESDEvent *ESDEvent)
-{
- // checks of which type a event is:
- // beam-beam interaction (I), beam-gas (A/C), empty (E)
- //
-
- TString firedTriggerClasses = ESDEvent->GetFiredTriggerClasses();
-
- if (firedTriggerClasses.Contains("CINT1A-ABCE-NOPF-ALL")) { // A
- return AliCDMesonBase::kBinEventA;
- }
- if (firedTriggerClasses.Contains("CINT1C-ABCE-NOPF-ALL")) { // C
- return AliCDMesonBase::kBinEventC;
- }
- if (firedTriggerClasses.Contains("CINT1B-ABCE-NOPF-ALL")) { // I
- return AliCDMesonBase::kBinEventI;
- }
- if (firedTriggerClasses.Contains("CINT1-E-NOPF-ALL")) { // E
- return AliCDMesonBase::kBinEventE;
- }
- if (firedTriggerClasses.Contains("CDG5-E")) { // E
- return AliCDMesonBase::kBinEventE;
- }
- if (firedTriggerClasses.Contains("CDG5-I")) { // I
- return AliCDMesonBase::kBinEventI;
- }
- if (firedTriggerClasses.Contains("CDG5-AC")) { // AC
- return AliCDMesonBase::kBinEventAC;
- }
- return AliCDMesonBase::kBinEventUnknown;
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::SwapTrack(const AliVTrack* trks[])
-{
- //
- // swap two esdtracks, needed since only the properties of one a stored and
- // AliRoot sorts them
- //
-
- TRandom3 tmprd(0);
- if(tmprd.Rndm()>0.5)
- return;
-
- const AliVTrack* tmpt = trks[0];
- trks[0]=trks[1];
- trks[1]=tmpt;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetCombCh(const AliVTrack * trks[])
-{
- //
- // get combination of charges
- //
-
- const Double_t ch1 = trks[0]->Charge();
- const Double_t ch2 = trks[1]->Charge();
-
- if(ch1*ch2<0){
- return AliCDMesonBase::kBinPM;
- }
- else
- return AliCDMesonBase::kBinPPMM;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetCombPID(AliPIDResponse* pid, const AliVTrack* trks[],
- Int_t mode, TH2* comb2trkPID /*= 0x0 */)
-{
- //
- // Get PID for a two track system (data)
- //
-
- if (!pid){
- return AliCDMesonBase::kBinPIDUnknown;
- }
-
- // get PID for the single tracks
- Int_t kpid[2];
- for(Int_t ii=0; ii<2; ii++){
- kpid[ii] = GetPID(pid, trks[ii], mode);
- }
-
- // PID QA histogram
- if (comb2trkPID) {
- comb2trkPID->Fill(kpid[0], kpid[1]);
- }
-
- return CombinePID(kpid);
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetCombPID(const TParticle* particles[],
- TH2 *comb2trkPID /* = 0x0 */)
-{
- //
- // Get PID for a two track system (MC)
- //
-
- // get PID for the single tracks
- Int_t kpid[2];
- for(Int_t ii=0; ii<2; ii++){
- kpid[ii] = GetPID(particles[ii]->GetPdgCode());
- }
-
- // PID QA histogram
- if (comb2trkPID) {
- comb2trkPID->Fill(kpid[0], kpid[1]);
- }
-
- return CombinePID(kpid);
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetSPDTrackletMult(const AliESDEvent *ESDEvent,
- Int_t& sum, Int_t& forwardA,
- Int_t& forwardC, Int_t& central)
-{
- // obtain the multiplicity for eta < -.9 (forwardC), eta > 0.9 (fowardA), in
- // the central barrel (central) and its sum from SPD tracklets with the
- // AliMultiplicity class
- // code simular to PWGPP/ITS/AliAnalysisTaskSPD.cxx
-
- sum = forwardA = forwardC = central = 0; // initialize values
-
- const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
- for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
- iTracklet++) {
- float_t eta = mult->GetEta(iTracklet);
- if (eta < -0.9) {
- forwardC++;
- }
- else if (eta < 0.9) {
- central++;
- }
- else {
- forwardA++;
- }
- sum++;
- }
-}
-
-
-//------------------------------------------------------------------------------
-Bool_t AliCDMesonUtils::CutEvent(const AliESDEvent *ESDEvent, TH1 *hspd,
- TH1 *hpriv, TH1 *hpriVtxPos, TH1 *hpriVtxDist,
- TH2 *hfo, TH1* hfochans, Int_t &kfo,
- Int_t &nip, Int_t &nop, TH1 *hpriVtxX,
- TH1 *hpriVtxY, TH1 *hpriVtxZ)
-{
- //
- // CutEvent
- //
-
- AliTriggerAnalysis triggerAnalysis;
- /*
- //printf("AliCDMesonUtils active triggers: %s\n",
- // ESDEvent->GetHeader()->GetActiveTriggerInputs().Data());
- //printf("AliCDMesonUtils fired triggers: %s\n",
- // ESDEvent->GetHeader()->GetFiredTriggerInputs().Data());
- //*/
- //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
- //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
- // returns the number of fired chips in the SPD
- // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
- // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
-
-
- // SPD FastOR cut
- const Int_t fastORhw = triggerAnalysis.SPDFiredChips(ESDEvent, 1);
- if (hspd) hspd->Fill(fastORhw);
- /*
- if(fastORhw<1)
- return kFALSE;
- */
-
- if (hfochans) {
- const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
-
- for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
- if (mult->TestFastOrFiredChips(iChipKey)) {
- hfochans->Fill((Double_t)iChipKey);
- }
- }
- }
-
- if(kfo){ // spd trigger study
- Int_t nfoctr[10];
- GetNFO(ESDEvent, "[1.0]", nfoctr, 0x0, 0x0);
- nip = nfoctr[kInnerPixel];
- nop = nfoctr[kOuterPixel];
-
- if(AliCDMesonBase::GetGapBin("V0", GetV0(ESDEvent)) ==
- AliCDMesonBase::kBinDG) {
- hfo->Fill(nip, nop);
- }
-
- if(nop<2)
- kfo = 1;
- else
- kfo = nop;
-
- if(kfo>=10)
- kfo=9;
- }
-
- // collision vertex cut
- // A cut in XY is implicitly done during the reconstruction by constraining
- // the vertex to the beam diamond.
-
- // Primary vertex
- Bool_t kpr0 = kTRUE;
- const AliESDVertex *vertex = ESDEvent->GetPrimaryVertexTracks();
- if(vertex->GetNContributors()<1) {
- // SPD vertex
- vertex = ESDEvent->GetPrimaryVertexSPD();
- if(vertex->GetNContributors()<1) {
- // NO VERTEX, SKIP EVENT
- kpr0 = kFALSE;
- }
- }
- const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ()) < 10.);
- // 10 is the common value, unit: cm
- if (hpriv) hpriv->Fill(kpriv);
- if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
- if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
- if(!kpriv)
- return kFALSE;
-
- if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
- if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
- if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
- return kTRUE;
-}
-
-//------------------------------------------------------------------------------
-Bool_t AliCDMesonUtils::CutEvent(const AliAODEvent *AODEvent, TH1 *hpriv,
- TH1 *hpriVtxX, TH1 *hpriVtxY, TH1 *hpriVtxZ,
- TH1 *hpriVtxPos, TH1 *hpriVtxDist)
-{
- //
- // Cut Event for AOD Events, to be combined with the ESD Track Cut
- //
-
- // TODO: no idea about fast or yet, to be thought of
-
- // Primary vertex
- Bool_t kpr0 = kTRUE;
- const AliAODVertex *vertex = AODEvent->GetPrimaryVertex();
- if(vertex->GetNContributors()<1) {
- // SPD vertex
- vertex = AODEvent->GetPrimaryVertexSPD();
- if(vertex->GetNContributors()<1) {
- // NO GOOD VERTEX, SKIP EVENT
- kpr0 = kFALSE;
- }
- }
- const Bool_t kpriv = kpr0 && (fabs(vertex->GetZ())<10.);
- // 10 is the common value, unit: cm
- hpriv->Fill(kpriv);
- if (hpriVtxDist) hpriVtxDist->Fill(vertex->GetZ());
- if (hpriVtxPos && kpr0) hpriVtxPos->Fill(!kpriv);
- if(!kpriv)
- return kFALSE;
-
- if(hpriVtxX) hpriVtxX->Fill(vertex->GetX());
- if(hpriVtxY) hpriVtxY->Fill(vertex->GetY());
- if(hpriVtxZ) hpriVtxZ->Fill(vertex->GetZ());
- return kTRUE;
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::DoVZEROStudy(const AliESDEvent *ESDEvent,
- TObjArray* hists, Int_t run)
-{
- //
- //
- // IMPORTANT: the order of the histograms here and in
- // AliCDMesonBase::GetHistVZEROStudies(..) has to match
-
- const AliESDVZERO* esdV0 = ESDEvent->GetVZEROData();
- if (!esdV0) {
- Printf("ERROR: esd V0 not available");
- return;
- }
-
- // determine trigger decision
- AliTriggerAnalysis triggerAnalysis;
- //const Bool_t khw = kFALSE;
-
- //Float_t v0a = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
- // AliTriggerAnalysis::kV0BB) ? 1. : 0.;
- //Float_t v0c = (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
- // AliTriggerAnalysis::kV0BB) ? 1. : 0.;
-
- // obtain OCDB objects
- AliCDBManager *man = AliCDBManager::Instance();
- TString cdbpath;
- if (man->IsDefaultStorageSet()) {
- const AliCDBStorage *dsto = man->GetDefaultStorage();
- cdbpath = TString(dsto->GetBaseFolder());
- }
- else { //should not be used!
- man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
- cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
- }
- man->SetSpecificStorage("VZERO/Trigger/Data",cdbpath);
- man->SetSpecificStorage("VZERO/Calib/Data",cdbpath);
- man->SetRun(run);
-
- AliCDBEntry *ent1 = man->Get("VZERO/Trigger/Data");
- if (!ent1) {
- printf("AliCDMesonUtils failed loading VZERO trigger entry from OCDB\n");
- return;
- }
- AliVZEROTriggerData *trigData = (AliVZEROTriggerData*)ent1->GetObject();
- if (!trigData) {
- printf("AliCDMesonUtils failed loading VZERO trigger data from OCDB\n");
- return;
- }
-
- AliCDBEntry *ent2 = man->Get("VZERO/Calib/Data");
- if (!ent2) {
- printf("AliCDMesonUtils failed loading VZERO calib entry from OCDB\n");
- return;
- }
- AliVZEROCalibData *calData = (AliVZEROCalibData*)ent2->GetObject();
- if (!calData) {
- printf("AliCDMesonUtils failed loading VZERO calibration data from OCDB\n");
- return;
- }
-
- // fill histograms
- Int_t pmtHist = 0;
- Int_t iPMT = 0;
- for (iPMT = 0; iPMT < 64; ++iPMT) {
- Int_t board = AliVZEROCalibData::GetBoardNumber(iPMT);
- Int_t channel = AliVZEROCalibData::GetFEEChannelNumber(iPMT);
- ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
- trigData->GetPedestalCut(0, board, channel));
- // calData->GetDiscriThr(iPMT));
- }
- pmtHist = iPMT;
- for (iPMT = 0; iPMT < 64; ++iPMT) {
- ((TH2*)hists->At(iPMT+pmtHist))->Fill(esdV0->GetAdc(iPMT),
- esdV0->GetMultiplicity(iPMT));
- }
- /*
- // not used in 2010 for pp
- ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeA(), v0a);
- ((TH2*)hists->At(pmtHist++))->Fill((Float_t)esdV0->GetTriggerChargeC(), v0c);
- */
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetGapConfig(const AliESDEvent *ESDEvent,
- TH2 *hitMapSPDinner, TH2 *hitMapSPDouter,
- TH2 *hitMapSPDtrklt, TH2 *hitMapFMDa,
- TH2 *hitMapFMDc, TH1 **fmdSums,
- TH2 *TPCGapDCAaSide, TH2 *TPCGapDCAcSide)
-{
- //
- // GetGapConfigAndTracks
- //
- // retrieves the gap configuration of a track and returns it as
- // an bit vector
- // kBaseLine ensures, that this event is valid
- // + is equivalent to | in this case
- Int_t gapConfig = AliCDMesonBase::kBitBaseLine + GetV0(ESDEvent)
- + GetFMD(ESDEvent, hitMapFMDa, hitMapFMDc, fmdSums)
- + GetSPD(ESDEvent, hitMapSPDinner, hitMapSPDouter, hitMapSPDtrklt);
- if (gapConfig == AliCDMesonBase::kBitBaseLine) {
- gapConfig += GetTPC(ESDEvent, TPCGapDCAaSide, TPCGapDCAcSide);
- }
- else {
- gapConfig += GetTPC(ESDEvent, 0x0, 0x0);
- }
- if (GetFastORmultiplicity(ESDEvent) > 0) {
- gapConfig += AliCDMesonBase::kBitCentAct;
- }
- return gapConfig; // + GetZDC(ESDEvent);
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::FillEtaPhiMap(const AliVEvent *event,
- const AliCDMesonTracks* tracks, TH2 *map,
- TH2 *map_c)
-{
- //
- // Fills the eta phi information about all tracks and about the tracks surving
- // the track cuts (CutTrack) into two separate histograms
- //
-
- // all tracks
- for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); ++itrk) {
- if (AliVParticle *trk = event->GetTrack(itrk)) {
- if (map) {
- map->Fill(trk->Eta(), trk->Phi());
- }
- }
- }
-
- // tracks that survived the cuts
- for (Int_t itrk = 0; itrk < tracks->GetCombinedTracks(); ++itrk) {
- if (map_c) {
- map_c->Fill(tracks->GetTrack(itrk)->Eta(), tracks->GetTrack(itrk)->Phi());
- }
- }
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetMultFMD(const AliESDEvent *ESDEvent, Int_t& fmdA,
- Int_t& fmdC, Float_t *fmdSums)
-{
- //
- // Multiplicity seen by FMD
- //
- // WARNING: this function is only working with a modified AliRoot so far
-
-#ifdef STD_ALIROOT
- fmdA = FMDHitCombinations(ESDEvent, 0);
- fmdC = FMDHitCombinations(ESDEvent, 1);
-#else
- AliTriggerAnalysis triggerAnalysis;
- triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
- triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide, &fmdA);
- triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide, &fmdC);
-#endif
-
- if (fmdSums) {
- const AliESDFMD* fmdData =
- (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
- for (UShort_t det=1; det<=3; det++) {
- Int_t nRings = (det == 1 ? 1 : 2);
- for (UShort_t ir = 0; ir < nRings; ir++) {
- Char_t ring = (ir == 0 ? 'I' : 'O');
- UShort_t nsec = (ir == 0 ? 20 : 40);
- UShort_t nstr = (ir == 0 ? 512 : 256);
- for (UShort_t sec =0; sec < nsec; sec++) {
- for (UShort_t strip = 0; strip < nstr; strip++) {
- Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
- if (mult == AliESDFMD::kInvalidMult) continue;
-
- if (det == 1 && ring == 'I') fmdSums[0] += mult;
- else if (det == 2 && ring == 'I') fmdSums[1] += mult;
- else if (det == 2 && ring == 'O') fmdSums[2] += mult;
- else if (det == 3 && ring == 'I') fmdSums[3] += mult;
- else if (det == 3 && ring == 'O') fmdSums[4] += mult;
- }
- }
- }
- }
- }
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetMultSPD(const AliESDEvent * ESDEvent, Int_t& spdIA,
- Int_t& spdIC, Int_t& spdOA, Int_t& spdOC)
-{
- //
- // Retrieves the multiplicity seen by the SPD FastOR
- //
-
- Int_t nfoctr[10];
- GetNFO(ESDEvent, "]0.9[", nfoctr, 0x0, 0x0);
-
- spdIA = nfoctr[kIPA];
- spdIC = nfoctr[kIPC];
- spdOA = nfoctr[kOPA];
- spdOC = nfoctr[kOPC];
-}
-
-
-//==============================================================================
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetV0(const AliESDEvent * ESDEvent)
-{
- //
- //GetV0
- //
-
- AliTriggerAnalysis triggerAnalysis;
- const Bool_t khw = kFALSE;
- const Bool_t v0A =
- (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kASide, khw) ==
- AliTriggerAnalysis::kV0BB);
- const Bool_t v0C =
- (triggerAnalysis.V0Trigger(ESDEvent, AliTriggerAnalysis::kCSide, khw) ==
- AliTriggerAnalysis::kV0BB);
-
- return v0A * AliCDMesonBase::kBitV0A + v0C * AliCDMesonBase::kBitV0C;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetFMD(const AliESDEvent *ESDEvent, TH2 *hitMapFMDa,
- TH2 *hitMapFMDc, TH1 **fmdSums)
-{
- //
- // GetFMD
- //
-
- AliTriggerAnalysis triggerAnalysis;
- triggerAnalysis.SetFMDThreshold(0.3, 0.5); // parameters got from FMD
- const Bool_t fmdA =
- triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kASide);
- const Bool_t fmdC =
- triggerAnalysis.FMDTrigger(ESDEvent, AliTriggerAnalysis::kCSide);
-
- //printf("FR - GetFMD\n");
-
- // prepartions for a charge summation algorithm
- Bool_t hitMaps = (Bool_t)(hitMapFMDa && hitMapFMDc);
- Bool_t calcSum = fmdSums ?
- (fmdSums[0] && fmdSums[1] && fmdSums[2] && fmdSums[3] && fmdSums[4]) :
- kFALSE; // are the histograms defined?
- Float_t sum[] = { 0., 0., 0., 0., 0. };
- // summed multiplicity in the FMD detectors
-
- // hit map generation
- if (hitMaps || calcSum) {
- const AliESDFMD* fmdData =
- (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
-
- for (UShort_t det=1; det<=3;det++) {
- Int_t nRings = (det == 1 ? 1 : 2);
- for (UShort_t ir = 0; ir < nRings; ir++) {
- Char_t ring = (ir == 0 ? 'I' : 'O');
- UShort_t nsec = (ir == 0 ? 20 : 40);
- UShort_t nstr = (ir == 0 ? 512 : 256);
- for (UShort_t sec =0; sec < nsec; sec++) {
- for (UShort_t strip = 0; strip < nstr; strip++) {
- Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
- if (mult == AliESDFMD::kInvalidMult) continue;
-
- if (calcSum) {
- if (det == 1 && ring == 'I') sum[0] += mult;
- else if (det == 2 && ring == 'I') sum[1] += mult;
- else if (det == 2 && ring == 'O') sum[2] += mult;
- else if (det == 3 && ring == 'I') sum[3] += mult;
- else if (det == 3 && ring == 'O') sum[4] += mult;
- }
-
- if (hitMaps) { // care about hit map specific information
- const Float_t eta = fmdData->Eta(det,ring,sec,strip);
- const Float_t phi =
- fmdData->Phi(det,ring,sec,strip) / 180. * TMath::Pi();
- //printf("FR - GetFMD: %f %f %f\n", eta, phi, mult);
- if (eta != AliESDFMD::kInvalidEta) {
- if ((-3.5 < eta) && (eta < -1.5)) {
- hitMapFMDc->Fill(eta, phi, mult); // use mult as weight
- }
- else if ((1.5 < eta) && (eta < 5.5)) {
- hitMapFMDa->Fill(eta, phi, mult); // use mult as weight
- }
- }
- }
- }
- }
- }
- }
- }
-
- if (calcSum) {
- //printf("DEBUG -- SUM(%f,%f,%f,%f,%f)\n", sum[0], sum[1], sum[2], sum[3],
- // sum[4]);
- for (UInt_t i = 0; i < 5; i++) { //
- fmdSums[i]->Fill(sum[i]);
- }
- }
-
- return fmdA * AliCDMesonBase::kBitFMDA + fmdC * AliCDMesonBase::kBitFMDC;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetSPD(const AliESDEvent *ESDEvent, TH2 *hitMapSPDinner,
- TH2 *hitMapSPDouter, TH2 *hitMapSPDtrklt)
-{
- //
- // GetSPD
- //
-
- Int_t nfoctr[10];
- GetNFO(ESDEvent, "]0.9[", nfoctr, hitMapSPDinner, hitMapSPDouter);
- // get multiplicity from fastOR and fill corresponding hit maps
-
- if (hitMapSPDtrklt) FillSPDtrkltMap(ESDEvent, hitMapSPDtrklt);
- // fill tracklet hit map
-
- const Int_t ipA = nfoctr[kIPA]; // inner layer A side
- const Int_t ipC = nfoctr[kIPC]; // inner layer C side
- const Int_t opA = nfoctr[kOPA]; // outer layer A side
- const Int_t opC = nfoctr[kOPC]; // outer layer C side
-
- const Bool_t spdA = ipA + opA; // A side hit?
- const Bool_t spdC = ipC + opC; // C side hit?
-
- return spdA * AliCDMesonBase::kBitSPDA + spdC * AliCDMesonBase::kBitSPDC;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetTPC(const AliESDEvent * ESDEvent, TH2 *TPCGapDCAaSide,
- TH2 *TPCGapDCAcSide)
-{
- //
- //GetTPC
- //
-
- const Double_t etacut = 0.9;
- Int_t nA = 0;
- Int_t nC = 0;
-
- AliESDtrackCuts cuts;
- cuts.SetMaxDCAToVertexXY(0.1);
- cuts.SetMaxDCAToVertexZ(2.);
-
- for(Int_t itrack = 0; itrack < ESDEvent->GetNumberOfTracks(); itrack++){
- const AliESDtrack* esdtrack = ESDEvent->GetTrack(itrack);
- Float_t b[2];
- Float_t bCov[3];
- esdtrack->GetImpactParameters(b,bCov);
-
- if (bCov[0]<=0 || bCov[2]<=0) {
- printf("AliCDMesonUtils - Estimated b resolution lower or equal zero!\n");
- bCov[0]=0;
- bCov[2]=0;
- }
-
- Float_t dcaToVertexXY = b[0];
- Float_t dcaToVertexZ = b[1];
- if (esdtrack->Eta() > etacut) {
- if (cuts.AcceptTrack(esdtrack)) nA++;
- if (TPCGapDCAaSide) {
- TPCGapDCAaSide->Fill(dcaToVertexXY, dcaToVertexZ);
- }
- }
- else if (esdtrack->Eta() < -etacut) {
- if (cuts.AcceptTrack(esdtrack)) nC++;
- if (TPCGapDCAcSide) {
- TPCGapDCAcSide->Fill(dcaToVertexXY, dcaToVertexZ);
- }
- }
- }
-
- const Bool_t tpcA = nA;
- const Bool_t tpcC = nC;
-
- return tpcA * AliCDMesonBase::kBitTPCA + tpcC * AliCDMesonBase::kBitTPCC;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetZDC(const AliESDEvent * ESDEvent)
-{
- //
- //GetZDC
- //
-
- const Int_t qa = ESDEvent->GetESDZDC()->GetESDQuality();
- Bool_t zdcA = kFALSE, zdcC = kFALSE;
- for(Int_t ii=0; ii<6; ii++){
- if(qa & (1<<ii)){
- if(ii<4) zdcA = kTRUE;
- else zdcC = kTRUE;
- }
- }
-
- return zdcA * AliCDMesonBase::kBitZDCA + zdcC * AliCDMesonBase::kBitZDCC;
-}
-
-
-//==============================================================================
-//------------------------------------------------------------------------------
-#ifdef STD_ALIROOT
-Int_t AliCDMesonUtils::FMDHitCombinations(const AliESDEvent* ESDEvent,
- Int_t side)
-{
- //
- // copy of the FMDHitCombinations function originating from AliTriggerAnalysis
- //
- // side == 0 -> A side, side > 0 -> C side
-
- // workaround for AliESDEvent::GetFMDData is not const!
- const AliESDFMD* fmdData =
- (const_cast<AliESDEvent*>(ESDEvent))->GetFMDData();
- if (!fmdData)
- {
- puts("AliESDFMD not available");
- return -1;
- }
-
- Int_t detFrom = (side == 0) ? 1 : 3;
- Int_t detTo = (side == 0) ? 2 : 3;
-
- Int_t triggers = 0;
- Float_t totalMult = 0;
- for (UShort_t det=detFrom;det<=detTo;det++) {
- Int_t nRings = (det == 1 ? 1 : 2);
- for (UShort_t ir = 0; ir < nRings; ir++) {
- Char_t ring = (ir == 0 ? 'I' : 'O');
- UShort_t nsec = (ir == 0 ? 20 : 40);
- UShort_t nstr = (ir == 0 ? 512 : 256);
- for (UShort_t sec =0; sec < nsec; sec++) {
- for (UShort_t strip = 0; strip < nstr; strip++) {
- Float_t mult = fmdData->Multiplicity(det,ring,sec,strip);
- if (mult == AliESDFMD::kInvalidMult) continue;
-
- if (mult > 0.3) // fFMDLowCut
- totalMult = totalMult + mult;
- else {
- if (totalMult > 0.5) // fFMDHitCut
- triggers++;
- }
- }
- }
- }
- }
- return triggers;
-}
-#endif
-
-
-//==============================================================================
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::SPDLoadGeom(Int_t run)
-{
- // method to get the gGeomanager
- // it is called at the CreatedOutputObject stage
- // to comply with the CAF environment
-
- AliCDBManager *man = AliCDBManager::Instance();
-
- TString cdbpath;
- if (man->IsDefaultStorageSet()) {
- const AliCDBStorage *dsto = man->GetDefaultStorage();
- cdbpath = TString(dsto->GetBaseFolder());
- //printf("default was set to: %s\n", cdbpath.Data());
- }
- else { //should not be used!
- // man->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
- // would be needed on grid
- man->SetDefaultStorage(gSystem->Getenv("TRAIN_CDB_PATH"));
- cdbpath = TString(gSystem->Getenv("TRAIN_CDB_PATH"));
- }
-
- man->SetSpecificStorage("ITS/Align/Data",cdbpath);
- man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
- man->SetRun(run);
-
- AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
- if (!obj) {
- printf("AliCDMesonUtils failed loading geometry object\n");
- return;
- }
- AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
- AliGeomManager::ApplyAlignObjsFromCDB("ITS");
-}
-
-//------------------------------------------------------------------------------
-Bool_t AliCDMesonUtils::SPDLoc2Glo(Int_t id, const Double_t *loc,
- Double_t *glo)
-{
- //
- //SPDLoc2Glo, do not touch
- //
-
- static TGeoHMatrix mat;
- Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
- if (vid<0) {
- printf("AliCDMesonUtils Did not find module with such ID %d\n",id);
- return kFALSE;
- }
- AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
- mat.LocalToMaster(loc,glo);
- return kTRUE;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::CheckChipEta(Int_t chipKey, TString scut,
- const Double_t vtxPos[],
- TH2 *hitMapSPDinner, TH2 *hitMapSPDouter)
-{
- //
- //CheckChipEta
- //
-
- // retrieves the position in eta for a given chip and applies the cut
- // results:
- // 0 <= out of range
- // -1 <= negative pseudo-rapidity position, in range (C-Side)
- // 1 <= positive pseudo-rapidity position, in range (A-Side)
- //
- // scut: "[0.9" or "]0.9", only 3 digits for the value!!
-
-
- const Bool_t kincl = (scut[0] == '[');
- const TString cutval = scut(1,3);
- const Double_t etacut = fabs(cutval.Atof());
-
- //no eta cut, save time
- if(kincl && etacut>=2)
- return kTRUE;
-
- Int_t etaside = 1;
- //------------------------------- NOT TO TOUCH ------------------------>>
- UInt_t module=999, offchip=999;
- AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
- UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
- if(hs<2) offchip = 4 - offchip; // inversion in the inner layer...
-
- const Int_t col[]={
- hs<2? 0 : 31,
- hs<2? 31 : 0,
- hs<2? 31 : 0,
- hs<2? 0 : 31};
- const Int_t aa[]={0, 0, 255, 255};
- const AliITSsegmentationSPD seg;
-
- for(Int_t ic=0; ic<4; ic++){
- Float_t localchip[3]={0.,0.,0.};
- seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]);
- // local coordinate of the chip center
- //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
- const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
- Double_t glochip[3]={0.,0.,0.};
- if(!SPDLoc2Glo(module,local,glochip)){
- return kFALSE;
- }
-
- //-------------------------------------------------------------------<<
-
- const TVector3 pos(glochip[0]-vtxPos[0], glochip[1]-vtxPos[1],
- glochip[2]-vtxPos[2]);
- //pos.Print();
-
- if (chipKey < 400) { // inner SPD layer
- if (hitMapSPDinner) {
- Double_t phi = pos.Phi(); // output in the range -Pi +Pi
- if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
- const Double_t eta = pos.Eta();
- hitMapSPDinner->Fill(eta, phi);
- }
- }
- else {
- if (hitMapSPDouter) { // outer SPD layer
- Double_t phi = pos.Phi(); // output in the range -Pi +Pi
- if (phi < 0.) phi += TMath::TwoPi(); // remap to the interval [0, TwoPi)
- const Double_t eta = pos.Eta();
- hitMapSPDouter->Fill(eta, phi);
- }
- }
-
- if( kincl && fabs(pos.Eta()) > etacut)
- return kFALSE;
-
- if(!kincl){
- if(fabs(pos.Eta()) < etacut)
- return kFALSE;
- else if(pos.Eta()<0)
- etaside = -1;
- else
- etaside = 1;
- }
- }
-
- return etaside;
-}
-
-
-//------------------------------------------------------------------------------
-void AliCDMesonUtils::GetNFO(const AliESDEvent *ESDEvent, TString etacut,
- Int_t ctr[], TH2 *hitMapSPDinner,
- TH2 *hitMapSPDouter)
-{
- //
- // GetNFO
- //
- // analyzes the SPD fastOR for a given eta range and returns
- // an array with the number of hits in:
-
- Int_t ninner=0; // inner layer
- Int_t nouter=0; // outer layer
- Int_t ipA = 0; // inner layer A side
- Int_t ipC = 0; // inner layer C side
- Int_t opA = 0; // outer layer A side
- Int_t opC = 0; // outer layer C side
-
- const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
-
- // position of the primary vertex
- Double_t tmp[3] = { 0., 0., 0. };
- ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
- Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
-
-
- for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
- if(mult->TestFastOrFiredChips(iChipKey)){
- // here you check if the FastOr bit is 1 or 0
- const Int_t iseta = CheckChipEta(iChipKey, etacut, vtxPos, hitMapSPDinner,
- hitMapSPDouter);
- if(iseta==0)
- continue;
-
- if(iChipKey<400) {
- ninner++; // here you count the FastOr bits in the inner layer
- if(iseta>0)
- ipA ++;
- else
- ipC ++;
- }
- else {
- nouter++; // here you count the FastOr bits in the outer layer
- if(iseta>0)
- opA ++;
- else
- opC ++;
- }
- }
- }
-
- ctr[kInnerPixel]= ninner;
- ctr[kOuterPixel]= nouter;
- ctr[kIPA]= ipA;
- ctr[kIPC]= ipC;
- ctr[kOPA]= opA;
- ctr[kOPC]= opC;
-
- return;
-}
-
-
-//--------------------------------------------------------------------------
-Int_t AliCDMesonUtils::GetFastORmultiplicity(const AliESDEvent* ESDEvent)
-{
- // determine the number of fired fastOR chips in both layers within
- // -0.9 < eta < 0.9
- //
-
- const AliMultiplicity *mult = ESDEvent->GetMultiplicity();
-
- // position of the primary vertex
- Double_t tmp[3] = { 0., 0., 0. };
- ESDEvent->GetPrimaryVertex()->GetXYZ(tmp);
- Double_t vtxPos[3] = { tmp[0], tmp[1], tmp[2] };
-
- Int_t multiplicity = 0;
-
- for (Int_t iChipKey=0; iChipKey < 1200; iChipKey++) {
- if(mult->TestFastOrFiredChips(iChipKey)){
- // here you check if the FastOr bit is 1 or 0
- const Int_t iseta = CheckChipEta(iChipKey, "[0.9]", vtxPos, 0x0, 0x0);
- if(iseta==0)
- continue;
- else
- ++multiplicity;
- }
- }
-
- return multiplicity;
-}
-
-
-//--------------------------------------------------------------------------
-void AliCDMesonUtils::FillSPDtrkltMap(const AliVEvent* event,
- TH2 *hitMapSPDtrklt)
-{
- //
- // fill eta phi map of SPD tracklets
- //
-
- if (hitMapSPDtrklt) {
- if (TString(event->ClassName()) == "AliESDEvent") {
- const AliMultiplicity *mult = ((AliESDEvent*)event)->GetMultiplicity();
- for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
- iTracklet++) {
- Double_t eta = mult->GetEta(iTracklet);
- Double_t phi = mult->GetPhi(iTracklet);
- hitMapSPDtrklt->Fill(eta, phi);
- }
- }
- else if (TString(event->ClassName()) == "AliAODEvent") {
- const AliAODTracklets* mult = ((AliAODEvent*)event)->GetTracklets();
- for (Int_t iTracklet = 0; iTracklet < mult->GetNumberOfTracklets();
- iTracklet++) {
- Double_t eta = -TMath::Log(TMath::Tan(mult->GetTheta(iTracklet)/2.));
- Double_t phi = mult->GetPhi(iTracklet);
- hitMapSPDtrklt->Fill(eta, phi);
- }
- }
- }
-}
-
-
-//==========================================================================
-Int_t AliCDMesonUtils::GetPID(AliPIDResponse *pid, const AliVTrack *trk,
- Int_t mode /* = 0 */)
-{
- // determines PID for ESDs and AODs
- //
- //
-
- if (!pid) return AliCDMesonBase::kBinPIDUnknown; // no pid available
-
- Double_t tpcpion = -999., tpckaon = -999., tpcproton = -999.,
- tpcelectron = -999.;
- Double_t tofpion = -999., tofkaon = -999., tofproton = -999.,
- tofelectron = -999.;
- Double_t itspion = -999., itskaon = -999., itsproton = -999.,
- itselectron = -999.;
-
-
- // check whether the track was pure ITS standalone track (has not left TPC)
- const Bool_t kits = !(trk->GetStatus() & AliESDtrack::kTPCout);
-
- if (kits) { // do ITS pid
- const Double_t nin = 3.;
-
- itspion = pid->NumberOfSigmasITS( trk, AliPID::kPion );
- itskaon = pid->NumberOfSigmasITS( trk, AliPID::kKaon );
- itsproton = pid->NumberOfSigmasITS( trk, AliPID::kProton );
- itselectron = pid->NumberOfSigmasITS( trk, AliPID::kElectron );
-
- if(fabs(itspion) < nin) // inclusion only
- return AliCDMesonBase::kBinPion;
- else if(fabs(itskaon) < nin)
- return AliCDMesonBase::kBinKaon;
- else if(fabs(itsproton) < nin)
- return AliCDMesonBase::kBinProton;
- else if(fabs(itselectron) < nin)
- return AliCDMesonBase::kBinElectron;
- else{
- return AliCDMesonBase::kBinPIDUnknown;
- }
- }
-
- tpcpion = pid->NumberOfSigmasTPC( trk, AliPID::kPion );
- tpckaon = pid->NumberOfSigmasTPC( trk, AliPID::kKaon );
- tpcproton = pid->NumberOfSigmasTPC( trk, AliPID::kProton );
- tpcelectron = pid->NumberOfSigmasTPC( trk, AliPID::kElectron );
-
- // derive information, whether tof pid is available
- const Bool_t ka = !(trk->GetStatus() & AliESDtrack::kTOFmismatch);
- const Bool_t kb = (trk->GetStatus() & AliESDtrack::kTOFpid);
- const Bool_t ktof = ka && kb;
-
- // retrieve TOF information if available
- if(ktof){
- tofpion = pid->NumberOfSigmasTOF(trk, AliPID::kPion);
- tofkaon = pid->NumberOfSigmasTOF(trk, AliPID::kKaon);
- tofproton = pid->NumberOfSigmasTOF(trk, AliPID::kProton);
- tofelectron = pid->NumberOfSigmasTOF(trk, AliPID::kElectron);
- }
-
- if (mode == 0) {
- // 2010 cuts for PID
- const Double_t nin = 3.;
-
- // inclusion + exclusion (either TPC or TOF)
- if((fabs(tpcpion) < nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
- && fabs(tpcelectron) > nin) ||
- (fabs(tofpion) < nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
- && fabs(tofelectron) > nin))
- return AliCDMesonBase::kBinPionE;
- else if((fabs(tpcpion) > nin && fabs(tpckaon) < nin && fabs(tpcproton) > nin
- && fabs(tpcelectron) > nin) ||
- (fabs(tofpion) > nin && fabs(tofkaon) < nin && fabs(tofproton) > nin
- && fabs(tofelectron) > nin))
- return AliCDMesonBase::kBinKaonE;
- else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) < nin
- && fabs(tpcelectron) > nin) ||
- (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) < nin
- && fabs(tofelectron) > nin))
- return AliCDMesonBase::kBinProtonE;
- else if((fabs(tpcpion) > nin && fabs(tpckaon) > nin && fabs(tpcproton) > nin
- && fabs(tpcelectron) < nin) ||
- (fabs(tofpion) > nin && fabs(tofkaon) > nin && fabs(tofproton) > nin
- && fabs(tofelectron) < nin))
- return AliCDMesonBase::kBinElectronE;
- else if(fabs(tpcpion) < nin && fabs(tofpion) < nin) // inclusion (TPC + TOF)
- return AliCDMesonBase::kBinPion;
- else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
- return AliCDMesonBase::kBinKaon;
- else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
- return AliCDMesonBase::kBinProton;
- else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
- return AliCDMesonBase::kBinElectron;
- else{
- return AliCDMesonBase::kBinPIDUnknown;
- }
- }
- else if (mode == 1) {
- // 2011 cuts for PID in LHC11f
- // TPC: [-3,5] sigma (pion)
- // TOF: 3 sigma for all,
- // only Pion is tuned!
- // ONLY INCLUSION CUTS NO EXCLUSION
- const Double_t nin = 3.;
-
- if(tpcpion < 4. && tpcpion > -2. && fabs(tofpion) < -3. )
- return AliCDMesonBase::kBinPion;
- else if(fabs(tpckaon) < nin && fabs(tofkaon) < nin)
- return AliCDMesonBase::kBinKaon;
- else if(fabs(tpcproton) < nin && fabs(tofproton) < nin)
- return AliCDMesonBase::kBinProton;
- else if(fabs(tpcelectron) < nin && fabs(tofelectron) < nin)
- return AliCDMesonBase::kBinElectron;
- else{
- return AliCDMesonBase::kBinPIDUnknown;
- }
- }
- return AliCDMesonBase::kBinPIDUnknown;
-}
-
-
-//==========================================================================
-Int_t AliCDMesonUtils::GetPID(Int_t pdgCode)
-{
- //
- // determine particle type based on PDG code
- //
-
- if (TMath::Abs(pdgCode) == 211) return AliCDMesonBase::kBinPionE;
- else if (TMath::Abs(pdgCode) == 321) return AliCDMesonBase::kBinKaonE;
- else if (TMath::Abs(pdgCode) == 2212) return AliCDMesonBase::kBinProtonE;
- else if (TMath::Abs(pdgCode) == 11) return AliCDMesonBase::kBinElectronE;
- else return AliCDMesonBase::kBinPIDUnknown;
-}
-
-
-//------------------------------------------------------------------------------
-Int_t AliCDMesonUtils::CombinePID(const Int_t pid[])
-{
- //
- // combine the PID result
- //
-
- // determine return value
- if (pid[0] == pid[1]) { // same result for both tracks
- return pid[0];
- }
- // one track identified with exclusion the other only without
- else if ((pid[0] == AliCDMesonBase::kBinPionE &&
- pid[1] == AliCDMesonBase::kBinPion) ||
- (pid[1] == AliCDMesonBase::kBinPionE &&
- pid[0] == AliCDMesonBase::kBinPion)) {
- return AliCDMesonBase::kBinPion;
- }
- else if ((pid[0] == AliCDMesonBase::kBinKaonE &&
- pid[1] == AliCDMesonBase::kBinKaon) ||
- (pid[1] == AliCDMesonBase::kBinKaonE &&
- pid[0] == AliCDMesonBase::kBinKaon)) {
- return AliCDMesonBase::kBinKaon;
- }
- else if ((pid[0] == AliCDMesonBase::kBinProtonE &&
- pid[1] == AliCDMesonBase::kBinProton) ||
- (pid[1] == AliCDMesonBase::kBinProtonE &&
- pid[0] == AliCDMesonBase::kBinProton)) {
- return AliCDMesonBase::kBinProton;
- }
- else if ((pid[0] == AliCDMesonBase::kBinElectronE &&
- pid[1] == AliCDMesonBase::kBinElectron) ||
- (pid[1] == AliCDMesonBase::kBinElectronE &&
- pid[0] == AliCDMesonBase::kBinElectron)) {
- return AliCDMesonBase::kBinElectron;
- }
- // one track identified and one not
- else if (((pid[0] == AliCDMesonBase::kBinPionE ||
- pid[0] == AliCDMesonBase::kBinPion) &&
- pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
- (pid[1] == AliCDMesonBase::kBinPion &&
- pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
- return AliCDMesonBase::kBinSinglePion;
- }
- else if (((pid[0] == AliCDMesonBase::kBinKaonE ||
- pid[0] == AliCDMesonBase::kBinKaon)&&
- pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
- (pid[1] == AliCDMesonBase::kBinKaon &&
- pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
- return AliCDMesonBase::kBinSingleKaon;
- }
- else if (((pid[0] == AliCDMesonBase::kBinProtonE ||
- pid[0] == AliCDMesonBase::kBinProton) &&
- pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
- (pid[1] == AliCDMesonBase::kBinProton &&
- pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
- return AliCDMesonBase::kBinSingleProton;
- }
- else if (((pid[0] == AliCDMesonBase::kBinElectronE ||
- pid[0] == AliCDMesonBase::kBinElectron) &&
- pid[1] == AliCDMesonBase::kBinPIDUnknown) ||
- (pid[1] == AliCDMesonBase::kBinElectron &&
- pid[0] == AliCDMesonBase::kBinPIDUnknown)) {
- return AliCDMesonBase::kBinSingleElectron;
- }
- else
- return AliCDMesonBase::kBinPIDUnknown;
-}
-
-
-//------------------------------------------------------------------------------
-TLorentzVector AliCDMesonUtils::GetKinematics(const Double_t *pa,
- const Double_t *pb,
- Double_t ma,
- Double_t mb, Double_t& cts)
-{
- //
- //get kinematics, cts = cos(theta_{#star})
- //
-
- TLorentzVector va, vb;
- va.SetXYZM(pa[0], pa[1], pa[2], ma);
- vb.SetXYZM(pb[0], pb[1], pb[2], mb);
- const TLorentzVector sumv = va+vb;
-
- const TVector3 bv = -sumv.BoostVector();
-
- va.Boost(bv);
- vb.Boost(bv);
-
- // 3-vectors in the restframe of the mother particle
- const TVector3 pra = va.Vect();
- const TVector3 prb = vb.Vect();
- const TVector3 diff = pra - prb; // their difference
-
- cts = (diff.Mag2()-pra.Mag2()-prb.Mag2()) / (-2.*pra.Mag()*prb.Mag());
- // cosine theta star, calculated according to the law-of-cosines
-
- return sumv;
-}
-
-
-//------------------------------------------------------------------------------
-Double_t AliCDMesonUtils::GetOA(const Double_t *pa, const Double_t *pb)
-{
- //
- //cosOpeningAngle
- //
-
- TVector3 va, vb;
- va.SetXYZ(pa[0], pa[1], pa[2]);
- vb.SetXYZ(pb[0], pb[1], pb[2]);
-
- const TVector3 ua = va.Unit();
- const TVector3 ub = vb.Unit();
-
- return ua.Dot(ub);
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff