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

#include <Riostream.h>

#include <TChain.h>
#include <TH1.h>
#include <TH2.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


ClassImp(AliAnalysisTaskResonanceQA)

//________________________________________________________________________
AliAnalysisTaskResonanceQA::AliAnalysisTaskResonanceQA(const char *name) :
   AliAnalysisTaskSE(name), 
   fT0(AliESDpid::kTOF_T0),
   fPrimaryThr(1E-6),
   fVz(10.0),
   fOutputList(0),
   fSelectedEvts(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::RsnName(ERsn type)
{
//
// 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 kXiStarM   : return "XiStarM";
      case kSigmaStarP: return "SigmaStarP";
      case kSigmaStar0: return "SigmaStar0";
      case kSigmaStarM: return "SigmaStarM";
      case kDeltaPP   : return "DeltaPP";
      default         : return "Unknown";
   }
}

//________________________________________________________________________
const char* AliAnalysisTaskResonanceQA::RsnSymbol(ERsn type)
{
//
// 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 kXiStarM   : return "#Xi^{*-}";
      case kSigmaStarP: return "#Sigma^{*+}";
      case kSigmaStar0: return "#Sigma^{*0}";
      case kSigmaStarM: return "#Sigma^{*-}";
      case kDeltaPP   : return "#Delta^{++}";
      default         : return "Unknown";
   }
}

//________________________________________________________________________
Int_t AliAnalysisTaskResonanceQA::RsnPDG(ERsn type)
{
//
// 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 kXiStarM   : return 3314;
      case kSigmaStarP: return 3224;
      case kSigmaStar0: return 3214;
      case kSigmaStarM: return 3114;
      case kDeltaPP   : return 2224;
      default         : return    0;
   }
}

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

   Int_t i;
   Int_t ESDdEdxBin = 1000;  Double_t ESDdEdxMin = 0;  Double_t ESDdEdxMax = 1000;
   Int_t ESDQPBin = 900;     Double_t ESDQPMin = -4.5; Double_t ESDQPMax = 4.5;
   Int_t PtBin = 500;        Double_t PtMin = 0.;      Double_t PtMax = 10.;
   Int_t YBin = 600;         Double_t YMin = -15.;     Double_t YMax = 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", 3, 0, 3);
   fdEdxTPC = new TH2F("fdEdxTPC", "dEdxTPC, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", ESDQPBin, ESDQPMin, ESDQPMax, ESDdEdxBin, ESDdEdxMin, ESDdEdxMax);
   fdEdxITS = new TH2F("fdEdxITS", "dEdxITS, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", ESDQPBin, ESDQPMin, ESDQPMax, ESDdEdxBin, ESDdEdxMin, ESDdEdxMax);
   fTOFpid = new TH2F("fTOFpid", "TOF PID;p (GeV/c);#beta;", ESDQPBin, ESDQPMin, ESDQPMax, 300, 0.1, 1.2);
   
   // output histograms for track quality
   fDCAXYvsPtBeforeCuts = new TH2F("DCAXYvsPtBeforeCuts", "DCAXYvsPtBeforeCuts;p_{t} (GeV/c);DCAXY", PtBin, PtMin, PtMax, 500, -10.0, 10.0);
   fDCAZvsPtBeforeCuts = new TH2F("DCAZvsPtBeforeCuts", "DCAZvsPtBeforeCuts;p_{t} (GeV/c);DCAZ", PtBin, PtMin, PtMax, 500, -10.0, 10.0);
   fNClusterPtBeforeCuts = new TH2F("NClusterPtBeforeCuts", "NClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", PtBin, PtMin, PtMax, 180, 0., 180.);
   fNFindableClusterPtBeforeCuts = new TH2F("NFindableClusterPtBeforeCuts", "NFindableClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", PtBin, PtMin, PtMax, 180, 0., 180.);
   fNCrossedRowsTPCPtBeforeCuts = new TH2F("NCrossedRowsTPCPtBeforeCuts", "NCrossedRowsTPCPtBeforeCuts;p_{t} (GeV/c);NCrossedRowsTPC", PtBin, PtMin, PtMax, 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 TH2F(Form("%s_all", RsnName(i)), Form("%s -- ALL;p_{t} (GeV/c);Rapidity", RsnName(i)), PtBin, PtMin, PtMax, YBin, YMin, YMax);
      fRsnYPt[1][i] = new TH2F(Form("%s_prim", RsnName(i)), Form("%s -- PRIMARY;p_{t} (GeV/c);Rapidity", RsnName(i)), PtBin, PtMin, PtMax, YBin, YMin, YMax);
      fProducedParticles->GetXaxis()->SetBinLabel(i + 1, RsnSymbol(i));
   }
   
   // add histogram to list
   fOutputList->Add(fSelectedEvts);
   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 selected events
   fSelectedEvts->Fill(1.1);
   
   // 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());
                  // counter for primaries is filled only if mother is less than 0 (primary)
                  if (dist < fPrimaryThr) fRsnYPt[1][irsn]->Fill(part->Pt(), part->Y());
                  // fill the global histogram
                  fProducedParticles->Fill(irsn + 1);
                  // if one PDG match was found, no need to continue loop
                  break;
               }
            }
         }
      }
   }
   
   // try to cast input event to ESD and loop on tracks
   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
   if (fESD) {
      // check primary vertex:
      // we accept only tracks/SPD with at least 1 contributor
      const AliESDVertex *v0 = fESD->GetPrimaryVertexTracks();
      if (!v0) return;
      if (v0->GetNContributors() < 1) {
         v0 = fESD->GetPrimaryVertexSPD();
         if (!v0) return;
         if (v0->GetNContributors() < 1) return;
      }
      if (!v0) return;
      if (TMath::Abs(v0->GetZv()) > fVz) return;
   
      // 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

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