Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskVtXY.cxx
CommitLineData
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"
42ClassImp(AliAnalysisTaskVtXY)
43
44//________________________________________________________________________
45AliAnalysisTaskVtXY::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//________________________________________________________________________
62void 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//________________________________________________________________________
86void 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//________________________________________________________________________
106void 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//________________________________________________________________________
172void 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);
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}