[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDv0Monitor.cxx

/**************************************************************************
* 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-commercialf 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.                  *
**************************************************************************/

/////////////////////////////////////////////////////
//
//  Monitor V0 for TRD
//
//  Authors:                                          
//  Markus Heide <mheide@uni-muenster.de> 
//////////////////////////////////////////////////////
     
#include "TPDGCode.h"
#include "TFile.h"
#include "TTree.h"
#include "TEventList.h"
#include "TObjArray.h"
#include "TH2.h"
#include "TH2F.h"
#include "TH1I.h"
#include "TList.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TLegend.h"

#include "AliLog.h"
#include "AliESDtrack.h"
#include "AliPID.h"
#include "AliTRDv0Monitor.h"
#include "info/AliTRDv0Info.h"


ClassImp(AliTRDv0Monitor)

//________________________________________________________________________
AliTRDv0Monitor::AliTRDv0Monitor() 
  :AliTRDrecoTask()
   ,fOutput(NULL)
   ,fhQualityReductions(NULL)
   ,fV0s(NULL)
   ,fData(NULL)
   ,fInfo(NULL)
   ,fP(-1.) 
{
  //
  // Default constructor
  //
  SetNameTitle("v0Monitor", "Reference V0 Monitor");
}

//________________________________________________________________________
AliTRDv0Monitor::AliTRDv0Monitor(const char *name, const char *title) 
  :AliTRDrecoTask(name, title)
   ,fOutput(NULL)
   ,fhQualityReductions(NULL)
   ,fV0s(NULL)
   ,fData(NULL)
   ,fInfo(NULL)
   ,fP(-1.)
{
  //
  // Default constructor
  //
  DefineInput(2, TObjArray::Class()); // v0 list
  DefineInput(3, TObjArray::Class()); // pid info list 
  DefineOutput(1, TList::Class());
  
}


//________________________________________________________________________
AliTRDv0Monitor::~AliTRDv0Monitor() 
{
  if(fOutput) delete fOutput;
}
//____________________________________________________________________
void AliTRDv0Monitor::MakeSummary(){
  TCanvas *cOut = new TCanvas("v0MonitorSummary1", "Summary 1 for task V0Monitor", 1024, 768);
  cOut->cd();
  //GetRefFigure(4);
  cOut->SaveAs("V0MonitorSummary.gif");

  cOut = new TCanvas("v0MonitorSummary2","Summary 2 for task V0Monitor", 1024, 768);
  cOut->cd();
  //GetRefFigure(5);
  cOut->SaveAs("V0MonitorSummary2.gif");
}

//________________________________________________________________________
Bool_t AliTRDv0Monitor::GetRefFigure(Int_t /*ifig*/)
{
 
  AliInfo("Implementation on going ...");
  return kTRUE;
}

//________________________________________________________________________
void AliTRDv0Monitor::UserCreateOutputObjects()
{
  // Create histograms
  // Called once

 
 fOutput = new TList;
 fOutput->SetName("V0Monitoring");
 

 const char *samplename[AliPID::kSPECIES] = {"electrons","muons","pions","kaons","protons"};
 const char *decayname[AliTRDv0Info::kNDecays] = {"gamma","K0s","Lambda","AntiLambda"};
 const char *detectorname[kNDets] = {"ITS","TPC","TOF"};

 fhQualityReductions = new TH1I(Form("fhQualityReductions"),Form("Number of tracks cut out by different quality cut steps"),11,-9,2);
 fOutput->Add(fhQualityReductions);

 for(Int_t ipart = 0;ipart < AliPID::kSPECIES; ipart++){
   fhCutReductions[ipart] = new TH1I(Form("fhCutReductions_%s",samplename[ipart]),Form("Number of tracks cut out by different cut steps for %s",samplename[ipart]),19,-17,2);
   fOutput->Add(fhCutReductions[ipart]);
   for(Int_t idetector = 0; idetector < kNDets; idetector++){
     fhDetPID[idetector][ipart] = new TH2F(Form("fhDetector_%s_%s",detectorname[idetector],samplename[ipart]),Form("%s Likelihood for %s vs. momentum",detectorname[idetector], samplename[ipart]),100,0.,10.,100, 0., 1.);

     fOutput->Add(fhDetPID[idetector][ipart]);
   }
   fhComPID[ipart] = new TH2F(Form("fhComPID_%s",samplename[ipart]),Form("Combined TPC/TOF PID: Likelihood for %s",samplename[ipart]),100,0.,10.,100,0.,1.);

     fOutput->Add(fhComPID[ipart]);

     for(Int_t cutstep = 0; cutstep < kNCutSteps; cutstep++){
       fhTPCdEdx[ipart][cutstep] = new TH2F(Form("fhTPCdEdx_%s_[%d]",samplename[ipart],cutstep),Form("TPC dE/dx for %s [%d]",samplename[ipart],cutstep),100,0.,10.,300,0.,300.);
       
       fOutput->Add(fhTPCdEdx[ipart][cutstep]);
     }
 }

 for(Int_t iDecay = 0; iDecay < AliTRDv0Info::kNDecays; iDecay++){   
   for(Int_t cutstep =0; cutstep < kNCutSteps; cutstep++){
     fhV0Chi2ndf[iDecay][cutstep] =  new TH2F(Form("fhV0Chi2ndf_%s_[%d]",decayname[iDecay],cutstep),Form("Chi2/NDF vs. momentum"),100,0.,10.,500, 0., 500.);
     
     fOutput->Add(fhV0Chi2ndf[iDecay][cutstep]);
     
     fhPsiPair[iDecay][cutstep] =  new TH2F(Form("fhV0PsiPair_%s_[%d]",decayname[iDecay],cutstep),Form("Psi_pair vs. momentum"),100,0.,10.,200, 0., 1.6);
     
     fOutput->Add(fhPsiPair[iDecay][cutstep]);

     fhPointAngle[iDecay][cutstep] =  new TH2F(Form("fhPointAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Pointing Angle vs. momentum"),100,0.,10.,500, 0., 1.6);     
     fOutput->Add(fhPointAngle[iDecay][cutstep]);

     fhDCA[iDecay][cutstep] =  new TH2F(Form("fhDCA_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Daughter DCA vs. momentum"),100,0.,10.,500, 0., 1.);
     
     fOutput->Add(fhDCA[iDecay][cutstep]);

     fhOpenAngle[iDecay][cutstep] =  new TH2F(Form("fhOpenAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Opening Angle vs. momentum"),100,0.,10.,500, 0., 1.6);
     
     fOutput->Add(fhOpenAngle[iDecay][cutstep]);

     fhRadius[iDecay][cutstep] =  new TH2F(Form("fhRadius_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Generation Radius vs. momentum"),100,0.,10.,500, 0., 150.);
     
     fOutput->Add(fhRadius[iDecay][cutstep]);

    
   }

   fhInvMass[iDecay] =  new TH2F(Form("fhInvMass_%s",decayname[iDecay]),Form("Invariant Mass vs. momentum"),100,0.,10.,500, 0., 2.);
     
   fOutput->Add(fhInvMass[iDecay]); 
 
 } 
 
 /*TH1F *hV0mcPID[AliPID::kSPECIES][AliPID::kSPECIES];
   Int_t nPBins = 200;
   
   
   for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++){
   for(Int_t iSample = 0; iSample < AliPID::kSPECIES; iSample++){
   
   fhV0mcPID[iSample][iSpecies] = new TH1F(Form("fhV0mcPID_%s_is_%s",name[iSample],name[iSpecies]),Form("%s contained in %s sample",name[iSpecies],name[iSample]), nPBins, 0.2, 13.);
   }
   }*/
}
//________________________________________________________________________
void AliTRDv0Monitor::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))) return;
  if(!(fV0s    = dynamic_cast<TObjArray*>(GetInputData(2)))) return;
  if(!(fInfo   = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
  
  
  AliTRDtrackInfo     *track = NULL;
  AliTRDv0Info *v0(NULL);

 
  for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
    track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
    for(Int_t iv(0); iv<fV0s->GetEntriesFast(); iv++){
      if(!(v0 = (AliTRDv0Info*)fV0s->At(iv))) continue;
      if(!v0->HasTrack(track)) continue;
      ULong_t status = track->GetStatus();
      if(!(status&AliESDtrack::kTRDpid)) continue;
      
      fhQualityReductions->Fill(v0->fQuality);//fills integer codes for tracks cut out by track/V0 quality cuts
      
      if(!(v0->fQuality == 1))continue;

      for(Int_t part = 0; part < AliPID::kSPECIES; part++){
	fhCutReductions[part]->Fill(v0->GetPID(part,track));//fill in numbers of tracks eliminated by different PID cuts
      }

     
      for(Int_t idecay(0), part(-1); idecay <  AliTRDv0Info::kNDecays; idecay++){//loop over decay types considered for reference data
	
	if(idecay ==  AliTRDv0Info::kLambda){ //protons and pions from Lambda
	  part = AliPID::kProton;
	} else if(idecay == AliTRDv0Info::kAntiLambda) { //antiprotons and pions from Anti-Lambda     
	  part = AliPID::kProton;
	} else if(idecay ==   AliTRDv0Info::kK0s) {//pions from K0s
	  part = AliPID::kPion;
	} else if(idecay ==  AliTRDv0Info::kGamma) {//electrons from conversions
	  part = AliPID::kElectron;
	} 
	
	//fill histograms with track/V0 quality cuts only
	fhPsiPair[idecay][0]->Fill(v0->fV0Momentum,v0->fPsiPair);//Angle between daughter momentum plane and plane perpendicular to magnetic field
	fhInvMass[idecay]->Fill(v0->fV0Momentum,v0->fInvMass[idecay]);//Invariant mass
	fhPointAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fPointingAngle);// = TMath::ACos(esdv0->GetV0CosineOfPointingAngle()); // Cosine of pointing angle
	fhOpenAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fOpenAngle);// opening angle between daughters
	fhDCA[idecay][0]->Fill(v0->fV0Momentum,v0->fDCA);// Distance of closest approach of daughter tracks	
	fhV0Chi2ndf[idecay][0]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);//Kalman Filter Chi2/NDF
	fhRadius[idecay][0]->Fill(v0->fV0Momentum,v0->fRadius);//distance of decay/conversion from primary vertex in x-y plane
	
	if(v0->HasTrack(track) == -1){	  
	  fhTPCdEdx[part][0]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);//TPC dE/dx for negative track	 
	}
	else if(v0->HasTrack(track) == 1){
	  fhTPCdEdx[part][0]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);//TPC dE/dx for positive track
	}
	//fill histograms after invariant mass cuts
	if((v0->fInvMass[idecay] < v0->fUpInvMass[idecay][0])&&(v0->fInvMass[idecay]> v0->fDownInvMass[idecay])){
	  fhV0Chi2ndf[idecay][1]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
	  fhPsiPair[idecay][1]->Fill(v0->fV0Momentum,v0->fPsiPair);
	  fhPointAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fPointingAngle);
	  fhOpenAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fOpenAngle);
	  fhDCA[idecay][1]->Fill(v0->fV0Momentum,v0->fDCA);
	  fhRadius[idecay][1]->Fill(v0->fV0Momentum,v0->fRadius);
	  if(v0->HasTrack(track) == -1)
	    fhTPCdEdx[part][1]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
	  else if(v0->HasTrack(track) == 1)
	    fhTPCdEdx[part][1]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);

	}

	//fill histograms after all reference selection cuts
	if(v0->GetPID(part,track)==1){
	  fhV0Chi2ndf[idecay][2]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
	  fhPsiPair[idecay][2]->Fill(v0->fV0Momentum,v0->fPsiPair);
	  fhPointAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fPointingAngle);
	  fhOpenAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fOpenAngle);
	  fhDCA[idecay][2]->Fill(v0->fV0Momentum,v0->fDCA);
	  fhRadius[idecay][2]->Fill(v0->fV0Momentum,v0->fRadius);
	  if(v0->HasTrack(track) == -1)
	    fhTPCdEdx[part][2]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
	  else if(v0->HasTrack(track) == 1)
	    fhTPCdEdx[part][2]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);

	}
      }
    }
  }
  PostData(1, fOutput);
}
//_______________________________________________________________________

/*Int_t AliTRDv0Monitor::GetPDG(Int_t index)
{//yet to come


}*/
//_______________________________________________________________________
Commit	Line	Data
0ed6f095	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercialf purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/////////////////////////////////////////////////////
	17	//
	18	// Monitor V0 for TRD
	19	//
	20	// Authors:
	21	// Markus Heide <mheide@uni-muenster.de>
	22	//////////////////////////////////////////////////////
	23
	24	#include "TPDGCode.h"
	25	#include "TFile.h"
	26	#include "TTree.h"
	27	#include "TEventList.h"
	28	#include "TObjArray.h"
	29	#include "TH2.h"
	30	#include "TH2F.h"
	31	#include "TH1I.h"
	32	#include "TList.h"
fa6441b5	33	#include "TCanvas.h"
	34	#include "TPad.h"
	35	#include "TLegend.h"
0ed6f095	36
	37	#include "AliLog.h"
	38	#include "AliESDtrack.h"
	39	#include "AliPID.h"
	40	#include "AliTRDv0Monitor.h"
	41	#include "info/AliTRDv0Info.h"
	42
	43
	44	ClassImp(AliTRDv0Monitor)
	45
	46	//________________________________________________________________________
	47	AliTRDv0Monitor::AliTRDv0Monitor()
	48	:AliTRDrecoTask()
	49	,fOutput(NULL)
fa6441b5	50	,fhQualityReductions(NULL)
0ed6f095	51	,fV0s(NULL)
	52	,fData(NULL)
	53	,fInfo(NULL)
	54	,fP(-1.)
	55	{
	56	//
	57	// Default constructor
	58	//
	59	SetNameTitle("v0Monitor", "Reference V0 Monitor");
	60	}
	61
	62	//________________________________________________________________________
	63	AliTRDv0Monitor::AliTRDv0Monitor(const char name, const char title)
	64	:AliTRDrecoTask(name, title)
	65	,fOutput(NULL)
fa6441b5	66	,fhQualityReductions(NULL)
0ed6f095	67	,fV0s(NULL)
	68	,fData(NULL)
	69	,fInfo(NULL)
	70	,fP(-1.)
	71	{
	72	//
	73	// Default constructor
	74	//
	75	DefineInput(2, TObjArray::Class()); // v0 list
	76	DefineInput(3, TObjArray::Class()); // pid info list
	77	DefineOutput(1, TList::Class());
	78
	79	}
	80
	81
	82	//________________________________________________________________________
	83	AliTRDv0Monitor::~AliTRDv0Monitor()
	84	{
	85	if(fOutput) delete fOutput;
	86	}
fa6441b5	87	//____________________________________________________________________
	88	void AliTRDv0Monitor::MakeSummary(){
	89	TCanvas *cOut = new TCanvas("v0MonitorSummary1", "Summary 1 for task V0Monitor", 1024, 768);
	90	cOut->cd();
	91	//GetRefFigure(4);
	92	cOut->SaveAs("V0MonitorSummary.gif");
	93
	94	cOut = new TCanvas("v0MonitorSummary2","Summary 2 for task V0Monitor", 1024, 768);
	95	cOut->cd();
	96	//GetRefFigure(5);
	97	cOut->SaveAs("V0MonitorSummary2.gif");
	98	}
129a3fe2	99
	100	//________________________________________________________________________
	101	Bool_t AliTRDv0Monitor::GetRefFigure(Int_t /ifig/)
	102	{
fa6441b5	103
129a3fe2	104	AliInfo("Implementation on going ...");
	105	return kTRUE;
	106	}
	107
0ed6f095	108	//________________________________________________________________________
	109	void AliTRDv0Monitor::UserCreateOutputObjects()
	110	{
	111	// Create histograms
	112	// Called once
	113
	114
	115	fOutput = new TList;
	116	fOutput->SetName("V0Monitoring");
	117
	118
	119	const char *samplename[AliPID::kSPECIES] = {"electrons","muons","pions","kaons","protons"};
	120	const char *decayname[AliTRDv0Info::kNDecays] = {"gamma","K0s","Lambda","AntiLambda"};
	121	const char *detectorname[kNDets] = {"ITS","TPC","TOF"};
	122
fa6441b5	123	fhQualityReductions = new TH1I(Form("fhQualityReductions"),Form("Number of tracks cut out by different quality cut steps"),11,-9,2);
fa6441b5	124	fOutput->Add(fhQualityReductions);
0ed6f095	125
0ed6f095	126	for(Int_t ipart = 0;ipart < AliPID::kSPECIES; ipart++){
fa6441b5	127	fhCutReductions[ipart] = new TH1I(Form("fhCutReductions_%s",samplename[ipart]),Form("Number of tracks cut out by different cut steps for %s",samplename[ipart]),19,-17,2);
fa6441b5	128	fOutput->Add(fhCutReductions[ipart]);
0ed6f095	129	for(Int_t idetector = 0; idetector < kNDets; idetector++){
fa6441b5	130	fhDetPID[idetector][ipart] = new TH2F(Form("fhDetector_%s_%s",detectorname[idetector],samplename[ipart]),Form("%s Likelihood for %s vs. momentum",detectorname[idetector], samplename[ipart]),100,0.,10.,100, 0., 1.);
0ed6f095	131
fa6441b5	132	fOutput->Add(fhDetPID[idetector][ipart]);
0ed6f095	133	}
fa6441b5	134	fhComPID[ipart] = new TH2F(Form("fhComPID_%s",samplename[ipart]),Form("Combined TPC/TOF PID: Likelihood for %s",samplename[ipart]),100,0.,10.,100,0.,1.);
0ed6f095	135
fa6441b5	136	fOutput->Add(fhComPID[ipart]);
0ed6f095	137
0ed6f095	138	for(Int_t cutstep = 0; cutstep < kNCutSteps; cutstep++){
fa6441b5	139	fhTPCdEdx[ipart][cutstep] = new TH2F(Form("fhTPCdEdx_%s_[%d]",samplename[ipart],cutstep),Form("TPC dE/dx for %s [%d]",samplename[ipart],cutstep),100,0.,10.,300,0.,300.);
0ed6f095	140
fa6441b5	141	fOutput->Add(fhTPCdEdx[ipart][cutstep]);
0ed6f095	142	}
	143	}
	144
	145	for(Int_t iDecay = 0; iDecay < AliTRDv0Info::kNDecays; iDecay++){
	146	for(Int_t cutstep =0; cutstep < kNCutSteps; cutstep++){
fa6441b5	147	fhV0Chi2ndf[iDecay][cutstep] = new TH2F(Form("fhV0Chi2ndf_%s_[%d]",decayname[iDecay],cutstep),Form("Chi2/NDF vs. momentum"),100,0.,10.,500, 0., 500.);
0ed6f095	148
fa6441b5	149	fOutput->Add(fhV0Chi2ndf[iDecay][cutstep]);
0ed6f095	150
fa6441b5	151	fhPsiPair[iDecay][cutstep] = new TH2F(Form("fhV0PsiPair_%s_[%d]",decayname[iDecay],cutstep),Form("Psi_pair vs. momentum"),100,0.,10.,200, 0., 1.6);
0ed6f095	152
fa6441b5	153	fOutput->Add(fhPsiPair[iDecay][cutstep]);
0ed6f095	154
fa6441b5	155	fhPointAngle[iDecay][cutstep] = new TH2F(Form("fhPointAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Pointing Angle vs. momentum"),100,0.,10.,500, 0., 1.6);
fa6441b5	156	fOutput->Add(fhPointAngle[iDecay][cutstep]);
0ed6f095	157
fa6441b5	158	fhDCA[iDecay][cutstep] = new TH2F(Form("fhDCA_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Daughter DCA vs. momentum"),100,0.,10.,500, 0., 1.);
0ed6f095	159
fa6441b5	160	fOutput->Add(fhDCA[iDecay][cutstep]);
0ed6f095	161
fa6441b5	162	fhOpenAngle[iDecay][cutstep] = new TH2F(Form("fhOpenAngle_%s_[%d]",decayname[iDecay],cutstep),Form("Opening Angle vs. momentum"),100,0.,10.,500, 0., 1.6);
0ed6f095	163
fa6441b5	164	fOutput->Add(fhOpenAngle[iDecay][cutstep]);
0ed6f095	165
fa6441b5	166	fhRadius[iDecay][cutstep] = new TH2F(Form("fhRadius_%s_[%d]",decayname[iDecay],cutstep),Form("V0 Generation Radius vs. momentum"),100,0.,10.,500, 0., 150.);
0ed6f095	167
fa6441b5	168	fOutput->Add(fhRadius[iDecay][cutstep]);
0ed6f095	169
	170
	171
	172
	173	}
	174
fa6441b5	175	fhInvMass[iDecay] = new TH2F(Form("fhInvMass_%s",decayname[iDecay]),Form("Invariant Mass vs. momentum"),100,0.,10.,500, 0., 2.);
0ed6f095	176
fa6441b5	177	fOutput->Add(fhInvMass[iDecay]);
0ed6f095	178
	179	}
	180
	181	/TH1F hV0mcPID[AliPID::kSPECIES][AliPID::kSPECIES];
	182	Int_t nPBins = 200;
	183
	184
	185
	186	for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++){
	187	for(Int_t iSample = 0; iSample < AliPID::kSPECIES; iSample++){
	188
fa6441b5	189	fhV0mcPID[iSample][iSpecies] = new TH1F(Form("fhV0mcPID_%s_is_%s",name[iSample],name[iSpecies]),Form("%s contained in %s sample",name[iSpecies],name[iSample]), nPBins, 0.2, 13.);
0ed6f095	190	}
	191	}*/
	192	}
	193	//________________________________________________________________________
	194	void AliTRDv0Monitor::UserExec(Option_t *)
	195	{
	196	// Main loop
	197	// Called for each event
	198
	199	if(!(fTracks = dynamic_cast<TObjArray*>(GetInputData(1)))) return;
	200	if(!(fV0s = dynamic_cast<TObjArray*>(GetInputData(2)))) return;
	201	if(!(fInfo = dynamic_cast<TObjArray*>(GetInputData(3)))) return;
	202
	203
	204	AliTRDtrackInfo *track = NULL;
	205	AliTRDv0Info *v0(NULL);
	206
	207
	208
	209	for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
	210	track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
	211	for(Int_t iv(0); iv<fV0s->GetEntriesFast(); iv++){
	212	if(!(v0 = (AliTRDv0Info*)fV0s->At(iv))) continue;
	213	if(!v0->HasTrack(track)) continue;
	214	ULong_t status = track->GetStatus();
	215	if(!(status&AliESDtrack::kTRDpid)) continue;
	216
fa6441b5	217	fhQualityReductions->Fill(v0->fQuality);//fills integer codes for tracks cut out by track/V0 quality cuts
0ed6f095	218
	219	if(!(v0->fQuality == 1))continue;
	220
	221	for(Int_t part = 0; part < AliPID::kSPECIES; part++){
fa6441b5	222	fhCutReductions[part]->Fill(v0->GetPID(part,track));//fill in numbers of tracks eliminated by different PID cuts
0ed6f095	223	}
	224
	225
	226	for(Int_t idecay(0), part(-1); idecay < AliTRDv0Info::kNDecays; idecay++){//loop over decay types considered for reference data
	227
	228	if(idecay == AliTRDv0Info::kLambda){ //protons and pions from Lambda
	229	part = AliPID::kProton;
	230	} else if(idecay == AliTRDv0Info::kAntiLambda) { //antiprotons and pions from Anti-Lambda
	231	part = AliPID::kProton;
	232	} else if(idecay == AliTRDv0Info::kK0s) {//pions from K0s
	233	part = AliPID::kPion;
	234	} else if(idecay == AliTRDv0Info::kGamma) {//electrons from conversions
	235	part = AliPID::kElectron;
	236	}
	237
	238	//fill histograms with track/V0 quality cuts only
fa6441b5	239	fhPsiPair[idecay][0]->Fill(v0->fV0Momentum,v0->fPsiPair);//Angle between daughter momentum plane and plane perpendicular to magnetic field
	240	fhInvMass[idecay]->Fill(v0->fV0Momentum,v0->fInvMass[idecay]);//Invariant mass
	241	fhPointAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fPointingAngle);// = TMath::ACos(esdv0->GetV0CosineOfPointingAngle()); // Cosine of pointing angle
	242	fhOpenAngle[idecay][0]->Fill(v0->fV0Momentum,v0->fOpenAngle);// opening angle between daughters
	243	fhDCA[idecay][0]->Fill(v0->fV0Momentum,v0->fDCA);// Distance of closest approach of daughter tracks
	244	fhV0Chi2ndf[idecay][0]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);//Kalman Filter Chi2/NDF
	245	fhRadius[idecay][0]->Fill(v0->fV0Momentum,v0->fRadius);//distance of decay/conversion from primary vertex in x-y plane
0ed6f095	246
0ed6f095	247	if(v0->HasTrack(track) == -1){
fa6441b5	248	fhTPCdEdx[part][0]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);//TPC dE/dx for negative track
0ed6f095	249	}
0ed6f095	250	else if(v0->HasTrack(track) == 1){
fa6441b5	251	fhTPCdEdx[part][0]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);//TPC dE/dx for positive track
0ed6f095	252	}
	253	//fill histograms after invariant mass cuts
	254	if((v0->fInvMass[idecay] < v0->fUpInvMass[idecay][0])&&(v0->fInvMass[idecay]> v0->fDownInvMass[idecay])){
fa6441b5	255	fhV0Chi2ndf[idecay][1]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
	256	fhPsiPair[idecay][1]->Fill(v0->fV0Momentum,v0->fPsiPair);
	257	fhPointAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fPointingAngle);
	258	fhOpenAngle[idecay][1]->Fill(v0->fV0Momentum,v0->fOpenAngle);
	259	fhDCA[idecay][1]->Fill(v0->fV0Momentum,v0->fDCA);
	260	fhRadius[idecay][1]->Fill(v0->fV0Momentum,v0->fRadius);
0ed6f095	261	if(v0->HasTrack(track) == -1)
fa6441b5	262	fhTPCdEdx[part][1]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
0ed6f095	263	else if(v0->HasTrack(track) == 1)
fa6441b5	264	fhTPCdEdx[part][1]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);
0ed6f095	265
	266	}
	267
	268	//fill histograms after all reference selection cuts
	269	if(v0->GetPID(part,track)==1){
fa6441b5	270	fhV0Chi2ndf[idecay][2]->Fill(v0->fV0Momentum,v0->fChi2ndf[idecay]);
	271	fhPsiPair[idecay][2]->Fill(v0->fV0Momentum,v0->fPsiPair);
	272	fhPointAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fPointingAngle);
	273	fhOpenAngle[idecay][2]->Fill(v0->fV0Momentum,v0->fOpenAngle);
	274	fhDCA[idecay][2]->Fill(v0->fV0Momentum,v0->fDCA);
	275	fhRadius[idecay][2]->Fill(v0->fV0Momentum,v0->fRadius);
0ed6f095	276	if(v0->HasTrack(track) == -1)
fa6441b5	277	fhTPCdEdx[part][2]->Fill(v0->fTrackN->P(),v0->fTPCdEdx[AliTRDv0Info::kNeg]);
0ed6f095	278	else if(v0->HasTrack(track) == 1)
fa6441b5	279	fhTPCdEdx[part][2]->Fill(v0->fTrackP->P(),v0->fTPCdEdx[AliTRDv0Info::kPos]);
0ed6f095	280
	281	}
	282	}
	283	}
	284	}
	285	PostData(1, fOutput);
	286	}
	287	//_______________________________________________________________________
	288
	289	/*Int_t AliTRDv0Monitor::GetPDG(Int_t index)
	290	{//yet to come
	291
	292
	293
	294	}*/
	295	//_______________________________________________________________________
	296