[u/mrichter/AliRoot.git] / MFT / AliMFTClusterQA.cxx

#include "TObject.h"
#include "AliRunLoader.h"
#include "AliRun.h"
#include "AliLoader.h"
#include "AliMFT.h"
#include "TClonesArray.h"
#include "AliMFTCluster.h"
#include "AliMFTSegmentation.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "AliLog.h"
#include "TString.h"

#include "AliMFTClusterQA.h"

//====================================================================================================================================================
//
// Class for the analysis of the MFT clusters (a.k.a. rec points). Few QA histograms are created
//
// Contact author: antonio.uras@cern.ch
//
//====================================================================================================================================================

ClassImp(AliMFTClusterQA)

//====================================================================================================================================================

AliMFTClusterQA::AliMFTClusterQA():
  TObject(),
  fMFTLoader(0),
  fRunLoader(0),
  fMFT(0),
  fNPlanes(0),
  fNEvents(0),
  fEv(0),
  fFileOut(0),
  fReadDir(0),
  fOutDir(0)
{
  
  // default constructor

  for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
    fHistNClustersPerEvent[iPlane]     = 0;
    fHistNPixelsPerCluster[iPlane]     = 0;
    fHistClusterSizeX[iPlane]          = 0; 
    fHistClusterSizeY[iPlane]          = 0;
    fHistClusterRadialPosition[iPlane] = 0;
    fClusterScatterPlotXY[iPlane]      = 0;
  }

}

//====================================================================================================================================================

void AliMFTClusterQA::Init(Char_t *readDir, Char_t *outDir, Int_t nEventsToAnalyze) {

  fReadDir = readDir;
  fOutDir  = outDir;

  fRunLoader = AliRunLoader::Open(Form("%s/galice.root", fReadDir.Data()));
  gAlice = fRunLoader->GetAliRun();
  if (!gAlice) fRunLoader->LoadgAlice();
  fNEvents = fRunLoader->GetNumberOfEvents();
  if (nEventsToAnalyze>0) fNEvents = TMath::Min(fNEvents, nEventsToAnalyze);

  fMFT = (AliMFT*) gAlice->GetDetector("MFT"); 
  fNPlanes = fMFT->GetSegmentation()->GetNPlanes();

  BookHistos();

  fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
  fMFTLoader -> LoadRecPoints("READ");

}

//====================================================================================================================================================

Bool_t AliMFTClusterQA::LoadNextEvent() {

  if (fEv>=fNEvents) return kFALSE;
  AliDebug(1, Form("event %5d",fEv));
  
  fRunLoader->GetEvent(fEv);
  fEv++;

  if (!fMFTLoader->TreeR()->GetEvent()) return kTRUE;

  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
    Int_t nClusters = fMFT->GetRecPointsList(iPlane)->GetEntries();
    fHistNClustersPerEvent[iPlane] -> Fill(nClusters);
    fClusterScatterPlotXY[iPlane]  -> Fill(0., 0.);    // "scaler" bin
    AliDebug(1,Form("nClusters = %5d", nClusters));
    for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
      AliMFTCluster *cluster = (AliMFTCluster*) (fMFT->GetRecPointsList(iPlane))->At(iCluster);
      fHistNPixelsPerCluster[iPlane] -> Fill(cluster->GetSize());
      fHistClusterSizeX[iPlane]      -> Fill(cluster->GetErrX()*1.e4);   // converted in microns
      fHistClusterSizeY[iPlane]      -> Fill(cluster->GetErrY()*1.e4);   // converted in microns
      fHistClusterRadialPosition[iPlane] -> Fill(TMath::Sqrt(cluster->GetX()*cluster->GetX()+cluster->GetY()*cluster->GetY()));
      fClusterScatterPlotXY[iPlane]      -> Fill(cluster->GetX(), cluster->GetY());
    }
  }

  return kTRUE;

}  

//====================================================================================================================================================

void AliMFTClusterQA::BookHistos() {

  fFileOut = new TFile(Form("%s/MFT.RecPoints.QA.root",fOutDir.Data()), "recreate");

  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {

    fHistNClustersPerEvent[iPlane] = new TH1D(Form("fHistNClustersPerEvent_Pl%02d",iPlane), 
					      Form("Number of clusters per event in Plane%02d",iPlane),
					      25000, -0.5, 24999.5);

    fHistNPixelsPerCluster[iPlane] = new TH1D(Form("fHistNPixelsPerCluster_Pl%02d",iPlane), 
					      Form("Number of pixels per cluster in Plane%02d",iPlane),
					      15, -0.5, 14.5);

    fHistClusterSizeX[iPlane]      = new TH1D(Form("fHistClusterSizeX_Pl%02d",iPlane), 
					      Form("#Deltax for clusters in Plane%02d",iPlane), 
					      100, 0., 100.);
    
    fHistClusterSizeY[iPlane]      = new TH1D(Form("fHistClusterSizeY_Pl%02d",iPlane), 
					      Form("#Deltay for clusters in Plane%02d",iPlane), 
					      100, 0., 100.);
    
    fHistNClustersPerEvent[iPlane] -> SetXTitle("N_{clusters} per Event");
    fHistNPixelsPerCluster[iPlane] -> SetXTitle("N_{pixels} per Cluster");
    fHistClusterSizeX[iPlane]      -> SetXTitle("#Deltax  [#mum]");
    fHistClusterSizeY[iPlane]      -> SetXTitle("#Deltay  [#mum]");

    fHistNClustersPerEvent[iPlane] -> Sumw2();
    fHistNPixelsPerCluster[iPlane] -> Sumw2();
    fHistClusterSizeX[iPlane]      -> Sumw2();
    fHistClusterSizeY[iPlane]      -> Sumw2();

    //------------------------------------------------------------

    Int_t rMax = Int_t(10.*(fMFT->GetSegmentation()->GetPlane(iPlane)->GetRMaxSupport()));  // expressed in mm

    fHistClusterRadialPosition[iPlane] = new TH1D(Form("fHistClusterRadialPosition_Pl%02d",iPlane),
						  Form("Cluster radial position (Plane%02d)",iPlane),
						  rMax, 0, Double_t(rMax)/10.);

    fClusterScatterPlotXY[iPlane] = new TH2D(Form("fClusterScatterPlotXY_Pl%02d",iPlane),
					     Form("Cluster scatter plot (Plane%02d)",iPlane),
					     2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10., 2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10.);
    
    fHistClusterRadialPosition[iPlane] -> SetXTitle("R  [cm]");
    fClusterScatterPlotXY[iPlane]      -> SetXTitle("X  [cm]");
    fClusterScatterPlotXY[iPlane]      -> SetYTitle("Y  [cm]");

    fHistClusterRadialPosition[iPlane] -> Sumw2();
    fClusterScatterPlotXY[iPlane]      -> Sumw2();
    
  }
  
}

//====================================================================================================================================================

void AliMFTClusterQA::Terminate() {

  AliInfo("Writing QA histos...");

  // ---- equalize radial clusters distribution ----------------------

  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
    for (Int_t iBin=1; iBin<=fHistClusterRadialPosition[iPlane]->GetNbinsX(); iBin++) {
      Double_t rMin = fHistClusterRadialPosition[iPlane]->GetBinLowEdge(iBin);        // in cm
      Double_t rMax = fHistClusterRadialPosition[iPlane]->GetBinWidth(iBin) + rMin;   // in cm
      Double_t area = 100.*TMath::Pi()*(rMax*rMax - rMin*rMin);                       // in mm^2
      Double_t density = fHistClusterRadialPosition[iPlane]->GetBinContent(iBin)/area;
      fHistClusterRadialPosition[iPlane]->SetBinContent(iBin, density);
      fHistClusterRadialPosition[iPlane]->SetBinError(iBin, fHistClusterRadialPosition[iPlane]->GetBinError(iBin)/area);
    }
    fHistClusterRadialPosition[iPlane] -> SetBinContent(1, fEv);      // "scaler" bin
  }

  // -----------------------------------------------------------------

  fFileOut->cd();

  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
    fHistNClustersPerEvent[iPlane]     -> Write();
    fHistNPixelsPerCluster[iPlane]     -> Write();
    fHistClusterSizeX[iPlane]          -> Write();
    fHistClusterSizeY[iPlane]          -> Write();
    fHistClusterRadialPosition[iPlane] -> Write();
    fClusterScatterPlotXY[iPlane]      -> Write();
  }

  fFileOut -> Close();

}

//====================================================================================================================================================
Commit	Line	Data
820b4d9e	1	#include "TObject.h"
	2	#include "AliRunLoader.h"
	3	#include "AliRun.h"
	4	#include "AliLoader.h"
	5	#include "AliMFT.h"
	6	#include "TClonesArray.h"
	7	#include "AliMFTCluster.h"
	8	#include "AliMFTSegmentation.h"
	9	#include "TFile.h"
	10	#include "TH1D.h"
d4643a10	11	#include "TH2D.h"
820b4d9e	12	#include "AliLog.h"
	13	#include "TString.h"
	14
	15	#include "AliMFTClusterQA.h"
	16
	17	//====================================================================================================================================================
	18	//
	19	// Class for the analysis of the MFT clusters (a.k.a. rec points). Few QA histograms are created
	20	//
	21	// Contact author: antonio.uras@cern.ch
	22	//
	23	//====================================================================================================================================================
	24
	25	ClassImp(AliMFTClusterQA)
	26
	27	//====================================================================================================================================================
	28
	29	AliMFTClusterQA::AliMFTClusterQA():
	30	TObject(),
	31	fMFTLoader(0),
	32	fRunLoader(0),
	33	fMFT(0),
	34	fNPlanes(0),
	35	fNEvents(0),
	36	fEv(0),
	37	fFileOut(0),
	38	fReadDir(0),
	39	fOutDir(0)
	40	{
	41
	42	// default constructor
	43
78461b61	44	for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
bcaf50eb	45	fHistNClustersPerEvent[iPlane] = 0;
	46	fHistNPixelsPerCluster[iPlane] = 0;
	47	fHistClusterSizeX[iPlane] = 0;
	48	fHistClusterSizeY[iPlane] = 0;
	49	fHistClusterRadialPosition[iPlane] = 0;
	50	fClusterScatterPlotXY[iPlane] = 0;
820b4d9e	51	}
	52
	53	}
	54
	55	//====================================================================================================================================================
	56
	57	void AliMFTClusterQA::Init(Char_t readDir, Char_t outDir, Int_t nEventsToAnalyze) {
	58
	59	fReadDir = readDir;
	60	fOutDir = outDir;
	61
	62	fRunLoader = AliRunLoader::Open(Form("%s/galice.root", fReadDir.Data()));
	63	gAlice = fRunLoader->GetAliRun();
	64	if (!gAlice) fRunLoader->LoadgAlice();
	65	fNEvents = fRunLoader->GetNumberOfEvents();
	66	if (nEventsToAnalyze>0) fNEvents = TMath::Min(fNEvents, nEventsToAnalyze);
	67
	68	fMFT = (AliMFT*) gAlice->GetDetector("MFT");
	69	fNPlanes = fMFT->GetSegmentation()->GetNPlanes();
	70
	71	BookHistos();
	72
	73	fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
	74	fMFTLoader -> LoadRecPoints("READ");
	75
	76	}
	77
	78	//====================================================================================================================================================
	79
	80	Bool_t AliMFTClusterQA::LoadNextEvent() {
	81
	82	if (fEv>=fNEvents) return kFALSE;
	83	AliDebug(1, Form("event %5d",fEv));
	84
	85	fRunLoader->GetEvent(fEv);
	86	fEv++;
	87
	88	if (!fMFTLoader->TreeR()->GetEvent()) return kTRUE;
	89
	90	for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
	91	Int_t nClusters = fMFT->GetRecPointsList(iPlane)->GetEntries();
	92	fHistNClustersPerEvent[iPlane] -> Fill(nClusters);
d4643a10	93	fClusterScatterPlotXY[iPlane] -> Fill(0., 0.); // "scaler" bin
820b4d9e	94	AliDebug(1,Form("nClusters = %5d", nClusters));
	95	for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
	96	AliMFTCluster cluster = (AliMFTCluster) (fMFT->GetRecPointsList(iPlane))->At(iCluster);
	97	fHistNPixelsPerCluster[iPlane] -> Fill(cluster->GetSize());
	98	fHistClusterSizeX[iPlane] -> Fill(cluster->GetErrX()*1.e4); // converted in microns
	99	fHistClusterSizeY[iPlane] -> Fill(cluster->GetErrY()*1.e4); // converted in microns
bcaf50eb	100	fHistClusterRadialPosition[iPlane] -> Fill(TMath::Sqrt(cluster->GetX()cluster->GetX()+cluster->GetY()cluster->GetY()));
bcaf50eb	101	fClusterScatterPlotXY[iPlane] -> Fill(cluster->GetX(), cluster->GetY());
820b4d9e	102	}
	103	}
	104
	105	return kTRUE;
	106
	107	}
	108
	109	//====================================================================================================================================================
	110
	111	void AliMFTClusterQA::BookHistos() {
	112
	113	fFileOut = new TFile(Form("%s/MFT.RecPoints.QA.root",fOutDir.Data()), "recreate");
	114
	115	for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
	116
	117	fHistNClustersPerEvent[iPlane] = new TH1D(Form("fHistNClustersPerEvent_Pl%02d",iPlane),
	118	Form("Number of clusters per event in Plane%02d",iPlane),
b03d16a3	119	25000, -0.5, 24999.5);
820b4d9e	120
	121	fHistNPixelsPerCluster[iPlane] = new TH1D(Form("fHistNPixelsPerCluster_Pl%02d",iPlane),
	122	Form("Number of pixels per cluster in Plane%02d",iPlane),
	123	15, -0.5, 14.5);
	124
	125	fHistClusterSizeX[iPlane] = new TH1D(Form("fHistClusterSizeX_Pl%02d",iPlane),
	126	Form("#Deltax for clusters in Plane%02d",iPlane),
	127	100, 0., 100.);
	128
	129	fHistClusterSizeY[iPlane] = new TH1D(Form("fHistClusterSizeY_Pl%02d",iPlane),
	130	Form("#Deltay for clusters in Plane%02d",iPlane),
	131	100, 0., 100.);
	132
	133	fHistNClustersPerEvent[iPlane] -> SetXTitle("N_{clusters} per Event");
	134	fHistNPixelsPerCluster[iPlane] -> SetXTitle("N_{pixels} per Cluster");
	135	fHistClusterSizeX[iPlane] -> SetXTitle("#Deltax [#mum]");
	136	fHistClusterSizeY[iPlane] -> SetXTitle("#Deltay [#mum]");
	137
	138	fHistNClustersPerEvent[iPlane] -> Sumw2();
	139	fHistNPixelsPerCluster[iPlane] -> Sumw2();
	140	fHistClusterSizeX[iPlane] -> Sumw2();
	141	fHistClusterSizeY[iPlane] -> Sumw2();
d4643a10	142
	143	//------------------------------------------------------------
	144
bcaf50eb	145	Int_t rMax = Int_t(10.*(fMFT->GetSegmentation()->GetPlane(iPlane)->GetRMaxSupport())); // expressed in mm
	146
	147	fHistClusterRadialPosition[iPlane] = new TH1D(Form("fHistClusterRadialPosition_Pl%02d",iPlane),
	148	Form("Cluster radial position (Plane%02d)",iPlane),
	149	rMax, 0, Double_t(rMax)/10.);
	150
d4643a10	151	fClusterScatterPlotXY[iPlane] = new TH2D(Form("fClusterScatterPlotXY_Pl%02d",iPlane),
	152	Form("Cluster scatter plot (Plane%02d)",iPlane),
	153	2rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10., 2rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10.);
	154
bcaf50eb	155	fHistClusterRadialPosition[iPlane] -> SetXTitle("R [cm]");
	156	fClusterScatterPlotXY[iPlane] -> SetXTitle("X [cm]");
	157	fClusterScatterPlotXY[iPlane] -> SetYTitle("Y [cm]");
	158
	159	fHistClusterRadialPosition[iPlane] -> Sumw2();
	160	fClusterScatterPlotXY[iPlane] -> Sumw2();
820b4d9e	161
	162	}
	163
	164	}
	165
	166	//====================================================================================================================================================
	167
	168	void AliMFTClusterQA::Terminate() {
	169
	170	AliInfo("Writing QA histos...");
	171
bcaf50eb	172	// ---- equalize radial clusters distribution ----------------------
	173
	174	for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
	175	for (Int_t iBin=1; iBin<=fHistClusterRadialPosition[iPlane]->GetNbinsX(); iBin++) {
	176	Double_t rMin = fHistClusterRadialPosition[iPlane]->GetBinLowEdge(iBin); // in cm
	177	Double_t rMax = fHistClusterRadialPosition[iPlane]->GetBinWidth(iBin) + rMin; // in cm
	178	Double_t area = 100.TMath::Pi()(rMaxrMax - rMinrMin); // in mm^2
	179	Double_t density = fHistClusterRadialPosition[iPlane]->GetBinContent(iBin)/area;
	180	fHistClusterRadialPosition[iPlane]->SetBinContent(iBin, density);
	181	fHistClusterRadialPosition[iPlane]->SetBinError(iBin, fHistClusterRadialPosition[iPlane]->GetBinError(iBin)/area);
	182	}
	183	fHistClusterRadialPosition[iPlane] -> SetBinContent(1, fEv); // "scaler" bin
	184	}
	185
	186	// -----------------------------------------------------------------
	187
820b4d9e	188	fFileOut->cd();
	189
	190	for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
bcaf50eb	191	fHistNClustersPerEvent[iPlane] -> Write();
	192	fHistNPixelsPerCluster[iPlane] -> Write();
	193	fHistClusterSizeX[iPlane] -> Write();
	194	fHistClusterSizeY[iPlane] -> Write();
	195	fHistClusterRadialPosition[iPlane] -> Write();
	196	fClusterScatterPlotXY[iPlane] -> Write();
820b4d9e	197	}
	198
	199	fFileOut -> Close();
	200
	201	}
	202
	203	//====================================================================================================================================================