[u/mrichter/AliRoot.git] / ITS / UPGRADE / testITSU / AliTaskITSUPerf.cxx

/*************************************************************************
* Copyright(c) 1998-2008, 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.                  * 
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
// Class AliTaskITSUPerf                                                 //
// Analysis task to produce data and MC histos needed for tracklets      //
// dNdEta extraction in multiple bins in one go                          //
// Author:  ruben.shahoyan@cern.ch                                       //
///////////////////////////////////////////////////////////////////////////

#include "TFile.h"
#include "TChain.h"
#include "TTree.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TH2F.h" 
#include "TList.h"
#include "TObjArray.h"
#include "TGeoGlobalMagField.h"

#include "AliAnalysisManager.h"

#include "AliESDEvent.h"  
#include "AliESDInputHandler.h"
#include "AliESDInputHandlerRP.h"
#include "AliCDBPath.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliCDBStorage.h"
#include "AliGeomManager.h"
#include "AliMagF.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
#include "AliStack.h"
#include "AliGenEventHeader.h"
#include "AliCentrality.h"
#include "AliLog.h"
#include "AliPhysicsSelection.h"
#include "AliGenEventHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliESDtrackCuts.h"

#include "AliTaskITSUPerf.h"
#include "../ITS/UPGRADE/AliITSURecoDet.h"
#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
#include "../ITS/UPGRADE/AliITSUClusterPix.h"
#include "../ITS/UPGRADE/AliITSUTrackCond.h"
#include <TGeoManager.h>


ClassImp(AliTaskITSUPerf)


const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = {
  "ITSokTPCok"
  ,"ITSokTPCfk"
  ,"ITSfkTPCok"
  ,"ITSfkTPCfk"
  ,"ITSTPCmismatch"
  ,"ITSTPCnoMatch"
};

//________________________________________________________________________
/*//Default constructor
AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
  : AliAnalysisTaskSE(name),
*/  
//________________________________________________________________________
AliTaskITSUPerf::AliTaskITSUPerf(const char *name) 
  : AliAnalysisTaskSE(name)
    //
  ,fOutput(0)
  ,fHistosCentMCLb(0)
  ,fHistosCent(0)
  ,fRPTree(0)
  ,fStack(0)
  ,fMCEvent(0)
  ,fVtxSPD(0)
  ,fVtxTrc(0)
  ,fESDEvent(0)
  ,fUseSpecialOutput(kTRUE)
  ,fUseMC(kTRUE)
  ,fTrackingCond(0)
    //
  ,fGeom(0)
  ,fITS(0)
  ,fNSelTracksMC(0)
  ,fMCStatus(0)
  ,fNPtBins(20)
  ,fNResBins(100)
  ,fPtMin(0)
  ,fPtMax(10.)
  ,fEtaMin(-3.)
  ,fEtaMax( 3.)
  ,fZVertexMin(-15.)
  ,fZVertexMax(15.)
    //
  ,fCurrCentBin(0)
  ,fNCentBins(1)
  ,fUseCentralityVar(0)
{
  // Constructor

  DefineOutput(1, TList::Class());
  //
}

//________________________________________________________________________
AliTaskITSUPerf::~AliTaskITSUPerf()
{
  // Destructor
  // histograms are in the output list and deleted when the output
  // list is deleted by the TSelector dtor
  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
    printf("Deleteing output\n");
    delete fOutput;
    fOutput = 0;
  }
  //
  fHistosCentMCLb.Clear();
  fHistosCent.Clear();
  //
  delete fITS;
  delete fGeom;
  //
}

//________________________________________________________________________
void AliTaskITSUPerf::UserCreateOutputObjects() 
{
  //
  if (fUseSpecialOutput) OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner(); 
  //
  // init geometry
  if (!gGeoManager) {
    TString geom = "geometry.root";
    AliGeomManager::LoadGeometry(geom.Data());
    if (!gGeoManager) AliFatal("Failed to load geometry");
  }
  //
  // create ITS interface
  fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE);
  AliITSUClusterPix::SetGeom(fGeom);
  fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface");
  fITS->CreateClusterArrays();
  //
  // Create histograms
  for (int ib=0;ib<fNCentBins;ib++) {
    BookHistos(ib);
  }
  //
  int nhist = fOutput->GetEntries();
  for (int i=0;i<nhist;i++) {
    TObject* hst = fOutput->At(i);
    if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
    ((TH1*)hst)->Sumw2();
  }
  //
  PostData(1, fOutput);
  //
}

//________________________________________________________________________
void AliTaskITSUPerf::UserExec(Option_t *) 
{
  // Main loop
  static int evID = 0;
  AliInfo(Form("Event: %d",evID++));
  //
  AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
  fRPTree = 0;
  AliESDInputHandlerRP *handRP = (AliESDInputHandlerRP*)anMan->GetInputEventHandler();
  if (!handRP) { AliFatal("No RP handler"); return; }
  //
  fESDEvent  = handRP->GetEvent();
  if (!fESDEvent) { AliFatal("No AliESDEvent"); return; }
  //
  fRPTree = handRP->GetTreeR("ITS");
  if (!fRPTree) { AliFatal("Invalid ITS cluster tree"); return; }
  //
  AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
  if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; }
  fMCEvent = eventHandler->MCEvent();
  if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; }
  fStack = fMCEvent->Stack();
  if (!fStack) { AliFatal("Stack is not available"); return; }
  //
  // MC Generator info
  AliGenEventHeader* mcGenH = 0;
  mcGenH = fMCEvent->GenEventHeader();
  TArrayF vtmc(3);
  mcGenH->PrimaryVertex(vtmc);
  for (int i=3;i--;) fVtxMC[i] = vtmc[i];
  //
  printf("MagField: %f\n",fESDEvent->GetMagneticField());
  printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]);
  //
  fVtxSPD = fESDEvent->GetPrimaryVertexSPD();
  //
  if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e | Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors());
  else printf("VtxSPD: N/A\n");
  //
  fVtxTrc = fESDEvent->GetPrimaryVertexTracks();
  if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e | Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors());
  else printf("VtxTRC: N/A\n");
  //
  BuildMCInfo();
  CheckTracks();
  //
}      


//________________________________________________________________________
void AliTaskITSUPerf::Terminate(Option_t *) 
{
  Printf("Terminating...");
}

//_________________________________________________________________________
void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb)
{
  // book standard set of histos for specific centrality bin and MC status type, 
  // adding them to output list and management array
  //
  const double kMaxDPt = 0.3;
  const double kMaxDCA = 0.3;

  TH2* h2=0;
  // pt resolution
  h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"),
		Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"),
		fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt);
  AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) );
  //
  // transverse DCA resolution
  h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"),
		Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"),
		fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
  AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) );
  //
  // Z DCA resolution
  h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"),
		Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"),
		fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
  AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) );
  //
}

//_________________________________________________________________________
void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin)
{
  // book standard set of histos for specific centrality bin
  //
  TH2* h2=0;
  //
  // MC labels combination vs pt
  h2 = new TH2F(Form("B%d_%s",bin,"MCLabvsPT"),
		Form("B%d %s",bin,"MCLabComb vs p_{T}"),
		fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
  AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatus,-1,bin) );
  for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
  //  

}

//_________________________________________________________________________
void AliTaskITSUPerf::BookHistos(Int_t bin)
{
  // book standard set of histos, adding them to output list and management array
  //
  // standard histos
  for (int ilb=0;ilb<kNLabelTypes;ilb++) {  // create similar set of histos for all kind of MC labels
    BookStandardHistosCentMCLb(bin,ilb);
  }
  BookStandardHistosCent(bin);              // st.histos regardless of MClabels comb.
}

//_________________________________________________________________________
void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at)
{
  // add single histo to the set
  if (at>=0) array->AddAtAndExpand(h,at);
  else       array->Add(h);
  fOutput->Add(h);
}

//_________________________________________________________________________
Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const
{
  // generate ID for the histo. htype must be <100, mcStat<kNLabelTypes
  //
  if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes));
  if (centBin>=fNCentBins)  AliFatal(Form("CentBin  ID=%d > max.allowed %d",centBin,fNCentBins));
  //
  if (centBin>=0) {
    if (mcStat>=0) { // MCtruthness dependent histos
      if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC));
      return htype + kHNStdHistosCentMC*(mcStat + kNLabelTypes*centBin);
    }
    else {
      if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent));
      return htype + kHNStdHistosCent*centBin;
    }
  }
  // custom histo
  return htype;
  //
}

//_________________________________________________________________________
void AliTaskITSUPerf::CheckTracks()
{
  // build mc truth info for tracks 
  //
  int ntr = fESDEvent->GetNumberOfTracks();
  int ntrMC = fStack->GetNtrack();
  //
  double fldBz = fESDEvent->GetMagneticField();
  float dcaRZ[2];
  //
  for (int itr=0;itr<ntr;itr++) {
    AliESDtrack* trc = fESDEvent->GetTrack(itr);
    //
    // at the moment we consider only TPC/ITS tracks
    if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
    Int_t labMC    = trc->GetLabel();
    Int_t labMCabs = TMath::Abs(labMC);
    Int_t labMCTPC = trc->GetTPCLabel();
    Int_t labMCITS = trc->GetITSLabel();
    Int_t nClTPC   = trc->GetNcls(1);
    Int_t nClITS   = trc->GetNcls(0);
    //
    UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs];
    Bool_t reject = mcStatus & BIT(kTrCondFail);
    if (reject) continue;
    //
    int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS);
    //
    double pt  = trc->Pt();
    double eta = trc->Eta();
    //
    //    printf("#%3d pt:%5.2f eta:%+5.2f | Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d | Ntpc:%3d Nits:%3d\n",itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS);
    //
    GetHisto(&fHistosCent,kHMatchStatus)->Fill(pt,mcLabType);  // matching status
    //
    AliMCParticle *part  = 0;
    //
    if (labMCabs<ntrMC) part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs);
    //
    if (part) {
      //
      double ptMC = part->Pt();
      double etaMC = part->Eta();
      //
      if ( (mcStatus&BIT(kMCPrimBit)) ) {
	// compare MC vs reco track params
	if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC);
	//
	// DCA resolution
	trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ);
	GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]);
	GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]);
	//
      }
      //
    }
    //


  }
  //

}

//_________________________________________________________________________
void AliTaskITSUPerf::BuildMCInfo()
{
  // build mc truth info for tracks 
  //
  int ntrMC = fStack->GetNtrack();
  //
  if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC);
  fMCStatus.Reset();
  //
  for (int itr=0;itr<ntrMC;itr++) {
    //
    AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
    if (!part->Charge()) continue;
    //
    UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
    mcStatus = 0;
    // 
    if (fStack->IsPhysicalPrimary(itr)) mcStatus |= BIT(kMCPrimBit);
    //
    /*
    TParticle* tp = part->Particle();
    printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print();
    */
    //
  }
  //
  int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0;
  // fill MC clusters info (if available)
  if (fRPTree) {
    // load clusters
    fITS->LoadClusters(fRPTree);
    //
    for (int ilr=fITS->GetNLayersActive();ilr--;) {
      AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
      for (int icl=lr->GetNClusters();icl--;) {
	AliCluster* cl = lr->GetCluster(icl); 
	// check labels
	for (int ilb=0;ilb<3;ilb++) {
	  Int_t lab = cl->GetLabel(ilb);
	  if (lab<0||lab>=ntrMC) continue;
	  UInt_t& mcStatus = (UInt_t &)fMCStatus[lab];
	  mcStatus |= 0x1<<(ilr+kITSHitBits);
	}
      }
    }
    //
    if (nTrCond) {
      for (int itr=0;itr<ntrMC;itr++) {
	//
	UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
	Bool_t fail = kTRUE;
	UShort_t patt = (mcStatus&0xffff);
	for (int ic=0;ic<nTrCond;ic++) {
	  AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic);
	  if (cnd->CheckPattern(patt)) {fail=kFALSE; break;}
	}
	if (fail) mcStatus |= BIT(kTrCondFail);
      }
    }
    //
  }
  //
}


//_________________________________________________________________________
Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS)
{
  // build type of track (0: both ITS and TPC correct, 1: ITS corr. TPC fake, 2: ITS fake TPC corr., 3: ITS fake TPC fake, 4: mismatch
  //
  int tp = 0;
  if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch;
  if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch;
  else {
    if   (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1;
    else              tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1;
  }
  return tp;
  //
}

//_________________________________________________________________________
Int_t AliTaskITSUPerf::GetCentralityBin() const
{
  // calculate centrality bin
  return 0;
}
Commit	Line	Data
fa10255f	1	/*************************************************************************
	2	* Copyright(c) 1998-2008, 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-commercial 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	// Class AliTaskITSUPerf //
	18	// Analysis task to produce data and MC histos needed for tracklets //
	19	// dNdEta extraction in multiple bins in one go //
	20	// Author: ruben.shahoyan@cern.ch //
	21	///////////////////////////////////////////////////////////////////////////
	22
	23	#include "TFile.h"
	24	#include "TChain.h"
	25	#include "TTree.h"
	26	#include "TRandom.h"
	27	#include "TH1F.h"
	28	#include "TH2F.h"
	29	#include "TList.h"
	30	#include "TObjArray.h"
	31	#include "TGeoGlobalMagField.h"
	32
	33	#include "AliAnalysisManager.h"
	34
	35	#include "AliESDEvent.h"
	36	#include "AliESDInputHandler.h"
	37	#include "AliESDInputHandlerRP.h"
	38	#include "AliCDBPath.h"
	39	#include "AliCDBManager.h"
	40	#include "AliCDBEntry.h"
	41	#include "AliCDBStorage.h"
	42	#include "AliGeomManager.h"
	43	#include "AliMagF.h"
	44	#include "AliMCEventHandler.h"
	45	#include "AliMCEvent.h"
	46	#include "AliMCParticle.h"
	47	#include "AliStack.h"
	48	#include "AliGenEventHeader.h"
	49	#include "AliCentrality.h"
	50	#include "AliLog.h"
	51	#include "AliPhysicsSelection.h"
	52	#include "AliGenEventHeader.h"
	53	#include "AliGenHijingEventHeader.h"
	54	#include "AliGenDPMjetEventHeader.h"
	55	#include "AliGenCocktailEventHeader.h"
	56	#include "AliESDtrackCuts.h"
	57
	58	#include "AliTaskITSUPerf.h"
	59	#include "../ITS/UPGRADE/AliITSURecoDet.h"
	60	#include "../ITS/UPGRADE/AliITSUGeomTGeo.h"
	61	#include "../ITS/UPGRADE/AliITSUClusterPix.h"
	62	#include "../ITS/UPGRADE/AliITSUTrackCond.h"
	63	#include <TGeoManager.h>
	64
65
66	ClassImp(AliTaskITSUPerf)
67
68
69	const char* AliTaskITSUPerf::fgkLabelTypes[AliTaskITSUPerf::kNLabelTypes] = {
70	"ITSokTPCok"
71	,"ITSokTPCfk"
72	,"ITSfkTPCok"
73	,"ITSfkTPCfk"
74	,"ITSTPCmismatch"
75	,"ITSTPCnoMatch"
76	};
77
78	//________________________________________________________________________
79	/*//Default constructor
80	AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
81	: AliAnalysisTaskSE(name),
82	*/
83	//________________________________________________________________________
84	AliTaskITSUPerf::AliTaskITSUPerf(const char *name)
85	: AliAnalysisTaskSE(name)
86	//
87	,fOutput(0)
88	,fHistosCentMCLb(0)
89	,fHistosCent(0)
90	,fRPTree(0)
91	,fStack(0)
92	,fMCEvent(0)
93	,fVtxSPD(0)
94	,fVtxTrc(0)
95	,fESDEvent(0)
96	,fUseSpecialOutput(kTRUE)
97	,fUseMC(kTRUE)
98	,fTrackingCond(0)
99	//
100	,fGeom(0)
101	,fITS(0)
102	,fNSelTracksMC(0)
103	,fMCStatus(0)
104	,fNPtBins(20)
105	,fNResBins(100)
106	,fPtMin(0)
107	,fPtMax(10.)
108	,fEtaMin(-3.)
109	,fEtaMax( 3.)
110	,fZVertexMin(-15.)
111	,fZVertexMax(15.)
112	//
113	,fCurrCentBin(0)
114	,fNCentBins(1)
115	,fUseCentralityVar(0)
116	{
117	// Constructor
118
119	DefineOutput(1, TList::Class());
120	//
121	}
122
123	//________________________________________________________________________
124	AliTaskITSUPerf::~AliTaskITSUPerf()
125	{
126	// Destructor
127	// histograms are in the output list and deleted when the output
128	// list is deleted by the TSelector dtor
129	if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
130	printf("Deleteing output\n");
131	delete fOutput;
132	fOutput = 0;
133	}
134	//
135	fHistosCentMCLb.Clear();
136	fHistosCent.Clear();
137	//
138	delete fITS;
139	delete fGeom;
140	//
141	}
142
143	//________________________________________________________________________
144	void AliTaskITSUPerf::UserCreateOutputObjects()
145	{
146	//
147	if (fUseSpecialOutput) OpenFile(1);
148	fOutput = new TList();
149	fOutput->SetOwner();
150	//
151	// init geometry
152	if (!gGeoManager) {
153	TString geom = "geometry.root";
154	AliGeomManager::LoadGeometry(geom.Data());
155	if (!gGeoManager) AliFatal("Failed to load geometry");
156	}
157	//
158	// create ITS interface
159	fGeom = new AliITSUGeomTGeo(kTRUE,kTRUE);
160	AliITSUClusterPix::SetGeom(fGeom);
161	fITS = new AliITSURecoDet(fGeom,"ITSURecoInterface");
162	fITS->CreateClusterArrays();
163	//
164	// Create histograms
165	for (int ib=0;ib<fNCentBins;ib++) {
166	BookHistos(ib);
167	}
168	//
169	int nhist = fOutput->GetEntries();
170	for (int i=0;i<nhist;i++) {
171	TObject* hst = fOutput->At(i);
172	if (!hst \|\| !(hst->InheritsFrom(TH1::Class()))) continue;
173	((TH1*)hst)->Sumw2();
174	}
175	//
176	PostData(1, fOutput);
177	//
178	}
179
180	//________________________________________________________________________
181	void AliTaskITSUPerf::UserExec(Option_t *)
182	{
183	// Main loop
184	static int evID = 0;
185	AliInfo(Form("Event: %d",evID++));
186	//
187	AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
188	fRPTree = 0;
189	AliESDInputHandlerRP handRP = (AliESDInputHandlerRP)anMan->GetInputEventHandler();
190	if (!handRP) { AliFatal("No RP handler"); return; }
191	//
192	fESDEvent = handRP->GetEvent();
193	if (!fESDEvent) { AliFatal("No AliESDEvent"); return; }
194	//
195	fRPTree = handRP->GetTreeR("ITS");
196	if (!fRPTree) { AliFatal("Invalid ITS cluster tree"); return; }
197	//
198	AliMCEventHandler* eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
199	if (!eventHandler) { AliFatal("Could not retrieve MC event handler"); return; }
200	fMCEvent = eventHandler->MCEvent();
201	if (!fMCEvent) { AliFatal("Could not retrieve MC event"); return; }
202	fStack = fMCEvent->Stack();
203	if (!fStack) { AliFatal("Stack is not available"); return; }
204	//
205	// MC Generator info
206	AliGenEventHeader* mcGenH = 0;
207	mcGenH = fMCEvent->GenEventHeader();
208	TArrayF vtmc(3);
209	mcGenH->PrimaryVertex(vtmc);
210	for (int i=3;i--;) fVtxMC[i] = vtmc[i];
211	//
212	printf("MagField: %f\n",fESDEvent->GetMagneticField());
213	printf("VtxMC : %+e %+e %+e\n",fVtxMC[0],fVtxMC[1],fVtxMC[2]);
214	//
215	fVtxSPD = fESDEvent->GetPrimaryVertexSPD();
216	//
217	if (fVtxSPD->GetStatus()) printf("VtxSPD: %+e %+e %+e \| Nc=%d\n",fVtxSPD->GetX(),fVtxSPD->GetY(),fVtxSPD->GetZ(),fVtxSPD->GetNContributors());
218	else printf("VtxSPD: N/A\n");
219	//
220	fVtxTrc = fESDEvent->GetPrimaryVertexTracks();
221	if (fVtxTrc->GetStatus()) printf("VtxTRC: %+e %+e %+e \| Nc=%d\n",fVtxTrc->GetX(),fVtxTrc->GetY(),fVtxTrc->GetZ(),fVtxTrc->GetNContributors());
222	else printf("VtxTRC: N/A\n");
223	//
224	BuildMCInfo();
225	CheckTracks();
226	//
227	}
228
229
230	//________________________________________________________________________
231	void AliTaskITSUPerf::Terminate(Option_t *)
232	{
233	Printf("Terminating...");
234	}
235
236	//_________________________________________________________________________
237	void AliTaskITSUPerf::BookStandardHistosCentMCLb(Int_t bin, Int_t mcLb)
238	{
239	// book standard set of histos for specific centrality bin and MC status type,
240	// adding them to output list and management array
241	//
242	const double kMaxDPt = 0.3;
243	const double kMaxDCA = 0.3;
244
245	TH2* h2=0;
246	// pt resolution
247	h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"dPT2PTvsPT"),
248	Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"#Deltap_{T}/p_{T} vs p_{T}"),
249	fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDPt,kMaxDPt);
250	AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResPTvsPTMC,mcLb,bin) );
251	//
252	// transverse DCA resolution
253	h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCARvsPT"),
254	Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA R vs p_{T}"),
255	fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
256	AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCARvsPTMC,mcLb,bin) );
257	//
258	// Z DCA resolution
259	h2 = new TH2F(Form("B%d_%s_%s",bin,fgkLabelTypes[mcLb],"DCAZvsPT"),
260	Form("B%d %s %s",bin,fgkLabelTypes[mcLb],"DCA Z vs p_{T}"),
261	fNPtBins,fPtMin,fPtMax,fNResBins,-kMaxDCA,kMaxDCA);
262	AddHisto(&fHistosCentMCLb,h2, GetHistoID(kHResDCAZvsPTMC,mcLb,bin) );
263	//
264	}
265
266	//_________________________________________________________________________
267	void AliTaskITSUPerf::BookStandardHistosCent(Int_t bin)
268	{
269	// book standard set of histos for specific centrality bin
270	//
271	TH2* h2=0;
272	//
273	// MC labels combination vs pt
274	h2 = new TH2F(Form("B%d_%s",bin,"MCLabvsPT"),
275	Form("B%d %s",bin,"MCLabComb vs p_{T}"),
276	fNPtBins,fPtMin,fPtMax,kNLabelTypes,-0.5,kNLabelTypes-0.5);
277	AddHisto(&fHistosCent,h2, GetHistoID(kHMatchStatus,-1,bin) );
278	for (int i=0;i<kNLabelTypes;i++) h2->GetYaxis()->SetBinLabel(i+1,fgkLabelTypes[i]);
279	//
280
281	}
282
283	//_________________________________________________________________________
284	void AliTaskITSUPerf::BookHistos(Int_t bin)
285	{
286	// book standard set of histos, adding them to output list and management array
287	//
288	// standard histos
289	for (int ilb=0;ilb<kNLabelTypes;ilb++) { // create similar set of histos for all kind of MC labels
290	BookStandardHistosCentMCLb(bin,ilb);
291	}
292	BookStandardHistosCent(bin); // st.histos regardless of MClabels comb.
293	}
294
295	//_________________________________________________________________________
296	void AliTaskITSUPerf::AddHisto(TObjArray* array,TObject* h, Int_t at)
297	{
298	// add single histo to the set
299	if (at>=0) array->AddAtAndExpand(h,at);
300	else array->Add(h);
301	fOutput->Add(h);
302	}
303
304	//_________________________________________________________________________
305	Int_t AliTaskITSUPerf::GetHistoID(Int_t htype, Int_t mcStat, Int_t centBin) const
306	{
307	// generate ID for the histo. htype must be <100, mcStat<kNLabelTypes
308	//
309	if (mcStat>=kNLabelTypes) AliFatal(Form("MCStatus ID=%d > max.allowed %d",mcStat,kNLabelTypes));
310	if (centBin>=fNCentBins) AliFatal(Form("CentBin ID=%d > max.allowed %d",centBin,fNCentBins));
311	//
312	if (centBin>=0) {
313	if (mcStat>=0) { // MCtruthness dependent histos
314	if (htype>=kHNStdHistosCentMC) AliFatal(Form("StdCentMC histo ID=%d > max.allowed %d",htype,kHNStdHistosCentMC));
315	return htype + kHNStdHistosCentMC(mcStat + kNLabelTypescentBin);
316	}
317	else {
318	if (htype>=kHNStdHistosCent) AliFatal(Form("StdCent histo ID=%d > max.allowed %d",htype,kHNStdHistosCent));
319	return htype + kHNStdHistosCent*centBin;
320	}
321	}
322	// custom histo
323	return htype;
324	//
325	}
326
327	//_________________________________________________________________________
328	void AliTaskITSUPerf::CheckTracks()
329	{
330	// build mc truth info for tracks
331	//
332	int ntr = fESDEvent->GetNumberOfTracks();
333	int ntrMC = fStack->GetNtrack();
334	//
335	double fldBz = fESDEvent->GetMagneticField();
336	float dcaRZ[2];
337	//
338	for (int itr=0;itr<ntr;itr++) {
339	AliESDtrack* trc = fESDEvent->GetTrack(itr);
340	//
341	// at the moment we consider only TPC/ITS tracks
342	if (!trc->IsOn(AliESDtrack::kTPCin)) continue;
343	Int_t labMC = trc->GetLabel();
344	Int_t labMCabs = TMath::Abs(labMC);
345	Int_t labMCTPC = trc->GetTPCLabel();
346	Int_t labMCITS = trc->GetITSLabel();
347	Int_t nClTPC = trc->GetNcls(1);
348	Int_t nClITS = trc->GetNcls(0);
349	//
350	UInt_t& mcStatus = (UInt_t &)fMCStatus[labMCabs];
351	Bool_t reject = mcStatus & BIT(kTrCondFail);
352	if (reject) continue;
353	//
354	int mcLabType = GetMCLabType(labMCTPC,labMCITS,nClTPC,nClITS);
355	//
356	double pt = trc->Pt();
357	double eta = trc->Eta();
358	//
359	// printf("#%3d pt:%5.2f eta:%+5.2f \| Lb:%+6d LbTPC:%+6d LbITS:%+6d MCTp:%d \| Ntpc:%3d Nits:%3d\n",itr,pt,eta,labMC,labMCTPC,labMCITS,mcLabType, nClTPC,nClITS);
360	//
361	GetHisto(&fHistosCent,kHMatchStatus)->Fill(pt,mcLabType); // matching status
362	//
363	AliMCParticle *part = 0;
364	//
365	if (labMCabs<ntrMC) part = (AliMCParticle*)fMCEvent->GetTrack(labMCabs);
366	//
367	if (part) {
368	//
369	double ptMC = part->Pt();
370	double etaMC = part->Eta();
371	//
372	if ( (mcStatus&BIT(kMCPrimBit)) ) {
373	// compare MC vs reco track params
374	if (ptMC>0) GetHisto(&fHistosCentMCLb,kHResPTvsPTMC, mcLabType)->Fill(ptMC, (ptMC-pt)/ptMC);
375	//
376	// DCA resolution
377	trc->GetDZ(fVtxMC[0],fVtxMC[1],fVtxMC[2], fldBz, dcaRZ);
378	GetHisto(&fHistosCentMCLb,kHResDCARvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[0]);
379	GetHisto(&fHistosCentMCLb,kHResDCAZvsPTMC, mcLabType)->Fill(ptMC, dcaRZ[1]);
380	//
381	}
382	//
383	}
384	//
385
386
387	}
388	//
389
390	}
391
392	//_________________________________________________________________________
393	void AliTaskITSUPerf::BuildMCInfo()
394	{
395	// build mc truth info for tracks
396	//
397	int ntrMC = fStack->GetNtrack();
398	//
399	if (fMCStatus.GetSize()<ntrMC) fMCStatus.Set(ntrMC);
400	fMCStatus.Reset();
401	//
402	for (int itr=0;itr<ntrMC;itr++) {
403	//
404	AliMCParticle part = (AliMCParticle)fMCEvent->GetTrack(itr);
405	if (!part->Charge()) continue;
406	//
407	UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
408	mcStatus = 0;
409	//
410	if (fStack->IsPhysicalPrimary(itr)) mcStatus \|= BIT(kMCPrimBit);
411	//
412	/*
413	TParticle* tp = part->Particle();
414	printf("mc%4d IsPri=%d",itr,fStack->IsPhysicalPrimary(itr)); tp->Print();
415	*/
416	//
417	}
418	//
419	int nTrCond = fTrackingCond ? fTrackingCond->GetEntriesFast() : 0;
420	// fill MC clusters info (if available)
421	if (fRPTree) {
422	// load clusters
423	fITS->LoadClusters(fRPTree);
424	//
425	for (int ilr=fITS->GetNLayersActive();ilr--;) {
426	AliITSURecoLayer* lr = fITS->GetLayerActive(ilr);
427	for (int icl=lr->GetNClusters();icl--;) {
428	AliCluster* cl = lr->GetCluster(icl);
429	// check labels
430	for (int ilb=0;ilb<3;ilb++) {
431	Int_t lab = cl->GetLabel(ilb);
432	if (lab<0\|\|lab>=ntrMC) continue;
433	UInt_t& mcStatus = (UInt_t &)fMCStatus[lab];
434	mcStatus \|= 0x1<<(ilr+kITSHitBits);
435	}
436	}
437	}
438	//
439	if (nTrCond) {
440	for (int itr=0;itr<ntrMC;itr++) {
441	//
442	UInt_t& mcStatus = (UInt_t &)fMCStatus[itr];
443	Bool_t fail = kTRUE;
444	UShort_t patt = (mcStatus&0xffff);
445	for (int ic=0;ic<nTrCond;ic++) {
446	AliITSUTrackCond* cnd = (AliITSUTrackCond*) fTrackingCond->At(ic);
447	if (cnd->CheckPattern(patt)) {fail=kFALSE; break;}
448	}
449	if (fail) mcStatus \|= BIT(kTrCondFail);
450	}
451	}
452	//
453	}
454	//
455	}
456
457
458	//_________________________________________________________________________
459	Int_t AliTaskITSUPerf::GetMCLabType(Int_t labMCTPC,Int_t labMCITS, Int_t nClTPC, Int_t nClITS)
460	{
461	// build type of track (0: both ITS and TPC correct, 1: ITS corr. TPC fake, 2: ITS fake TPC corr., 3: ITS fake TPC fake, 4: mismatch
462	//
463	int tp = 0;
464	if (nClTPC>0 && nClITS<1) return kITSTPCNoMatch;
465	if (TMath::Abs(labMCTPC)!=TMath::Abs(labMCITS)) tp = kITSTPCMismatch;
466	else {
467	if (labMCITS<0) tp = labMCTPC<0 ? kITS0TPC0 : kITS0TPC1;
468	else tp = labMCTPC<0 ? kITS1TPC0 : kITS1TPC1;
469	}
470	return tp;
471	//
472	}
473
474	//_________________________________________________________________________
475	Int_t AliTaskITSUPerf::GetCentralityBin() const
476	{
477	// calculate centrality bin
478	return 0;
479	}
480