[u/mrichter/AliRoot.git] / PWG1 / TPC / AliPerformanceTPC.cxx

//------------------------------------------------------------------------------
// Implementation of AliPerformanceTPC class. It keeps information from 
// comparison of reconstructed and MC particle tracks. In addtion, 
// it keeps selection cuts used during comparison. The comparison 
// information is stored in the ROOT histograms. Analysis of these 
// histograms can be done by using Analyse() class function. The result of 
// the analysis (histograms/graphs) are stored in the folder which is
// a data member of AliPerformanceTPC.
//
// Author: J.Otwinowski 04/02/2008 
//------------------------------------------------------------------------------

/*
 
  // after running comparison task, read the file, and get component
  gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  LoadMyLibs();

  TFile f("Output.root");
  AliPerformanceTPC * compObj = (AliPerformanceTPC*)coutput->FindObject("AliPerformanceTPC");
 
  // analyse comparison data
  compObj->Analyse();

  // the output histograms/graphs will be stored in the folder "folderTPC" 
  compObj->GetAnalysisFolder()->ls("*");

  // user can save whole comparison object (or only folder with anlysed histograms) 
  // in the seperate output file (e.g.)
  TFile fout("Analysed_TPC.root","recreate");
  compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  fout.Close();

*/

#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TAxis.h"
#include "TPostScript.h"

#include "AliPerformanceTPC.h" 
#include "AliESDEvent.h" 
#include "AliESDVertex.h"
#include "AliESDtrack.h"
#include "AliLog.h" 
#include "AliMCEvent.h" 
#include "AliHeader.h" 
#include "AliGenEventHeader.h" 
#include "AliStack.h" 
#include "AliMCInfoCuts.h" 
#include "AliRecInfoCuts.h" 
#include "AliTracker.h" 
#include "AliTreeDraw.h" 
#include "AliTPCTransform.h" 
#include "AliTPCseed.h" 
#include "AliTPCcalibDB.h" 
#include "AliESDfriend.h" 
#include "AliESDfriendTrack.h" 
#include "AliTPCclusterMI.h" 

using namespace std;

ClassImp(AliPerformanceTPC)

//_____________________________________________________________________________
AliPerformanceTPC::AliPerformanceTPC():
  AliPerformanceObject("AliPerformanceTPC"),
  fTPCClustHisto(0),
  fTPCEventHisto(0),
  fTPCTrackHisto(0),

  // Cuts 
  fCutsRC(0),  
  fCutsMC(0),  

  // histogram folder 
  fAnalysisFolder(0),

  // select track vertex
  fUseTrackVertex(kFALSE) 
{
  Init();
}

//_____________________________________________________________________________
AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  AliPerformanceObject(name,title),
  fTPCClustHisto(0),
  fTPCEventHisto(0),
  fTPCTrackHisto(0),

  // Cuts 
  fCutsRC(0),  
  fCutsMC(0),  

  // histogram folder 
  fAnalysisFolder(0),

  // select track vertex
  fUseTrackVertex(kFALSE) 
{
  // named constructor	
  // 
  SetAnalysisMode(analysisMode);
  SetHptGenerator(hptGenerator);

  Init();
}

//_____________________________________________________________________________
AliPerformanceTPC::~AliPerformanceTPC()
{
  // destructor
   
  if(fTPCClustHisto) delete fTPCClustHisto; fTPCClustHisto=0;     
  if(fTPCEventHisto) delete fTPCEventHisto; fTPCEventHisto=0;     
  if(fTPCTrackHisto) delete fTPCTrackHisto; fTPCTrackHisto=0;     
  if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
}

//_____________________________________________________________________________
void AliPerformanceTPC::Init(){
  //
  // histogram bining
  //

  // set pt bins
  Int_t nPtBins = 50;
  Double_t ptMin = 1.e-2, ptMax = 10.;

  Double_t *binsPt = 0;
  if (IsHptGenerator())  { 
    nPtBins = 100; ptMax = 100.;
    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
  } else {
    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
  }

  /*
  Int_t nPtBins = 31;
  Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
  Double_t ptMin = 0., ptMax = 10.; 

  if(IsHptGenerator() == kTRUE) {
    nPtBins = 100;
    ptMin = 0.; ptMax = 100.; 
  }
  */
  // 
  //gclX:gclY:TPCSide
  Int_t binsTPCClustHisto[3]=  { 500,   500,  2 };
  Double_t minTPCClustHisto[3]={-250., -250., 0.};
  Double_t maxTPCClustHisto[3]={ 250.,  250., 2.};

  fTPCClustHisto = new THnSparseF("fTPCClustHisto","gclX:gclY:TPCSide",3,binsTPCClustHisto,minTPCClustHisto,maxTPCClustHisto);
  fTPCClustHisto->GetAxis(0)->SetTitle("gclX (cm)");
  fTPCClustHisto->GetAxis(1)->SetTitle("gclY (cm)");
  fTPCClustHisto->GetAxis(2)->SetTitle("TPCSide");
  fTPCClustHisto->Sumw2();
 

  // Xv:Yv:Zv:mult:multP:multN:vertStatus
  Int_t binsTPCEventHisto[7]=  {100,  100,   100,  151,   151,   151, 2   };
  Double_t minTPCEventHisto[7]={-10., -10., -30.,  -0.5,  -0.5,  -0.5, 0.  };
  Double_t maxTPCEventHisto[7]={ 10.,  10.,  30.,  150.5, 150.5, 150.5, 2. };

  fTPCEventHisto = new THnSparseF("fTPCEventHisto","Xv:Yv:Zv:mult:multP:multN:vertStatus",7,binsTPCEventHisto,minTPCEventHisto,maxTPCEventHisto);
  fTPCEventHisto->GetAxis(0)->SetTitle("Xv (cm)");
  fTPCEventHisto->GetAxis(1)->SetTitle("Yv (cm)");
  fTPCEventHisto->GetAxis(2)->SetTitle("Zv (cm)");
  fTPCEventHisto->GetAxis(3)->SetTitle("mult");
  fTPCEventHisto->GetAxis(4)->SetTitle("multP");
  fTPCEventHisto->GetAxis(5)->SetTitle("multN");
  fTPCEventHisto->GetAxis(6)->SetTitle("vertStatus");
  fTPCEventHisto->Sumw2();


  // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:DCAr:DCAz:eta:phi:pt:charge
  Int_t binsTPCTrackHisto[9]=  { 160,  50,  60,  100, 100,  30,   144,             nPtBins,  3 };
  Double_t minTPCTrackHisto[9]={ 0.,   0.,  0., -10,  -10., -1.5, 0.,             ptMin,   -1.5};
  Double_t maxTPCTrackHisto[9]={ 160., 10., 1.2, 10,   10.,  1.5, 2.*TMath::Pi(), ptMax,    1.5};

  fTPCTrackHisto = new THnSparseF("fTPCTrackHisto","nClust:chi2PerClust:nClust/nFindableClust:DCAr:DCAz:eta:phi:pt:charge",9,binsTPCTrackHisto,minTPCTrackHisto,maxTPCTrackHisto);
  fTPCTrackHisto->SetBinEdges(7,binsPt);

  fTPCTrackHisto->GetAxis(0)->SetTitle("nClust");
  fTPCTrackHisto->GetAxis(1)->SetTitle("chi2PerClust");
  fTPCTrackHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
  fTPCTrackHisto->GetAxis(3)->SetTitle("DCAr (cm)");
  fTPCTrackHisto->GetAxis(4)->SetTitle("DCAz (cm)");
  fTPCTrackHisto->GetAxis(5)->SetTitle("#eta");
  fTPCTrackHisto->GetAxis(6)->SetTitle("#phi (rad)");
  fTPCTrackHisto->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
  fTPCTrackHisto->GetAxis(8)->SetTitle("charge");
  fTPCTrackHisto->Sumw2();

  // Init cuts 
  if(!fCutsMC) 
    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
  if(!fCutsRC) 
    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");

  // init folder
  fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
}

//_____________________________________________________________________________
void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
{
  if(!esdTrack) return;

  // Fill TPC only resolution comparison information 
  const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
  if(!track) return;

  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
  esdTrack->GetImpactParametersTPC(dca,cov);

  Float_t q = esdTrack->Charge();
  Float_t pt = track->Pt();
  Float_t eta = track->Eta();
  Float_t phi = track->Phi();
  Int_t nClust = esdTrack->GetTPCclusters(0);
  Int_t nFindableClust = esdTrack->GetTPCNclsF();

  Float_t chi2PerCluster = 0.;
  if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);

  Float_t clustPerFindClust = 0.;
  if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
  
  //
  // select primaries
  //
  Double_t dcaToVertex = -1;
  if( fCutsRC->GetDCAToVertex2D() ) 
  {
      dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ()); 
  }
  if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
  if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
  if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;

  Double_t vTPCTrackHisto[9] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,q};
  fTPCTrackHisto->Fill(vTPCTrackHisto); 
 
  //
  // Fill rec vs MC information
  //
  if(!stack) return;

}

//_____________________________________________________________________________
void AliPerformanceTPC::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
{
  // Fill comparison information (TPC+ITS) 
  AliDebug(AliLog::kWarning, "Warning: Not implemented");
}
 
//_____________________________________________________________________________
void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
{
  // Fill comparison information (constarained parameters) 
  AliDebug(AliLog::kWarning, "Warning: Not implemented");
}
 
//_____________________________________________________________________________
void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
{
  // Process comparison information 
  //
  if(!esdEvent) 
  {
    Error("Exec","esdEvent not available");
    return;
  }
  AliHeader* header = 0;
  AliGenEventHeader* genHeader = 0;
  AliStack* stack = 0;
  TArrayF vtxMC(3);
  
  if(bUseMC)
  {
    if(!mcEvent) {
      Error("Exec","mcEvent not available");
      return;
    }
    // get MC event header
    header = mcEvent->Header();
    if (!header) {
      Error("Exec","Header not available");
      return;
    }
    // MC particle stack
    stack = mcEvent->Stack();
    if (!stack) {
      Error("Exec","Stack not available");
      return;
    }
    // get MC vertex
    genHeader = header->GenEventHeader();
    if (!genHeader) {
      Error("Exec","Could not retrieve genHeader from Header");
      return;
    }
    genHeader->PrimaryVertex(vtxMC);
  } 
  
  // trigger
  if(!bUseMC &&GetTriggerClass()) {
    Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
    if(!isEventTriggered) return; 
  }

  // get TPC event vertex
  const AliESDVertex *vtxESD = NULL; 
  if(fUseTrackVertex) {
    vtxESD = esdEvent->GetPrimaryVertexTracks();
  } else {
    vtxESD = esdEvent->GetPrimaryVertexTPC();
  }
  if(!vtxESD) return;

  //  events with rec. vertex
  Int_t mult=0; Int_t multP=0; Int_t multN=0;
  if(vtxESD->GetStatus() >0)
  {
  //  Process ESD events
  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
  { 
    AliESDtrack *track = esdEvent->GetTrack(iTrack);
    if(!track) continue;

    if(bUseESDfriend && esdFriend) 
    {
      AliESDfriendTrack *friendTrack=esdFriend->GetTrack(iTrack);
      if(friendTrack) 
      {
        TObject *calibObject=0;
        AliTPCseed *seed=0;
        if (!friendTrack) continue;
          for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j) {
	    if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) {
	    break;
	  }
        }

        //AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
	for (Int_t irow=0;irow<159;irow++) {
	if(!seed) continue;
	  
	  AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
	  if (!cluster) continue;

	      Float_t gclf[3];
	      cluster->GetGlobalXYZ(gclf);

	      //Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
	      //Int_t i[1]={cluster->GetDetector()};
              //transform->Transform(x,i,0,1);
	      //printf("gx %f gy  %f  gz %f \n", cluster->GetX(), cluster->GetY(),cluster->GetZ());
	      //printf("gclf[0] %f gclf[1]  %f  gclf[2] %f \n", gclf[0], gclf[1],  gclf[2]);
     
             Int_t TPCside; 
	     if(gclf[2]>0.) TPCside=0; // A side 
	     else TPCside=1;

             Double_t vTPCClust[3] = { gclf[0], gclf[1],  TPCside };
             fTPCClustHisto->Fill(vTPCClust);
         }
       }
    }

    if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
    else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
    else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
    else {
      printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
      return;
    }

   // TPC only
   AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
   if(!tpcTrack) continue;

   // track selection
   if( fCutsRC->AcceptTrack(tpcTrack) ) { 
     mult++;
     if(tpcTrack->Charge()>0.) multP++;
     if(tpcTrack->Charge()<0.) multN++;
   }

   if(tpcTrack) delete tpcTrack;
  }
  }
  //
  
  Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
  fTPCEventHisto->Fill(vTPCEvent);
}

//_____________________________________________________________________________
void AliPerformanceTPC::Analyse() {
  //
  // Analyse comparison information and store output histograms
  // in the folder "folderTPC"
  //
  TH1::AddDirectory(kFALSE);
  TH1F *h=0;
  TH2D *h2D=0;
  TObjArray *aFolderObj = new TObjArray;
  char name[256];
  char title[256];

  //
  // Cluster histograms
  //
  fTPCClustHisto->GetAxis(2)->SetRange(1,1); // A-side
  h2D = fTPCClustHisto->Projection(1,0);
  h2D->SetName("h_clust_A_side");
  h2D->SetTitle("gclX:gclY - A_side");
  aFolderObj->Add(h2D);

  fTPCClustHisto->GetAxis(2)->SetRange(2,2); // C-side
  h2D = fTPCClustHisto->Projection(1,0);
  h2D->SetName("h_clust_C_side");
  h2D->SetTitle("gclX:gclY - C_side");
  aFolderObj->Add(h2D);

  //
  // event histograms
  //
  for(Int_t i=0; i<6; i++) 
  {
      h = (TH1F*)fTPCEventHisto->Projection(i);
      sprintf(name,"h_tpc_event_%d",i);
      h->SetName(name);
      h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
      h->GetYaxis()->SetTitle("events");
      sprintf(title,"%s",fTPCEventHisto->GetAxis(i)->GetTitle());
      h->SetTitle(title);

      aFolderObj->Add(h);
  }

  // reconstructed vertex status > 0
  fTPCEventHisto->GetAxis(6)->SetRange(2,2);
  for(Int_t i=0; i<6; i++) 
  {
      h = (TH1F*)fTPCEventHisto->Projection(i);
      sprintf(name,"h_tpc_event_recVertex%d",i);
      h->SetName(name);
      h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
      h->GetYaxis()->SetTitle("events");
      sprintf(title,"%s rec. vertex",fTPCEventHisto->GetAxis(i)->GetTitle());
      h->SetTitle(title);

      aFolderObj->Add(h);
  }

  //
  // Track histograms
  //
  for(Int_t i=0; i<9; i++) 
  {
      h = (TH1F*)fTPCTrackHisto->Projection(i);
      sprintf(name,"h_tpc_track_%d",i);
      h->SetName(name);
      h->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
      h->GetYaxis()->SetTitle("tracks");
      sprintf(title,"%s",fTPCTrackHisto->GetAxis(i)->GetTitle());
      h->SetTitle(title);

      if(i==7) h->Scale(1,"width");
      aFolderObj->Add(h);
  }

  //
  for(Int_t i=0; i<8; i++) 
  {
    for(Int_t j=i+1; j<9; j++) 
    {
      h2D = fTPCTrackHisto->Projection(i,j);
      sprintf(name,"h_tpc_track_%d_vs_%d",i,j);
      h2D->SetName(name);
      h2D->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(j)->GetTitle());
      h2D->GetYaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
      sprintf(title,"%s vs %s",fTPCTrackHisto->GetAxis(j)->GetTitle(),fTPCTrackHisto->GetAxis(i)->GetTitle());
      h2D->SetTitle(title);

      if(j==7) h2D->SetBit(TH1::kLogX);
      aFolderObj->Add(h2D);
    }  
  }

  // export objects to analysis folder
  fAnalysisFolder = ExportToFolder(aFolderObj);

  // delete only TObjArray
  if(aFolderObj) delete aFolderObj;
}

//_____________________________________________________________________________
TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array) 
{
  // recreate folder avery time and export objects to new one
  //
  AliPerformanceTPC * comp=this;
  TFolder *folder = comp->GetAnalysisFolder();

  TString name, title;
  TFolder *newFolder = 0;
  Int_t i = 0;
  Int_t size = array->GetSize();

  if(folder) { 
     // get name and title from old folder
     name = folder->GetName();  
     title = folder->GetTitle();  

	 // delete old one
     delete folder;

	 // create new one
     newFolder = CreateFolder(name.Data(),title.Data());
     newFolder->SetOwner();

	 // add objects to folder
     while(i < size) {
	   newFolder->Add(array->At(i));
	   i++;
	 }
  }

return newFolder;
}

//_____________________________________________________________________________
Long64_t AliPerformanceTPC::Merge(TCollection* const list) 
{
  // Merge list of objects (needed by PROOF)

  if (!list)
  return 0;

  if (list->IsEmpty())
  return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj = 0;

  // collection of generated histograms
  Int_t count=0;
  while((obj = iter->Next()) != 0) 
  {
    AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
    if (entry == 0) continue; 

    fTPCClustHisto->Add(entry->fTPCClustHisto);
    fTPCEventHisto->Add(entry->fTPCEventHisto);
    fTPCTrackHisto->Add(entry->fTPCTrackHisto);

    count++;
  }

return count;
}

//_____________________________________________________________________________
TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) { 
// create folder for analysed histograms
//
TFolder *folder = 0;
  folder = new TFolder(name.Data(),title.Data());

  return folder;
}
Commit	Line	Data
7cc34f08	1	//------------------------------------------------------------------------------
	2	// Implementation of AliPerformanceTPC class. It keeps information from
	3	// comparison of reconstructed and MC particle tracks. In addtion,
	4	// it keeps selection cuts used during comparison. The comparison
	5	// information is stored in the ROOT histograms. Analysis of these
	6	// histograms can be done by using Analyse() class function. The result of
	7	// the analysis (histograms/graphs) are stored in the folder which is
	8	// a data member of AliPerformanceTPC.
	9	//
	10	// Author: J.Otwinowski 04/02/2008
	11	//------------------------------------------------------------------------------
	12
	13	/*
	14
	15	// after running comparison task, read the file, and get component
	16	gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
	17	LoadMyLibs();
	18
	19	TFile f("Output.root");
	20	AliPerformanceTPC * compObj = (AliPerformanceTPC*)coutput->FindObject("AliPerformanceTPC");
	21
	22	// analyse comparison data
	23	compObj->Analyse();
	24
	25	// the output histograms/graphs will be stored in the folder "folderTPC"
	26	compObj->GetAnalysisFolder()->ls("*");
	27
	28	// user can save whole comparison object (or only folder with anlysed histograms)
	29	// in the seperate output file (e.g.)
	30	TFile fout("Analysed_TPC.root","recreate");
	31	compObj->Write(); // compObj->GetAnalysisFolder()->Write();
	32	fout.Close();
	33
	34	*/
	35
	36	#include "TCanvas.h"
	37	#include "TH1.h"
	38	#include "TH2.h"
	39	#include "TAxis.h"
	40	#include "TPostScript.h"
	41
	42	#include "AliPerformanceTPC.h"
	43	#include "AliESDEvent.h"
	44	#include "AliESDVertex.h"
	45	#include "AliESDtrack.h"
	46	#include "AliLog.h"
	47	#include "AliMCEvent.h"
	48	#include "AliHeader.h"
	49	#include "AliGenEventHeader.h"
	50	#include "AliStack.h"
	51	#include "AliMCInfoCuts.h"
	52	#include "AliRecInfoCuts.h"
	53	#include "AliTracker.h"
	54	#include "AliTreeDraw.h"
e6a60a90	55	#include "AliTPCTransform.h"
	56	#include "AliTPCseed.h"
	57	#include "AliTPCcalibDB.h"
	58	#include "AliESDfriend.h"
	59	#include "AliESDfriendTrack.h"
	60	#include "AliTPCclusterMI.h"
7cc34f08	61
	62	using namespace std;
	63
	64	ClassImp(AliPerformanceTPC)
	65
	66	//_____________________________________________________________________________
	67	AliPerformanceTPC::AliPerformanceTPC():
	68	AliPerformanceObject("AliPerformanceTPC"),
e6a60a90	69	fTPCClustHisto(0),
a13f4653	70	fTPCEventHisto(0),
a13f4653	71	fTPCTrackHisto(0),
7cc34f08	72
	73	// Cuts
	74	fCutsRC(0),
	75	fCutsMC(0),
	76
	77	// histogram folder
177e2662	78	fAnalysisFolder(0),
	79
	80	// select track vertex
	81	fUseTrackVertex(kFALSE)
7cc34f08	82	{
	83	Init();
	84	}
	85
	86	//_____________________________________________________________________________
	87	AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
	88	AliPerformanceObject(name,title),
e6a60a90	89	fTPCClustHisto(0),
a13f4653	90	fTPCEventHisto(0),
a13f4653	91	fTPCTrackHisto(0),
7cc34f08	92
	93	// Cuts
	94	fCutsRC(0),
	95	fCutsMC(0),
	96
	97	// histogram folder
177e2662	98	fAnalysisFolder(0),
	99
	100	// select track vertex
	101	fUseTrackVertex(kFALSE)
7cc34f08	102	{
	103	// named constructor
	104	//
	105	SetAnalysisMode(analysisMode);
	106	SetHptGenerator(hptGenerator);
	107
	108	Init();
	109	}
	110
	111	//_____________________________________________________________________________
	112	AliPerformanceTPC::~AliPerformanceTPC()
	113	{
	114	// destructor
	115
e6a60a90	116	if(fTPCClustHisto) delete fTPCClustHisto; fTPCClustHisto=0;
a13f4653	117	if(fTPCEventHisto) delete fTPCEventHisto; fTPCEventHisto=0;
a13f4653	118	if(fTPCTrackHisto) delete fTPCTrackHisto; fTPCTrackHisto=0;
7cc34f08	119	if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
	120	}
	121
	122	//_____________________________________________________________________________
	123	void AliPerformanceTPC::Init(){
	124	//
	125	// histogram bining
	126	//
	127
	128	// set pt bins
	129	Int_t nPtBins = 50;
	130	Double_t ptMin = 1.e-2, ptMax = 10.;
	131
	132	Double_t *binsPt = 0;
	133	if (IsHptGenerator()) {
	134	nPtBins = 100; ptMax = 100.;
	135	binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
	136	} else {
	137	binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
	138	}
	139
	140	/*
	141	Int_t nPtBins = 31;
	142	Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
	143	Double_t ptMin = 0., ptMax = 10.;
	144
	145	if(IsHptGenerator() == kTRUE) {
	146	nPtBins = 100;
	147	ptMin = 0.; ptMax = 100.;
	148	}
	149	*/
e6a60a90	150	//
	151	//gclX:gclY:TPCSide
	152	Int_t binsTPCClustHisto[3]= { 500, 500, 2 };
	153	Double_t minTPCClustHisto[3]={-250., -250., 0.};
	154	Double_t maxTPCClustHisto[3]={ 250., 250., 2.};
	155
	156	fTPCClustHisto = new THnSparseF("fTPCClustHisto","gclX:gclY:TPCSide",3,binsTPCClustHisto,minTPCClustHisto,maxTPCClustHisto);
	157	fTPCClustHisto->GetAxis(0)->SetTitle("gclX (cm)");
	158	fTPCClustHisto->GetAxis(1)->SetTitle("gclY (cm)");
	159	fTPCClustHisto->GetAxis(2)->SetTitle("TPCSide");
	160	fTPCClustHisto->Sumw2();
	161
7cc34f08	162
886bf4d3	163	// Xv:Yv:Zv:mult:multP:multN:vertStatus
	164	Int_t binsTPCEventHisto[7]= {100, 100, 100, 151, 151, 151, 2 };
	165	Double_t minTPCEventHisto[7]={-10., -10., -30., -0.5, -0.5, -0.5, 0. };
	166	Double_t maxTPCEventHisto[7]={ 10., 10., 30., 150.5, 150.5, 150.5, 2. };
a13f4653	167
886bf4d3	168	fTPCEventHisto = new THnSparseF("fTPCEventHisto","Xv:Yv:Zv:mult:multP:multN:vertStatus",7,binsTPCEventHisto,minTPCEventHisto,maxTPCEventHisto);
a13f4653	169	fTPCEventHisto->GetAxis(0)->SetTitle("Xv (cm)");
	170	fTPCEventHisto->GetAxis(1)->SetTitle("Yv (cm)");
	171	fTPCEventHisto->GetAxis(2)->SetTitle("Zv (cm)");
	172	fTPCEventHisto->GetAxis(3)->SetTitle("mult");
	173	fTPCEventHisto->GetAxis(4)->SetTitle("multP");
	174	fTPCEventHisto->GetAxis(5)->SetTitle("multN");
886bf4d3	175	fTPCEventHisto->GetAxis(6)->SetTitle("vertStatus");
a13f4653	176	fTPCEventHisto->Sumw2();
	177
	178
	179	// nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:DCAr:DCAz:eta:phi:pt:charge
e6a60a90	180	Int_t binsTPCTrackHisto[9]= { 160, 50, 60, 100, 100, 30, 144, nPtBins, 3 };
a13f4653	181	Double_t minTPCTrackHisto[9]={ 0., 0., 0., -10, -10., -1.5, 0., ptMin, -1.5};
	182	Double_t maxTPCTrackHisto[9]={ 160., 10., 1.2, 10, 10., 1.5, 2.*TMath::Pi(), ptMax, 1.5};
	183
	184	fTPCTrackHisto = new THnSparseF("fTPCTrackHisto","nClust:chi2PerClust:nClust/nFindableClust:DCAr:DCAz:eta:phi:pt:charge",9,binsTPCTrackHisto,minTPCTrackHisto,maxTPCTrackHisto);
	185	fTPCTrackHisto->SetBinEdges(7,binsPt);
	186
	187	fTPCTrackHisto->GetAxis(0)->SetTitle("nClust");
	188	fTPCTrackHisto->GetAxis(1)->SetTitle("chi2PerClust");
	189	fTPCTrackHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
	190	fTPCTrackHisto->GetAxis(3)->SetTitle("DCAr (cm)");
	191	fTPCTrackHisto->GetAxis(4)->SetTitle("DCAz (cm)");
	192	fTPCTrackHisto->GetAxis(5)->SetTitle("#eta");
	193	fTPCTrackHisto->GetAxis(6)->SetTitle("#phi (rad)");
	194	fTPCTrackHisto->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
	195	fTPCTrackHisto->GetAxis(8)->SetTitle("charge");
	196	fTPCTrackHisto->Sumw2();
7cc34f08	197
	198	// Init cuts
	199	if(!fCutsMC)
	200	AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
	201	if(!fCutsRC)
	202	AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
	203
	204	// init folder
	205	fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
	206	}
	207
	208	//_____________________________________________________________________________
	209	void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
	210	{
	211	if(!esdTrack) return;
	212
	213	// Fill TPC only resolution comparison information
	214	const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
	215	if(!track) return;
	216
	217	Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
	218	esdTrack->GetImpactParametersTPC(dca,cov);
	219
a13f4653	220	Float_t q = esdTrack->Charge();
7cc34f08	221	Float_t pt = track->Pt();
	222	Float_t eta = track->Eta();
	223	Float_t phi = track->Phi();
	224	Int_t nClust = esdTrack->GetTPCclusters(0);
	225	Int_t nFindableClust = esdTrack->GetTPCNclsF();
	226
	227	Float_t chi2PerCluster = 0.;
	228	if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
	229
	230	Float_t clustPerFindClust = 0.;
	231	if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
	232
a13f4653	233	//
	234	// select primaries
	235	//
	236	Double_t dcaToVertex = -1;
	237	if( fCutsRC->GetDCAToVertex2D() )
	238	{
177e2662	239	dcaToVertex = TMath::Sqrt(dca[0]dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
7cc34f08	240	}
a13f4653	241	if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
	242	if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
	243	if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
	244
	245	Double_t vTPCTrackHisto[9] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,q};
	246	fTPCTrackHisto->Fill(vTPCTrackHisto);
7cc34f08	247
	248	//
	249	// Fill rec vs MC information
	250	//
	251	if(!stack) return;
	252
	253	}
	254
	255	//_____________________________________________________________________________
	256	void AliPerformanceTPC::ProcessTPCITS(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
	257	{
	258	// Fill comparison information (TPC+ITS)
	259	AliDebug(AliLog::kWarning, "Warning: Not implemented");
	260	}
	261
	262	//_____________________________________________________________________________
	263	void AliPerformanceTPC::ProcessConstrained(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
	264	{
	265	// Fill comparison information (constarained parameters)
	266	AliDebug(AliLog::kWarning, "Warning: Not implemented");
	267	}
	268
	269	//_____________________________________________________________________________
	270	void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent const esdEvent, AliESDfriend const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
	271	{
	272	// Process comparison information
	273	//
	274	if(!esdEvent)
	275	{
	276	Error("Exec","esdEvent not available");
	277	return;
	278	}
	279	AliHeader* header = 0;
	280	AliGenEventHeader* genHeader = 0;
	281	AliStack* stack = 0;
	282	TArrayF vtxMC(3);
	283
	284	if(bUseMC)
	285	{
	286	if(!mcEvent) {
	287	Error("Exec","mcEvent not available");
	288	return;
	289	}
	290	// get MC event header
	291	header = mcEvent->Header();
	292	if (!header) {
	293	Error("Exec","Header not available");
	294	return;
	295	}
	296	// MC particle stack
	297	stack = mcEvent->Stack();
	298	if (!stack) {
	299	Error("Exec","Stack not available");
	300	return;
	301	}
	302	// get MC vertex
	303	genHeader = header->GenEventHeader();
	304	if (!genHeader) {
	305	Error("Exec","Could not retrieve genHeader from Header");
	306	return;
	307	}
	308	genHeader->PrimaryVertex(vtxMC);
	309	}
	310
e6a60a90	311	// trigger
046b4280	312	if(!bUseMC &&GetTriggerClass()) {
4bb5f8a9	313	Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
	314	if(!isEventTriggered) return;
	315	}
e6a60a90	316
e6a60a90	317	// get TPC event vertex
177e2662	318	const AliESDVertex *vtxESD = NULL;
	319	if(fUseTrackVertex) {
	320	vtxESD = esdEvent->GetPrimaryVertexTracks();
	321	} else {
	322	vtxESD = esdEvent->GetPrimaryVertexTPC();
	323	}
9ca899e5	324	if(!vtxESD) return;
e6a60a90	325
e6a60a90	326	// events with rec. vertex
a13f4653	327	Int_t mult=0; Int_t multP=0; Int_t multN=0;
e6a60a90	328	if(vtxESD->GetStatus() >0)
	329	{
	330	// Process ESD events
7cc34f08	331	for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
7cc34f08	332	{
a13f4653	333	AliESDtrack *track = esdEvent->GetTrack(iTrack);
7cc34f08	334	if(!track) continue;
7cc34f08	335
4e59b45c	336	if(bUseESDfriend && esdFriend)
	337	{
	338	AliESDfriendTrack *friendTrack=esdFriend->GetTrack(iTrack);
177e2662	339	if(friendTrack)
	340	{
	341	TObject *calibObject=0;
	342	AliTPCseed *seed=0;
	343	if (!friendTrack) continue;
	344	for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j) {
	345	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) {
	346	break;
	347	}
e6a60a90	348	}
	349
	350	//AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
	351	for (Int_t irow=0;irow<159;irow++) {
	352	if(!seed) continue;
	353
	354	AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
	355	if (!cluster) continue;
	356
	357	Float_t gclf[3];
	358	cluster->GetGlobalXYZ(gclf);
	359
	360	//Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
	361	//Int_t i[1]={cluster->GetDetector()};
	362	//transform->Transform(x,i,0,1);
	363	//printf("gx %f gy %f gz %f \n", cluster->GetX(), cluster->GetY(),cluster->GetZ());
	364	//printf("gclf[0] %f gclf[1] %f gclf[2] %f \n", gclf[0], gclf[1], gclf[2]);
	365
	366	Int_t TPCside;
	367	if(gclf[2]>0.) TPCside=0; // A side
	368	else TPCside=1;
	369
	370	Double_t vTPCClust[3] = { gclf[0], gclf[1], TPCside };
	371	fTPCClustHisto->Fill(vTPCClust);
	372	}
177e2662	373	}
e6a60a90	374	}
177e2662	375
7cc34f08	376	if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
	377	else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
	378	else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
	379	else {
	380	printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
	381	return;
	382	}
a13f4653	383
	384	// TPC only
	385	AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
	386	if(!tpcTrack) continue;
	387
	388	// track selection
	389	if( fCutsRC->AcceptTrack(tpcTrack) ) {
	390	mult++;
	391	if(tpcTrack->Charge()>0.) multP++;
	392	if(tpcTrack->Charge()<0.) multN++;
	393	}
	394
	395	if(tpcTrack) delete tpcTrack;
7cc34f08	396	}
e6a60a90	397	}
a13f4653	398	//
e6a60a90	399
886bf4d3	400	Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
a13f4653	401	fTPCEventHisto->Fill(vTPCEvent);
7cc34f08	402	}
	403
	404	//_____________________________________________________________________________
	405	void AliPerformanceTPC::Analyse() {
	406	//
	407	// Analyse comparison information and store output histograms
	408	// in the folder "folderTPC"
	409	//
	410	TH1::AddDirectory(kFALSE);
a13f4653	411	TH1F *h=0;
e6a60a90	412	TH2D *h2D=0;
7cc34f08	413	TObjArray *aFolderObj = new TObjArray;
7cc34f08	414	char name[256];
7cc34f08	415	char title[256];
a13f4653	416
e6a60a90	417	//
	418	// Cluster histograms
	419	//
	420	fTPCClustHisto->GetAxis(2)->SetRange(1,1); // A-side
60446b2a	421	h2D = fTPCClustHisto->Projection(1,0);
e6a60a90	422	h2D->SetName("h_clust_A_side");
	423	h2D->SetTitle("gclX:gclY - A_side");
	424	aFolderObj->Add(h2D);
	425
	426	fTPCClustHisto->GetAxis(2)->SetRange(2,2); // C-side
60446b2a	427	h2D = fTPCClustHisto->Projection(1,0);
e6a60a90	428	h2D->SetName("h_clust_C_side");
	429	h2D->SetTitle("gclX:gclY - C_side");
	430	aFolderObj->Add(h2D);
	431
a13f4653	432	//
	433	// event histograms
	434	//
	435	for(Int_t i=0; i<6; i++)
	436	{
	437	h = (TH1F*)fTPCEventHisto->Projection(i);
	438	sprintf(name,"h_tpc_event_%d",i);
	439	h->SetName(name);
	440	h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
	441	h->GetYaxis()->SetTitle("events");
	442	sprintf(title,"%s",fTPCEventHisto->GetAxis(i)->GetTitle());
	443	h->SetTitle(title);
	444
	445	aFolderObj->Add(h);
	446	}
	447
886bf4d3	448	// reconstructed vertex status > 0
	449	fTPCEventHisto->GetAxis(6)->SetRange(2,2);
	450	for(Int_t i=0; i<6; i++)
	451	{
	452	h = (TH1F*)fTPCEventHisto->Projection(i);
	453	sprintf(name,"h_tpc_event_recVertex%d",i);
	454	h->SetName(name);
	455	h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
	456	h->GetYaxis()->SetTitle("events");
	457	sprintf(title,"%s rec. vertex",fTPCEventHisto->GetAxis(i)->GetTitle());
	458	h->SetTitle(title);
	459
	460	aFolderObj->Add(h);
	461	}
	462
a13f4653	463	//
	464	// Track histograms
	465	//
	466	for(Int_t i=0; i<9; i++)
	467	{
	468	h = (TH1F*)fTPCTrackHisto->Projection(i);
	469	sprintf(name,"h_tpc_track_%d",i);
	470	h->SetName(name);
	471	h->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
	472	h->GetYaxis()->SetTitle("tracks");
	473	sprintf(title,"%s",fTPCTrackHisto->GetAxis(i)->GetTitle());
	474	h->SetTitle(title);
	475
	476	if(i==7) h->Scale(1,"width");
	477	aFolderObj->Add(h);
	478	}
	479
	480	//
	481	for(Int_t i=0; i<8; i++)
7cc34f08	482	{
a13f4653	483	for(Int_t j=i+1; j<9; j++)
7cc34f08	484	{
e6a60a90	485	h2D = fTPCTrackHisto->Projection(i,j);
a13f4653	486	sprintf(name,"h_tpc_track_%d_vs_%d",i,j);
7cc34f08	487	h2D->SetName(name);
a13f4653	488	h2D->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(j)->GetTitle());
	489	h2D->GetYaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
	490	sprintf(title,"%s vs %s",fTPCTrackHisto->GetAxis(j)->GetTitle(),fTPCTrackHisto->GetAxis(i)->GetTitle());
7cc34f08	491	h2D->SetTitle(title);
7cc34f08	492
a13f4653	493	if(j==7) h2D->SetBit(TH1::kLogX);
7cc34f08	494	aFolderObj->Add(h2D);
	495	}
	496	}
	497
	498	// export objects to analysis folder
	499	fAnalysisFolder = ExportToFolder(aFolderObj);
	500
	501	// delete only TObjArray
	502	if(aFolderObj) delete aFolderObj;
	503	}
	504
	505	//_____________________________________________________________________________
	506	TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array)
	507	{
	508	// recreate folder avery time and export objects to new one
	509	//
	510	AliPerformanceTPC * comp=this;
	511	TFolder *folder = comp->GetAnalysisFolder();
	512
	513	TString name, title;
	514	TFolder *newFolder = 0;
	515	Int_t i = 0;
	516	Int_t size = array->GetSize();
	517
	518	if(folder) {
	519	// get name and title from old folder
	520	name = folder->GetName();
	521	title = folder->GetTitle();
	522
	523	// delete old one
	524	delete folder;
	525
	526	// create new one
	527	newFolder = CreateFolder(name.Data(),title.Data());
	528	newFolder->SetOwner();
	529
	530	// add objects to folder
	531	while(i < size) {
	532	newFolder->Add(array->At(i));
	533	i++;
	534	}
	535	}
	536
	537	return newFolder;
	538	}
	539
	540	//_____________________________________________________________________________
	541	Long64_t AliPerformanceTPC::Merge(TCollection* const list)
	542	{
	543	// Merge list of objects (needed by PROOF)
	544
	545	if (!list)
	546	return 0;
	547
	548	if (list->IsEmpty())
	549	return 1;
	550
	551	TIterator* iter = list->MakeIterator();
	552	TObject* obj = 0;
	553
	554	// collection of generated histograms
	555	Int_t count=0;
	556	while((obj = iter->Next()) != 0)
	557	{
a13f4653	558	AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
a13f4653	559	if (entry == 0) continue;
7cc34f08	560
e6a60a90	561	fTPCClustHisto->Add(entry->fTPCClustHisto);
a13f4653	562	fTPCEventHisto->Add(entry->fTPCEventHisto);
a13f4653	563	fTPCTrackHisto->Add(entry->fTPCTrackHisto);
7cc34f08	564
a13f4653	565	count++;
7cc34f08	566	}
	567
	568	return count;
	569	}
	570
	571	//_____________________________________________________________________________
	572	TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) {
	573	// create folder for analysed histograms
	574	//
	575	TFolder *folder = 0;
	576	folder = new TFolder(name.Data(),title.Data());
	577
	578	return folder;
	579	}