[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskVtXY.cxx

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

//--------------------------------------------------------------------------
//    Implementation of the beam interaction spot location and size estimate
//    via VertexerTracksNoConstraint resolution extrapolation
//    to the infinit multiplicity alias zero resolution
//
// Origin: AliAnalysisTaskVtXY
//         A. Jacholkowski, Catania 
//         adam.jacholkowski@cern.ch
//-----------------------------------------------------------------
#include "TChain.h"
#include "TTree.h"
#include "TStyle.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TCanvas.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"

#include "AliAnalysisTaskVtXY.h"

#include "AliESDVertex.h"
#include "AliVertexerTracks.h"
ClassImp(AliAnalysisTaskVtXY)

//________________________________________________________________________
AliAnalysisTaskVtXY::AliAnalysisTaskVtXY(const char *name) 
  : AliAnalysisTask(name, ""),
  fESD(0), 
  fList(0),
  fHistVtx(0),
  fHistVty(0)
{
  // Constructor

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  DefineOutput(0, TList::Class());
  // Output slot #0 writes into a TList container
}

//________________________________________________________________________
void AliAnalysisTaskVtXY::ConnectInputData(Option_t *) 
{
  // Connect ESD or AOD here
  // Called once

  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } else {
    // Disable all branches and enable only the needed ones
    // The next two lines are different when data produced as AliESDEvent is read
    //tree->SetBranchStatus("*", kFALSE);
    //tree->SetBranchStatus("Tracks.*", kTRUE);

    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());

    if (!esdH) {
      Printf("ERROR: Could not get ESDInputHandler");
    } else
      fESD = esdH->GetEvent();
  }
}

//________________________________________________________________________
void AliAnalysisTaskVtXY::CreateOutputObjects()
{
  // Create histograms
  // Called once

   fHistVtx = new TProfile("fHistVtx","Xvert-RMS vs sigmaX(NC)",100,0.,0.05,-1.,1.,"s");
   fHistVty = new TProfile("fHistVty","Yvert-RMS vs sigmaY(NC)",100,0.,0.05,-1.,1.,"s");
   fHistVtx->GetXaxis()->SetTitle("sigmaX(NContributors)");
   fHistVtx->GetYaxis()->SetTitle("Xv-av&spread");
   fHistVty->GetXaxis()->SetTitle("sigmaY(NContributors)");
   fHistVty->GetYaxis()->SetTitle("Yv-av&spread");
   fHistVtx->SetMarkerStyle(kFullCircle);
   fHistVty->SetMarkerStyle(kOpenCircle);
   //
   fList = new TList();
   fList->Add(fHistVtx);
   fList->Add(fHistVty);
}

//________________________________________________________________________
void AliAnalysisTaskVtXY::Exec(Option_t *) 
{
  // Main loop
  // Called for each event
  // 0 condition / existence of ESD
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }
  Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
  Int_t Ngood = 0;
  // Track loop kept in case of need
  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
    AliESDtrack* track = fESD->GetTrack(iTracks);
    if (!track) {
      Printf("ERROR: Could not receive track %d", iTracks);
      continue;
    }
    Ngood++;
  } //end track loop 
  // redo the ESD vertex withour constraint
  AliVertexerTracks vertexer(fESD->GetMagneticField());
  vertexer.SetITSMode();
  AliESDVertex *vertex = vertexer.FindPrimaryVertex(fESD);
  // ***************************************************
  //  const AliESDVertex *vtxESD = 0;
  // if(fESD) {vtxESD = fESD->GetPrimaryVertexTracks();}
  Double_t ESDvtx[3] ={999.,999.,9999.};
  Double_t XRes,YRes ;
  Double_t sig[3];
  // Cond 1.vertex has to exist 
   if (!vertex) {
     Printf("ERROR: vertex not available");
     return;
    }
   // Cond 2 - vertexer has not failed
   if(!vertex->GetStatus()){
     Printf("WARNING: vertexer has failed");
       return;
   }
     TString vtitle = vertex->GetTitle();
   // Cond 3 --> vertexer type check
     if(!vtitle.Contains("VertexerTracksNoConstraint")){
       Printf("WARNING: not VertexerTracksNoConstraint");
      return;
   }
   ULong64_t TrigMask=0;
   TrigMask = fESD->GetTriggerMask();
   // Cond.4 Trigger mask eventually - not activated
   if(TrigMask==0){
     Printf("Trigger Mask = 0");
     // return;
    }
    if(vertex){
    vertex->GetXYZ(ESDvtx);
    XRes = vertex->GetXRes();
    YRes = vertex->GetYRes();
    vertex->GetSigmaXYZ(sig);
    fHistVtx->Fill(XRes,ESDvtx[0]);
    fHistVty->Fill(YRes,ESDvtx[1]);
    PostData(0,fList);
    }
  //
}      

//________________________________________________________________________
void AliAnalysisTaskVtXY::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fList = dynamic_cast<TList*>(GetOutputData(0));
  if(!fList){
    Printf("ERROR: fList not available");
    return;
  }
  fHistVtx = (TProfile *)fList->At(0);
  fHistVty = (TProfile *)fList->At(1);
  TCanvas *c1 = new TCanvas("AliAnalysisTaskVtXY","Vtx analysis",10,10,800,800);
  c1->SetFillColor(10); c1->SetHighLightColor(10);
  c1->Divide(1,2);
  c1->cd(1);
  fHistVtx->DrawCopy("");
  c1->cd(2);
  fHistVty->DrawCopy("");
  // derive histos with fits from the above profiles
  Double_t spread = 0.;
  Double_t content = 0.;
  Double_t entries = 0.;
  Double_t error = 0.;
  TH1D * hsx = new TH1D("hsx","Xvtx spreads",100,0.,0.050);
  //fill
  for(Int_t i=0; i<100; i++){
    spread = fHistVtx->GetBinError(i);
    content = fHistVtx->GetBinContent(i);
    entries = fHistVtx->GetBinEntries(i);
    error = 0.;
    if(entries<10){content = 0.; spread = 0.;}
    if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
    hsx->SetBinContent(i,spread);
    hsx->SetBinError(i,error);
  }
  hsx->GetXaxis()->SetTitle("sigma-Xvert");
  hsx->GetYaxis()->SetTitle("RMS(Xv) [cm]");
  hsx->GetYaxis()->SetTitleOffset(-0.3);
  hsx->GetYaxis()->SetRangeUser(0.,0.15);
  hsx->SetMarkerColor(3);
  hsx->SetMarkerStyle(20);
  /*TCanvas *cx = */new TCanvas("cx","",50,50,800,800);
  gStyle->SetOptFit(111);
  TPad *px = new TPad("px","",0,0,1,1);
  px->Draw();
  px->cd();
  px->SetFillColor(42);
  px->SetFrameFillColor(10);
  hsx->Fit("pol3","R","",0.001,0.02);
  hsx->Draw();
  // VtY
  TH1D * hsy = new TH1D("hsy","Yvtx spreads",100,0.,0.050);
  //fill
  for(Int_t i=0; i<100; i++){
    spread = fHistVty->GetBinError(i);
    content = fHistVty->GetBinContent(i);
    entries = fHistVty->GetBinEntries(i);
    error = 0.;
    if(entries<10){content = 0.; spread = 0.;}
    if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
    hsy->SetBinContent(i,spread);
    hsy->SetBinError(i,error);
  }
  hsy->GetXaxis()->SetTitle("sigma-Yvert");
  hsy->GetYaxis()->SetTitle("RMS(Yv) [cm]");
  hsy->GetYaxis()->SetTitleOffset(-0.3);
  hsy->GetYaxis()->SetRangeUser(0.,0.15);
  hsy->SetMarkerColor(2);
  hsy->SetMarkerStyle(20);
  /*TCanvas *cy = */new TCanvas("cy","",100,100,800,800);
  TPad *py = new TPad("py","",0,0,1,1);
  py->Draw();
  py->cd();
  py->SetFillColor(42);
  py->SetFrameFillColor(10);
  hsy->Fit("pol3","R","",0.001,0.02);
  hsy->Draw();
}
Commit	Line	Data
b514f4be	1	/**************************************************************************
	2	* Copyright(c) 1998-2003, 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	// Implementation of the beam interaction spot location and size estimate
	18	// via VertexerTracksNoConstraint resolution extrapolation
	19	// to the infinit multiplicity alias zero resolution
	20	//
	21	// Origin: AliAnalysisTaskVtXY
	22	// A. Jacholkowski, Catania
	23	// adam.jacholkowski@cern.ch
	24	//-----------------------------------------------------------------
	25	#include "TChain.h"
	26	#include "TTree.h"
	27	#include "TStyle.h"
	28	#include "TH2F.h"
	29	#include "TProfile.h"
	30	#include "TCanvas.h"
	31
	32	#include "AliAnalysisTask.h"
	33	#include "AliAnalysisManager.h"
	34
	35	#include "AliESDEvent.h"
	36	#include "AliESDInputHandler.h"
	37
	38	#include "AliAnalysisTaskVtXY.h"
	39
	40	#include "AliESDVertex.h"
	41	#include "AliVertexerTracks.h"
	42	ClassImp(AliAnalysisTaskVtXY)
	43
	44	//________________________________________________________________________
	45	AliAnalysisTaskVtXY::AliAnalysisTaskVtXY(const char *name)
	46	: AliAnalysisTask(name, ""),
	47	fESD(0),
	48	fList(0),
	49	fHistVtx(0),
	50	fHistVty(0)
	51	{
	52	// Constructor
	53
	54	// Define input and output slots here
	55	// Input slot #0 works with a TChain
	56	DefineInput(0, TChain::Class());
	57	DefineOutput(0, TList::Class());
	58	// Output slot #0 writes into a TList container
	59	}
	60
	61	//________________________________________________________________________
	62	void AliAnalysisTaskVtXY::ConnectInputData(Option_t *)
	63	{
	64	// Connect ESD or AOD here
65	// Called once
66
67	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
68	if (!tree) {
69	Printf("ERROR: Could not read chain from input slot 0");
70	} else {
71	// Disable all branches and enable only the needed ones
72	// The next two lines are different when data produced as AliESDEvent is read
73	//tree->SetBranchStatus("*", kFALSE);
74	//tree->SetBranchStatus("Tracks.*", kTRUE);
75
76	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
77
78	if (!esdH) {
79	Printf("ERROR: Could not get ESDInputHandler");
80	} else
81	fESD = esdH->GetEvent();
82	}
83	}
84
85	//________________________________________________________________________
86	void AliAnalysisTaskVtXY::CreateOutputObjects()
87	{
88	// Create histograms
89	// Called once
90
91	fHistVtx = new TProfile("fHistVtx","Xvert-RMS vs sigmaX(NC)",100,0.,0.05,-1.,1.,"s");
92	fHistVty = new TProfile("fHistVty","Yvert-RMS vs sigmaY(NC)",100,0.,0.05,-1.,1.,"s");
93	fHistVtx->GetXaxis()->SetTitle("sigmaX(NContributors)");
94	fHistVtx->GetYaxis()->SetTitle("Xv-av&spread");
95	fHistVty->GetXaxis()->SetTitle("sigmaY(NContributors)");
96	fHistVty->GetYaxis()->SetTitle("Yv-av&spread");
97	fHistVtx->SetMarkerStyle(kFullCircle);
98	fHistVty->SetMarkerStyle(kOpenCircle);
99	//
100	fList = new TList();
101	fList->Add(fHistVtx);
102	fList->Add(fHistVty);
103	}
104
105	//________________________________________________________________________
106	void AliAnalysisTaskVtXY::Exec(Option_t *)
107	{
108	// Main loop
109	// Called for each event
110	// 0 condition / existence of ESD
111	if (!fESD) {
112	Printf("ERROR: fESD not available");
113	return;
114	}
115	Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
116	Int_t Ngood = 0;
117	// Track loop kept in case of need
118	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
119	AliESDtrack* track = fESD->GetTrack(iTracks);
120	if (!track) {
121	Printf("ERROR: Could not receive track %d", iTracks);
122	continue;
123	}
124	Ngood++;
125	} //end track loop
126	// redo the ESD vertex withour constraint
127	AliVertexerTracks vertexer(fESD->GetMagneticField());
128	vertexer.SetITSMode();
129	AliESDVertex *vertex = vertexer.FindPrimaryVertex(fESD);
130	// ***************************************************
131	// const AliESDVertex *vtxESD = 0;
132	// if(fESD) {vtxESD = fESD->GetPrimaryVertexTracks();}
133	Double_t ESDvtx[3] ={999.,999.,9999.};
134	Double_t XRes,YRes ;
135	Double_t sig[3];
136	// Cond 1.vertex has to exist
137	if (!vertex) {
138	Printf("ERROR: vertex not available");
139	return;
140	}
141	// Cond 2 - vertexer has not failed
142	if(!vertex->GetStatus()){
143	Printf("WARNING: vertexer has failed");
144	return;
145	}
146	TString vtitle = vertex->GetTitle();
147	// Cond 3 --> vertexer type check
148	if(!vtitle.Contains("VertexerTracksNoConstraint")){
149	Printf("WARNING: not VertexerTracksNoConstraint");
150	return;
151	}
152	ULong64_t TrigMask=0;
153	TrigMask = fESD->GetTriggerMask();
154	// Cond.4 Trigger mask eventually - not activated
155	if(TrigMask==0){
156	Printf("Trigger Mask = 0");
157	// return;
158	}
159	if(vertex){
160	vertex->GetXYZ(ESDvtx);
161	XRes = vertex->GetXRes();
162	YRes = vertex->GetYRes();
163	vertex->GetSigmaXYZ(sig);
164	fHistVtx->Fill(XRes,ESDvtx[0]);
165	fHistVty->Fill(YRes,ESDvtx[1]);
166	PostData(0,fList);
167	}
168	//
169	}
170
171	//________________________________________________________________________
172	void AliAnalysisTaskVtXY::Terminate(Option_t *)
173	{
174	// Draw result to the screen
175	// Called once at the end of the query
176	fList = dynamic_cast<TList*>(GetOutputData(0));
177	if(!fList){
178	Printf("ERROR: fList not available");
179	return;
180	}
d837bbcf	181	fHistVtx = (TProfile *)fList->At(0);
d837bbcf	182	fHistVty = (TProfile *)fList->At(1);
b514f4be	183	TCanvas *c1 = new TCanvas("AliAnalysisTaskVtXY","Vtx analysis",10,10,800,800);
	184	c1->SetFillColor(10); c1->SetHighLightColor(10);
	185	c1->Divide(1,2);
	186	c1->cd(1);
	187	fHistVtx->DrawCopy("");
	188	c1->cd(2);
	189	fHistVty->DrawCopy("");
	190	// derive histos with fits from the above profiles
	191	Double_t spread = 0.;
	192	Double_t content = 0.;
	193	Double_t entries = 0.;
	194	Double_t error = 0.;
	195	TH1D * hsx = new TH1D("hsx","Xvtx spreads",100,0.,0.050);
	196	//fill
	197	for(Int_t i=0; i<100; i++){
	198	spread = fHistVtx->GetBinError(i);
	199	content = fHistVtx->GetBinContent(i);
	200	entries = fHistVtx->GetBinEntries(i);
	201	error = 0.;
	202	if(entries<10){content = 0.; spread = 0.;}
	203	if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
	204	hsx->SetBinContent(i,spread);
	205	hsx->SetBinError(i,error);
	206	}
	207	hsx->GetXaxis()->SetTitle("sigma-Xvert");
	208	hsx->GetYaxis()->SetTitle("RMS(Xv) [cm]");
	209	hsx->GetYaxis()->SetTitleOffset(-0.3);
	210	hsx->GetYaxis()->SetRangeUser(0.,0.15);
	211	hsx->SetMarkerColor(3);
	212	hsx->SetMarkerStyle(20);
	213	/TCanvas cx = */new TCanvas("cx","",50,50,800,800);
	214	gStyle->SetOptFit(111);
	215	TPad *px = new TPad("px","",0,0,1,1);
	216	px->Draw();
	217	px->cd();
	218	px->SetFillColor(42);
	219	px->SetFrameFillColor(10);
	220	hsx->Fit("pol3","R","",0.001,0.02);
	221	hsx->Draw();
	222	// VtY
	223	TH1D * hsy = new TH1D("hsy","Yvtx spreads",100,0.,0.050);
	224	//fill
	225	for(Int_t i=0; i<100; i++){
	226	spread = fHistVty->GetBinError(i);
	227	content = fHistVty->GetBinContent(i);
	228	entries = fHistVty->GetBinEntries(i);
	229	error = 0.;
	230	if(entries<10){content = 0.; spread = 0.;}
	231	if(entries>0){error = 0.0005 + spread/TMath::Sqrt(entries);}
	232	hsy->SetBinContent(i,spread);
	233	hsy->SetBinError(i,error);
	234	}
	235	hsy->GetXaxis()->SetTitle("sigma-Yvert");
	236	hsy->GetYaxis()->SetTitle("RMS(Yv) [cm]");
	237	hsy->GetYaxis()->SetTitleOffset(-0.3);
	238	hsy->GetYaxis()->SetRangeUser(0.,0.15);
	239	hsy->SetMarkerColor(2);
	240	hsy->SetMarkerStyle(20);
	241	/TCanvas cy = */new TCanvas("cy","",100,100,800,800);
	242	TPad *py = new TPad("py","",0,0,1,1);
	243	py->Draw();
	244	py->cd();
	245	py->SetFillColor(42);
	246	py->SetFrameFillColor(10);
247	hsy->Fit("pol3","R","",0.001,0.02);
248	hsy->Draw();
249	}