[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / extra / AliAnalysisTaskResonanceQA.cxx

#include <Riostream.h>

#include <TChain.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
#include <TParticle.h>
#include <TLorentzVector.h>

#include "AliESDpid.h"
#include "AliESDtrack.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrackCuts.h"

#include "AliPID.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"

#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliMCParticle.h"
#include "AliGenEventHeader.h"

#include "AliAnalysisTaskResonanceQA.h"

//
// analysis task creating basic QA plots for resonance particles
// Author: Ayben Karasu Uysal
// Computes some QA histograms useful for checking productions and data
// both for counting resonances inside MC, and for checking PID performances
//
// Revised by A. Pulvirenti
//


ClassImp(AliAnalysisTaskResonanceQA)

//________________________________________________________________________
AliAnalysisTaskResonanceQA::AliAnalysisTaskResonanceQA(const char *name) :
   AliAnalysisTaskSE(name), 
   fT0(AliESDpid::kTOF_T0),
   fPrimaryThr(1E-6),
   fVz(10.0),
   fOutputList(0),
   fSelectedEvts(0),
   fEventVz(0),
   fdEdxTPC(0),
   fdEdxITS(0),
   fTOFpid(0),
   fDCAXYvsPtBeforeCuts(0), 
   fDCAZvsPtBeforeCuts(0),
   fNClusterPtBeforeCuts(0), 
   fNFindableClusterPtBeforeCuts(0), 
   fNCrossedRowsTPCPtBeforeCuts(0),
   fProducedParticles(0), 
   fESD(0), 
   fESDpid(0),      
   fTrackCuts(0)
{
//   
// Constructor
//

   // Define input and output slots here
   // Input slot #0 works with a TChain
   DefineInput(0, TChain::Class());
   
   // Output slot #0 id reserved by the base class for AOD
   // Output slot #1 writes into a TH1 container
   DefineOutput(1, TList::Class());
   
   // setup to NULL all remaining histos
   Int_t i;
   for (i = 0; i < kResonances; i++) fRsnYPt[0][i] = fRsnYPt[1][i] = 0;
}

//________________________________________________________________________
const char* AliAnalysisTaskResonanceQA::EvtName(EEvtType type) const
{
//
// Return a short string with name of resonance
//

   switch(type) {
      case kBadNoVertex   : return "NoVertex";
      case kBadVertexOutZ : return "FarVertex";
      case kGoodEvent     : return "Good";
      default             : return "Unknown";
   }
}

//________________________________________________________________________
const char* AliAnalysisTaskResonanceQA::RsnName(ERsn type) const
{
//
// Return a short string with name of resonance
//

   switch(type) {
      case kPhi       : return "Phi";
      case kKStar0    : return "KStar0";
      case kRho       : return "Rho";
      case kLambdaStar: return "LambdaStar";
      case kXiStar0   : return "XiStar0";
      case kSigmaStarP: return "SigmaStarP";
      case kSigmaStarM: return "SigmaStarM";
      case kDeltaPP   : return "DeltaPP";
      default         : return "Unknown";
   }
}

//________________________________________________________________________
const char* AliAnalysisTaskResonanceQA::RsnSymbol(ERsn type) const
{
//
// Return a short string with name of resonance
//

   switch(type) {
      case kPhi       : return "#phi";
      case kKStar0    : return "K^{*0}";
      case kRho       : return "#rho";
      case kLambdaStar: return "#Lambda^{*}";
      case kXiStar0   : return "#Xi^{*0}";
      case kSigmaStarP: return "#Sigma^{*+}";
      case kSigmaStarM: return "#Sigma^{*-}";
      case kDeltaPP   : return "#Delta^{++}";
      default         : return "Unknown";
   }
}

//________________________________________________________________________
Int_t AliAnalysisTaskResonanceQA::RsnPDG(ERsn type) const
{
//
// Return PDG code of resonance
//

   switch(type) {
      case kPhi       : return  333;
      case kKStar0    : return  313;
      case kRho       : return  113;
      case kLambdaStar: return 3124;
      case kXiStar0   : return 3324;
      case kSigmaStarP: return 3224;
      case kSigmaStarM: return 3114;
      case kDeltaPP   : return 2224;
      default         : return    0;
   }
}

//________________________________________________________________________
void AliAnalysisTaskResonanceQA::UserCreateOutputObjects()
{
//
// Create histograms (called once)
//

   Int_t i;
   Int_t vESDdEdxBin = 1000;  Double_t vESDdEdxMin = 0;  Double_t vESDdEdxMax = 1000;
   Int_t vESDQPBin = 900;     Double_t vESDQPMin = -4.5; Double_t vESDQPMax = 4.5;
   Int_t vPtBin = 500;        Double_t vPtMin = 0.;      Double_t vPtMax = 10.;
   Int_t vYBin = 600;         Double_t vYMin = -15.;     Double_t vYMax = 15.;

   // standard objects
   fESDpid = new AliESDpid();
   fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
   fTrackCuts->SetPtRange(0.1, 1E20);
   fTrackCuts->SetEtaRange(-0.9, 0.9);

   // create output list
   fOutputList = new TList();
   fOutputList->SetOwner();

   // output histograms for PID signal
   fSelectedEvts = new TH1I("fSelectedEvts", "fSelectedEvts;fSelectedEvts", 4, 0, 4);
   fEventVz = new TH1F("fEventVz", "Vz distribution", 800, -40.0, 40.0);
   fdEdxTPC = new TH2F("fdEdxTPC", "dEdxTPC, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
   fdEdxITS = new TH2F("fdEdxITS", "dEdxITS, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
   fTOFpid = new TH2F("fTOFpid", "TOF PID;p (GeV/c);#beta;", vESDQPBin, vESDQPMin, vESDQPMax, 300, 0.1, 1.2);
   
   // output histograms for track quality
   fDCAXYvsPtBeforeCuts = new TH2F("DCAXYvsPtBeforeCuts", "DCAXYvsPtBeforeCuts;p_{t} (GeV/c);DCAXY", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
   fDCAZvsPtBeforeCuts = new TH2F("DCAZvsPtBeforeCuts", "DCAZvsPtBeforeCuts;p_{t} (GeV/c);DCAZ", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
   fNClusterPtBeforeCuts = new TH2F("NClusterPtBeforeCuts", "NClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
   fNFindableClusterPtBeforeCuts = new TH2F("NFindableClusterPtBeforeCuts", "NFindableClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
   fNCrossedRowsTPCPtBeforeCuts = new TH2F("NCrossedRowsTPCPtBeforeCuts", "NCrossedRowsTPCPtBeforeCuts;p_{t} (GeV/c);NCrossedRowsTPC", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
   
   // output histograms for resonances
   fProducedParticles = new TH1I("ProducedParticles", "ProducedParticles;Produced Particles", kResonances, 0, kResonances);
   for (i = 0; i < kResonances; i++) {
      fRsnYPt[0][i] = new TH3F(Form("%s_all" , RsnName(i)), Form("%s -- ALL;p_{t} (GeV/c);Rapidity"    , RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
      fRsnYPt[1][i] = new TH3F(Form("%s_prim", RsnName(i)), Form("%s -- PRIMARY;p_{t} (GeV/c);Rapidity", RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
      fProducedParticles->GetXaxis()->SetBinLabel(i + 1, RsnSymbol(i));
   }
   
   // add histogram to list
   fOutputList->Add(fSelectedEvts);
   fOutputList->Add(fEventVz);
   fOutputList->Add(fdEdxTPC);
   fOutputList->Add(fdEdxITS);
   fOutputList->Add(fTOFpid);
   fOutputList->Add(fDCAXYvsPtBeforeCuts);
   fOutputList->Add(fDCAZvsPtBeforeCuts);
   fOutputList->Add(fNClusterPtBeforeCuts);
   fOutputList->Add(fNFindableClusterPtBeforeCuts);
   fOutputList->Add(fNCrossedRowsTPCPtBeforeCuts);
   fOutputList->Add(fProducedParticles);
   for (i = 0; i < kResonances; i++) {
      fOutputList->Add(fRsnYPt[0][i]);
      fOutputList->Add(fRsnYPt[1][i]);
   }
   
   PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskResonanceQA::UserExec(Option_t *)
{
//
// Execute computations
//

   // first bin of this histogram counts processed events
   // second bin counts events passing physics selection
   // all events not passing physics selections are skipped
   fSelectedEvts->Fill(0.1);
   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
   if (!isSelected) return;
   
   // count CINT1B selected events
   fSelectedEvts->Fill(1.1);
   
   // try to cast input event to ESD and loop on tracks
   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
   
   // check if event is accepted w.r. to all additional selection criteria
   Bool_t hasVertex = kFALSE, goodVertex = kFALSE;
   Double_t vz = 1E20, vzsign = 1E20;
   if (fESD) {
      // check primary vertex:
      // we accept only tracks/SPD with at least 1 contributor
      const AliESDVertex *v0 = fESD->GetPrimaryVertexTracks();
      if (v0) {
         if (v0->GetNContributors() > 0) {
            vzsign = v0->GetZv();
            vz = TMath::Abs(vzsign);
            hasVertex = kTRUE;
         } else {
            v0 = fESD->GetPrimaryVertexSPD();
            if (v0) {
               if (v0->GetNContributors() > 0) {
                  vzsign = v0->GetZv();
                  vz = TMath::Abs(vzsign);
                  hasVertex = kTRUE;
               }
            }
         }
      }
      if (hasVertex) {
         goodVertex = (vz <= fVz);
         fSelectedEvts->Fill(2.1);
         fEventVz->Fill(vzsign);
      }
   }
   
   // determine event quality
   Bool_t eventAccepted = kFALSE;
   Int_t evtQuality = kEvtTypes;
   if (hasVertex) {
      if (goodVertex) {
         evtQuality = kGoodEvent;
         eventAccepted = kTRUE;
      } else {
         evtQuality = kBadVertexOutZ;
      }
   } else {
      evtQuality = kBadNoVertex;
   }
   
   // Message
   static Int_t evNum = -1;
   evNum++;
   switch (evtQuality) {
      case kBadNoVertex  : AliInfo(Form("Rejecting event #%5d because it has not a good vertex", evNum)); break;
      case kBadVertexOutZ: AliInfo(Form("Rejecting event #%5d because it is out of range: vz = %5.3f vs. %5.3f", evNum, vz, fVz)); break;
      case kGoodEvent    : AliDebugClass(1, Form("Accepting event #%5d", evNum)); break;
      default            : AliError("This should never appear!!!");
   }
   
   // if event is OK, loop on tracks to fill QA
   if (fESD && eventAccepted) {
      // if we arrive here, the event was processed successfully
      // then we fill last bin in first histogram
      fSelectedEvts->Fill(3.1);
      
      // settings for TOF time zero
      if (fESD->GetTOFHeader()) 
         fESDpid->SetTOFResponse(fESD, AliESDpid::kTOF_T0);
      else {
         Float_t t0spread[10];
         Float_t intrinsicTOFres = 100;
         for (Int_t i = 0; i < 10; i++) t0spread[i] = (TMath::Sqrt(fESD->GetSigma2DiamondZ())) / 0.03; //0.03 to convert from cm to ps
         fESDpid->GetTOFResponse().SetT0resolution(t0spread);
         fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
         fESDpid->SetTOFResponse(fESD, fT0);
      }
      fESDpid->MakePID(fESD, kFALSE, 0.);

      // use some variables to store useful info
      /*MISC   */ Int_t    iTrack, nTracks = fESD->GetNumberOfTracks();
      /*MISC   */ Double_t pt;
      /*QUALITY*/ Float_t  impactXY, impactZ, nCrossedRowsTPC;
      /*ITS    */ Double_t itsSignal = 0;
      /*TPC    */ Double_t ptotSE = 0, tpcSignal, signSE;
      /*TOF    */ Double_t momentum, length, time, beta; //, times[AliPID::kSPECIES], tzeroTrack = 0;

      // loop on tracks
      for (iTrack = 0; iTrack < nTracks; iTrack++) {
         
         // get track
         AliESDtrack* track = fESD->GetTrack(iTrack);
         if (!track) continue;
         
         // get transverse momentum (used for quality histograms)
         pt = track->Pt();
         
         // get charge and reject neutrals
         signSE = track->GetSign();
         if (signSE == 0) continue;

         // get DCA and TPC-related quality params
         track->GetImpactParameters(impactXY, impactZ);
         nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);

         // fill quality histograms when status flags are OK for TPC in+refit and ITS refit
         if (track->IsOn(AliESDtrack::kTPCin) && track->IsOn(AliESDtrack::kTPCrefit) && track->IsOn(AliESDtrack::kITSrefit)) {
            fNClusterPtBeforeCuts        ->Fill(pt, track->GetTPCNcls());
            fNFindableClusterPtBeforeCuts->Fill(pt, track->GetTPCNclsF());
            fNCrossedRowsTPCPtBeforeCuts ->Fill(pt, nCrossedRowsTPC);
            fDCAXYvsPtBeforeCuts         ->Fill(pt, impactXY);
            fDCAZvsPtBeforeCuts          ->Fill(pt, impactZ);
         }
         
         // from now on, work only with tracks passing standard cuts for 2010
         if (!fTrackCuts->AcceptTrack(track)) continue;

         // get momentum (used for PID histograms)
         momentum = track->P();

         // ITS
         itsSignal = track->GetITSsignal();
         fdEdxITS->Fill(signSE * momentum, itsSignal);
         
         // TPC
         if (track->GetInnerParam()) {
            AliExternalTrackParam trackIn1(*track->GetInnerParam());
            ptotSE = trackIn1.GetP();
            tpcSignal = track->GetTPCsignal();
            fdEdxTPC->Fill(signSE * ptotSE, tpcSignal);
         } // end TPC
         
         // TOF (requires checking some more status flags
         if (track->IsOn(AliESDtrack::kTOFpid) && track->IsOn(AliESDtrack::kTOFout) && track->IsOn(AliESDtrack::kTIME) && !track->IsOn(AliESDtrack::kTOFmismatch)) {
            length = track->GetIntegratedLength();
            
            // reject suspiciously small lengths/times
            if (track->GetIntegratedLength() < 350.) continue;
            if (track->GetTOFsignal() < 1E-6) continue;

            // thhese are maybe not needed
            // track->GetIntegratedTimes(times);
            // tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
            
            // get TOF measured time and compute beta
            time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(momentum);
            beta = length / (2.99792457999999984e-02 * time);
            fTOFpid->Fill(signSE * momentum, beta);
         } // end TOF
      } // end track loop
   } // end ESD check
   
   // loop on MC, if MC event and stack are available
   AliMCEvent *mcEvent = MCEvent();
   if (mcEvent) {
      // MC primary vertex
      TArrayF mcVertex(3);
      AliGenEventHeader *genEH = mcEvent->GenEventHeader();
      genEH->PrimaryVertex(mcVertex);
      // loop on particles
      AliStack *stack = mcEvent->Stack();
      if (stack) {
         Int_t iTrack, irsn, pdg, mother, motherPDG, nTracks = stack->GetNtrack();
         Double_t dx, dy, dz, dist;
         TLorentzVector vprod;
         for (iTrack = 0; iTrack < nTracks; iTrack++) {
            
            // get particle
            TParticle *part = stack->Particle(iTrack);
            if (!part) {
               AliError(Form("Could not receive track %d", iTrack));
               continue;
            }
            
            // get PDG code of particle and check physical primary
            pdg = part->GetPdgCode();
            mother = part->GetFirstMother();
            motherPDG = 0;
                        
            // loop over possible resonances
            for (irsn = 0; irsn < kResonances; irsn++) {
               if (TMath::Abs(pdg) == RsnPDG(irsn)) {
                  // debug check
                  if (mother >= 0) {
                     TParticle *p = stack->Particle(mother);
                     motherPDG = p->GetPdgCode();
                  }
                  part->ProductionVertex(vprod);
                  dx = mcVertex[0] - vprod.X();
                  dy = mcVertex[1] - vprod.Y();
                  dz = mcVertex[2] - vprod.Z();
                  dist = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
                  AliDebugClass(1, Form("PDG = %d -- mother ID = %d, pdg = %d -- dist. from MC vertex = %f -- prod. time = %f", pdg, mother, motherPDG, dist, vprod.T())); 
                  // global counter is always filled
                  fRsnYPt[0][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
                  // counter for primaries is filled only if mother is less than 0 (primary)
                  if (dist < fPrimaryThr) {
                     fRsnYPt[1][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
                  }
                  // fill the global histogram
                  fProducedParticles->Fill(irsn + 1);
                  // if one PDG match was found, no need to continue loop
                  break;
               }
            }
         }
      }
   }
   
   PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskResonanceQA::Terminate(Option_t *)
{
//
// Termination
// Quickly check that the output list is there
//
  
   fOutputList = dynamic_cast<TList*>(GetOutputData(1));
   if (!fOutputList) {
      AliError("Output list not found!!!");
      return;
   }
}
Commit	Line	Data
c865cb1d	1	#include <Riostream.h>
	2
	3	#include <TChain.h>
	4	#include <TH1.h>
	5	#include <TH2.h>
f4304ff3	6	#include <TH3.h>
c865cb1d	7	#include <TDatabasePDG.h>
	8	#include <TParticlePDG.h>
	9	#include <TParticle.h>
	10	#include <TLorentzVector.h>
	11
	12	#include "AliESDpid.h"
	13	#include "AliESDtrack.h"
	14	#include "AliESDEvent.h"
	15	#include "AliESDVertex.h"
	16	#include "AliESDtrackCuts.h"
	17
	18	#include "AliPID.h"
	19	#include "AliAnalysisTask.h"
	20	#include "AliAnalysisManager.h"
	21	#include "AliInputEventHandler.h"
	22
	23	#include "AliMCEvent.h"
	24	#include "AliStack.h"
	25	#include "AliMCParticle.h"
	26	#include "AliGenEventHeader.h"
	27
	28	#include "AliAnalysisTaskResonanceQA.h"
	29
035b8eff	30	//
c865cb1d	31	// analysis task creating basic QA plots for resonance particles
c865cb1d	32	// Author: Ayben Karasu Uysal
035b8eff	33	// Computes some QA histograms useful for checking productions and data
	34	// both for counting resonances inside MC, and for checking PID performances
	35	//
	36	// Revised by A. Pulvirenti
	37	//
c865cb1d	38
	39
	40	ClassImp(AliAnalysisTaskResonanceQA)
	41
	42	//________________________________________________________________________
	43	AliAnalysisTaskResonanceQA::AliAnalysisTaskResonanceQA(const char *name) :
	44	AliAnalysisTaskSE(name),
	45	fT0(AliESDpid::kTOF_T0),
	46	fPrimaryThr(1E-6),
e187bd70	47	fVz(10.0),
c865cb1d	48	fOutputList(0),
c865cb1d	49	fSelectedEvts(0),
f4304ff3	50	fEventVz(0),
c865cb1d	51	fdEdxTPC(0),
	52	fdEdxITS(0),
	53	fTOFpid(0),
	54	fDCAXYvsPtBeforeCuts(0),
	55	fDCAZvsPtBeforeCuts(0),
	56	fNClusterPtBeforeCuts(0),
	57	fNFindableClusterPtBeforeCuts(0),
	58	fNCrossedRowsTPCPtBeforeCuts(0),
	59	fProducedParticles(0),
	60	fESD(0),
	61	fESDpid(0),
	62	fTrackCuts(0)
	63	{
	64	//
	65	// Constructor
	66	//
	67
	68	// Define input and output slots here
	69	// Input slot #0 works with a TChain
	70	DefineInput(0, TChain::Class());
	71
	72	// Output slot #0 id reserved by the base class for AOD
	73	// Output slot #1 writes into a TH1 container
	74	DefineOutput(1, TList::Class());
	75
	76	// setup to NULL all remaining histos
	77	Int_t i;
	78	for (i = 0; i < kResonances; i++) fRsnYPt[0][i] = fRsnYPt[1][i] = 0;
f4304ff3	79	}
	80
	81	//________________________________________________________________________
	82	const char* AliAnalysisTaskResonanceQA::EvtName(EEvtType type) const
	83	{
	84	//
	85	// Return a short string with name of resonance
	86	//
	87
	88	switch(type) {
	89	case kBadNoVertex : return "NoVertex";
	90	case kBadVertexOutZ : return "FarVertex";
	91	case kGoodEvent : return "Good";
	92	default : return "Unknown";
	93	}
c865cb1d	94	}
	95
	96	//________________________________________________________________________
035b8eff	97	const char* AliAnalysisTaskResonanceQA::RsnName(ERsn type) const
c865cb1d	98	{
	99	//
	100	// Return a short string with name of resonance
	101	//
	102
	103	switch(type) {
	104	case kPhi : return "Phi";
	105	case kKStar0 : return "KStar0";
	106	case kRho : return "Rho";
	107	case kLambdaStar: return "LambdaStar";
	108	case kXiStar0 : return "XiStar0";
c865cb1d	109	case kSigmaStarP: return "SigmaStarP";
c865cb1d	110	case kSigmaStarM: return "SigmaStarM";
	111	case kDeltaPP : return "DeltaPP";
	112	default : return "Unknown";
	113	}
	114	}
	115
	116	//________________________________________________________________________
035b8eff	117	const char* AliAnalysisTaskResonanceQA::RsnSymbol(ERsn type) const
c865cb1d	118	{
	119	//
	120	// Return a short string with name of resonance
	121	//
	122
	123	switch(type) {
	124	case kPhi : return "#phi";
	125	case kKStar0 : return "K^{*0}";
	126	case kRho : return "#rho";
	127	case kLambdaStar: return "#Lambda^{*}";
	128	case kXiStar0 : return "#Xi^{*0}";
c865cb1d	129	case kSigmaStarP: return "#Sigma^{*+}";
c865cb1d	130	case kSigmaStarM: return "#Sigma^{*-}";
	131	case kDeltaPP : return "#Delta^{++}";
	132	default : return "Unknown";
	133	}
	134	}
	135
	136	//________________________________________________________________________
035b8eff	137	Int_t AliAnalysisTaskResonanceQA::RsnPDG(ERsn type) const
c865cb1d	138	{
	139	//
	140	// Return PDG code of resonance
	141	//
	142
	143	switch(type) {
	144	case kPhi : return 333;
	145	case kKStar0 : return 313;
	146	case kRho : return 113;
	147	case kLambdaStar: return 3124;
	148	case kXiStar0 : return 3324;
c865cb1d	149	case kSigmaStarP: return 3224;
c865cb1d	150	case kSigmaStarM: return 3114;
	151	case kDeltaPP : return 2224;
	152	default : return 0;
	153	}
	154	}
	155
	156	//________________________________________________________________________
	157	void AliAnalysisTaskResonanceQA::UserCreateOutputObjects()
	158	{
	159	//
	160	// Create histograms (called once)
	161	//
	162
	163	Int_t i;
035b8eff	164	Int_t vESDdEdxBin = 1000; Double_t vESDdEdxMin = 0; Double_t vESDdEdxMax = 1000;
	165	Int_t vESDQPBin = 900; Double_t vESDQPMin = -4.5; Double_t vESDQPMax = 4.5;
	166	Int_t vPtBin = 500; Double_t vPtMin = 0.; Double_t vPtMax = 10.;
	167	Int_t vYBin = 600; Double_t vYMin = -15.; Double_t vYMax = 15.;
c865cb1d	168
	169	// standard objects
	170	fESDpid = new AliESDpid();
	171	fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
	172	fTrackCuts->SetPtRange(0.1, 1E20);
	173	fTrackCuts->SetEtaRange(-0.9, 0.9);
	174
	175	// create output list
	176	fOutputList = new TList();
	177	fOutputList->SetOwner();
	178
	179	// output histograms for PID signal
f4304ff3	180	fSelectedEvts = new TH1I("fSelectedEvts", "fSelectedEvts;fSelectedEvts", 4, 0, 4);
f4304ff3	181	fEventVz = new TH1F("fEventVz", "Vz distribution", 800, -40.0, 40.0);
035b8eff	182	fdEdxTPC = new TH2F("fdEdxTPC", "dEdxTPC, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
	183	fdEdxITS = new TH2F("fdEdxITS", "dEdxITS, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
	184	fTOFpid = new TH2F("fTOFpid", "TOF PID;p (GeV/c);#beta;", vESDQPBin, vESDQPMin, vESDQPMax, 300, 0.1, 1.2);
c865cb1d	185
c865cb1d	186	// output histograms for track quality
035b8eff	187	fDCAXYvsPtBeforeCuts = new TH2F("DCAXYvsPtBeforeCuts", "DCAXYvsPtBeforeCuts;p_{t} (GeV/c);DCAXY", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
	188	fDCAZvsPtBeforeCuts = new TH2F("DCAZvsPtBeforeCuts", "DCAZvsPtBeforeCuts;p_{t} (GeV/c);DCAZ", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
	189	fNClusterPtBeforeCuts = new TH2F("NClusterPtBeforeCuts", "NClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
	190	fNFindableClusterPtBeforeCuts = new TH2F("NFindableClusterPtBeforeCuts", "NFindableClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
	191	fNCrossedRowsTPCPtBeforeCuts = new TH2F("NCrossedRowsTPCPtBeforeCuts", "NCrossedRowsTPCPtBeforeCuts;p_{t} (GeV/c);NCrossedRowsTPC", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
c865cb1d	192
	193	// output histograms for resonances
	194	fProducedParticles = new TH1I("ProducedParticles", "ProducedParticles;Produced Particles", kResonances, 0, kResonances);
	195	for (i = 0; i < kResonances; i++) {
f4304ff3	196	fRsnYPt[0][i] = new TH3F(Form("%s_all" , RsnName(i)), Form("%s -- ALL;p_{t} (GeV/c);Rapidity" , RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
f4304ff3	197	fRsnYPt[1][i] = new TH3F(Form("%s_prim", RsnName(i)), Form("%s -- PRIMARY;p_{t} (GeV/c);Rapidity", RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
c865cb1d	198	fProducedParticles->GetXaxis()->SetBinLabel(i + 1, RsnSymbol(i));
	199	}
	200
	201	// add histogram to list
	202	fOutputList->Add(fSelectedEvts);
f4304ff3	203	fOutputList->Add(fEventVz);
c865cb1d	204	fOutputList->Add(fdEdxTPC);
	205	fOutputList->Add(fdEdxITS);
	206	fOutputList->Add(fTOFpid);
	207	fOutputList->Add(fDCAXYvsPtBeforeCuts);
	208	fOutputList->Add(fDCAZvsPtBeforeCuts);
	209	fOutputList->Add(fNClusterPtBeforeCuts);
	210	fOutputList->Add(fNFindableClusterPtBeforeCuts);
	211	fOutputList->Add(fNCrossedRowsTPCPtBeforeCuts);
	212	fOutputList->Add(fProducedParticles);
	213	for (i = 0; i < kResonances; i++) {
	214	fOutputList->Add(fRsnYPt[0][i]);
	215	fOutputList->Add(fRsnYPt[1][i]);
	216	}
	217
	218	PostData(1, fOutputList);
	219	}
	220
	221	//________________________________________________________________________
	222	void AliAnalysisTaskResonanceQA::UserExec(Option_t *)
	223	{
	224	//
	225	// Execute computations
	226	//
	227
	228	// first bin of this histogram counts processed events
	229	// second bin counts events passing physics selection
	230	// all events not passing physics selections are skipped
	231	fSelectedEvts->Fill(0.1);
	232	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
	233	if (!isSelected) return;
	234
035b8eff	235	// count CINT1B selected events
c865cb1d	236	fSelectedEvts->Fill(1.1);
c865cb1d	237
c865cb1d	238	// try to cast input event to ESD and loop on tracks
c865cb1d	239	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
035b8eff	240
035b8eff	241	// check if event is accepted w.r. to all additional selection criteria
f4304ff3	242	Bool_t hasVertex = kFALSE, goodVertex = kFALSE;
f4304ff3	243	Double_t vz = 1E20, vzsign = 1E20;
c865cb1d	244	if (fESD) {
	245	// check primary vertex:
	246	// we accept only tracks/SPD with at least 1 contributor
	247	const AliESDVertex *v0 = fESD->GetPrimaryVertexTracks();
7248f715	248	if (v0) {
f4304ff3	249	if (v0->GetNContributors() > 0) {
	250	vzsign = v0->GetZv();
	251	vz = TMath::Abs(vzsign);
	252	hasVertex = kTRUE;
	253	} else {
7248f715	254	v0 = fESD->GetPrimaryVertexSPD();
7248f715	255	if (v0) {
f4304ff3	256	if (v0->GetNContributors() > 0) {
	257	vzsign = v0->GetZv();
	258	vz = TMath::Abs(vzsign);
	259	hasVertex = kTRUE;
	260	}
7248f715	261	}
7248f715	262	}
c865cb1d	263	}
f4304ff3	264	if (hasVertex) {
	265	goodVertex = (vz <= fVz);
	266	fSelectedEvts->Fill(2.1);
	267	fEventVz->Fill(vzsign);
	268	}
	269	}
	270
	271	// determine event quality
	272	Bool_t eventAccepted = kFALSE;
	273	Int_t evtQuality = kEvtTypes;
	274	if (hasVertex) {
	275	if (goodVertex) {
	276	evtQuality = kGoodEvent;
	277	eventAccepted = kTRUE;
	278	} else {
	279	evtQuality = kBadVertexOutZ;
	280	}
	281	} else {
	282	evtQuality = kBadNoVertex;
	283	}
	284
	285	// Message
	286	static Int_t evNum = -1;
	287	evNum++;
	288	switch (evtQuality) {
	289	case kBadNoVertex : AliInfo(Form("Rejecting event #%5d because it has not a good vertex", evNum)); break;
	290	case kBadVertexOutZ: AliInfo(Form("Rejecting event #%5d because it is out of range: vz = %5.3f vs. %5.3f", evNum, vz, fVz)); break;
	291	case kGoodEvent : AliDebugClass(1, Form("Accepting event #%5d", evNum)); break;
	292	default : AliError("This should never appear!!!");
035b8eff	293	}
c865cb1d	294
035b8eff	295	// if event is OK, loop on tracks to fill QA
	296	if (fESD && eventAccepted) {
	297	// if we arrive here, the event was processed successfully
	298	// then we fill last bin in first histogram
f4304ff3	299	fSelectedEvts->Fill(3.1);
035b8eff	300
c865cb1d	301	// settings for TOF time zero
	302	if (fESD->GetTOFHeader())
	303	fESDpid->SetTOFResponse(fESD, AliESDpid::kTOF_T0);
	304	else {
	305	Float_t t0spread[10];
	306	Float_t intrinsicTOFres = 100;
	307	for (Int_t i = 0; i < 10; i++) t0spread[i] = (TMath::Sqrt(fESD->GetSigma2DiamondZ())) / 0.03; //0.03 to convert from cm to ps
	308	fESDpid->GetTOFResponse().SetT0resolution(t0spread);
	309	fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
	310	fESDpid->SetTOFResponse(fESD, fT0);
	311	}
	312	fESDpid->MakePID(fESD, kFALSE, 0.);
	313
	314	// use some variables to store useful info
	315	/MISC / Int_t iTrack, nTracks = fESD->GetNumberOfTracks();
	316	/MISC / Double_t pt;
	317	/QUALITY/ Float_t impactXY, impactZ, nCrossedRowsTPC;
	318	/ITS / Double_t itsSignal = 0;
	319	/TPC / Double_t ptotSE = 0, tpcSignal, signSE;
	320	/TOF / Double_t momentum, length, time, beta; //, times[AliPID::kSPECIES], tzeroTrack = 0;
	321
	322	// loop on tracks
	323	for (iTrack = 0; iTrack < nTracks; iTrack++) {
	324
	325	// get track
	326	AliESDtrack* track = fESD->GetTrack(iTrack);
	327	if (!track) continue;
	328
	329	// get transverse momentum (used for quality histograms)
	330	pt = track->Pt();
	331
	332	// get charge and reject neutrals
	333	signSE = track->GetSign();
	334	if (signSE == 0) continue;
	335
	336	// get DCA and TPC-related quality params
	337	track->GetImpactParameters(impactXY, impactZ);
	338	nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
	339
	340	// fill quality histograms when status flags are OK for TPC in+refit and ITS refit
	341	if (track->IsOn(AliESDtrack::kTPCin) && track->IsOn(AliESDtrack::kTPCrefit) && track->IsOn(AliESDtrack::kITSrefit)) {
	342	fNClusterPtBeforeCuts ->Fill(pt, track->GetTPCNcls());
	343	fNFindableClusterPtBeforeCuts->Fill(pt, track->GetTPCNclsF());
	344	fNCrossedRowsTPCPtBeforeCuts ->Fill(pt, nCrossedRowsTPC);
	345	fDCAXYvsPtBeforeCuts ->Fill(pt, impactXY);
	346	fDCAZvsPtBeforeCuts ->Fill(pt, impactZ);
	347	}
	348
	349	// from now on, work only with tracks passing standard cuts for 2010
	350	if (!fTrackCuts->AcceptTrack(track)) continue;
	351
	352	// get momentum (used for PID histograms)
	353	momentum = track->P();
	354
	355	// ITS
	356	itsSignal = track->GetITSsignal();
	357	fdEdxITS->Fill(signSE * momentum, itsSignal);
	358
	359	// TPC
	360	if (track->GetInnerParam()) {
	361	AliExternalTrackParam trackIn1(*track->GetInnerParam());
	362	ptotSE = trackIn1.GetP();
	363	tpcSignal = track->GetTPCsignal();
	364	fdEdxTPC->Fill(signSE * ptotSE, tpcSignal);
365	} // end TPC
366
367	// TOF (requires checking some more status flags
368	if (track->IsOn(AliESDtrack::kTOFpid) && track->IsOn(AliESDtrack::kTOFout) && track->IsOn(AliESDtrack::kTIME) && !track->IsOn(AliESDtrack::kTOFmismatch)) {
369	length = track->GetIntegratedLength();
370
371	// reject suspiciously small lengths/times
372	if (track->GetIntegratedLength() < 350.) continue;
373	if (track->GetTOFsignal() < 1E-6) continue;
374
375	// thhese are maybe not needed
376	// track->GetIntegratedTimes(times);
377	// tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
378
379	// get TOF measured time and compute beta
380	time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(momentum);
381	beta = length / (2.99792457999999984e-02 * time);
382	fTOFpid->Fill(signSE * momentum, beta);
383	} // end TOF
384	} // end track loop
385	} // end ESD check
035b8eff	386
	387	// loop on MC, if MC event and stack are available
	388	AliMCEvent *mcEvent = MCEvent();
	389	if (mcEvent) {
	390	// MC primary vertex
	391	TArrayF mcVertex(3);
	392	AliGenEventHeader *genEH = mcEvent->GenEventHeader();
	393	genEH->PrimaryVertex(mcVertex);
	394	// loop on particles
	395	AliStack *stack = mcEvent->Stack();
	396	if (stack) {
	397	Int_t iTrack, irsn, pdg, mother, motherPDG, nTracks = stack->GetNtrack();
	398	Double_t dx, dy, dz, dist;
	399	TLorentzVector vprod;
	400	for (iTrack = 0; iTrack < nTracks; iTrack++) {
	401
	402	// get particle
	403	TParticle *part = stack->Particle(iTrack);
	404	if (!part) {
	405	AliError(Form("Could not receive track %d", iTrack));
	406	continue;
	407	}
	408
	409	// get PDG code of particle and check physical primary
	410	pdg = part->GetPdgCode();
	411	mother = part->GetFirstMother();
	412	motherPDG = 0;
	413
	414	// loop over possible resonances
	415	for (irsn = 0; irsn < kResonances; irsn++) {
	416	if (TMath::Abs(pdg) == RsnPDG(irsn)) {
	417	// debug check
	418	if (mother >= 0) {
	419	TParticle *p = stack->Particle(mother);
	420	motherPDG = p->GetPdgCode();
	421	}
	422	part->ProductionVertex(vprod);
	423	dx = mcVertex[0] - vprod.X();
	424	dy = mcVertex[1] - vprod.Y();
	425	dz = mcVertex[2] - vprod.Z();
	426	dist = TMath::Sqrt(dxdx + dydy + dz*dz);
	427	AliDebugClass(1, Form("PDG = %d -- mother ID = %d, pdg = %d -- dist. from MC vertex = %f -- prod. time = %f", pdg, mother, motherPDG, dist, vprod.T()));
	428	// global counter is always filled
f4304ff3	429	fRsnYPt[0][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
035b8eff	430	// counter for primaries is filled only if mother is less than 0 (primary)
035b8eff	431	if (dist < fPrimaryThr) {
f4304ff3	432	fRsnYPt[1][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
035b8eff	433	}
	434	// fill the global histogram
	435	fProducedParticles->Fill(irsn + 1);
	436	// if one PDG match was found, no need to continue loop
	437	break;
	438	}
	439	}
	440	}
	441	}
	442	}
c865cb1d	443
	444	PostData(1, fOutputList);
	445	}
	446
	447	//________________________________________________________________________
	448	void AliAnalysisTaskResonanceQA::Terminate(Option_t *)
	449	{
	450	//
	451	// Termination
	452	// Quickly check that the output list is there
	453	//
	454
	455	fOutputList = dynamic_cast<TList*>(GetOutputData(1));
	456	if (!fOutputList) {
	457	AliError("Output list not found!!!");
	458	return;
	459	}
	460	}