All FMD corrections are now divided into the analysis + adding new corrections
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskDensity.cxx
CommitLineData
3bb122c7 1
2#include <TROOT.h>
3#include <TSystem.h>
4#include <TInterpreter.h>
5#include <TChain.h>
6#include <TFile.h>
7#include <TList.h>
8#include <iostream>
9#include "TAxis.h"
10#include "TH2F.h"
f58a4769 11#include "TF1.h"
3bb122c7 12#include "AliFMDAnalysisTaskDensity.h"
13#include "AliAnalysisManager.h"
14#include "AliESDFMD.h"
15#include "AliESDEvent.h"
16#include "AliAODEvent.h"
17#include "AliAODHandler.h"
18#include "AliMCEventHandler.h"
19#include "AliStack.h"
20#include "AliESDVertex.h"
21#include "TMath.h"
22#include "AliFMDAnaParameters.h"
78f6f750 23//#include "AliFMDParameters.h"
24//#include "AliFMDGeometry.h"
25//#include "AliFMDRing.h"
3bb122c7 26
27ClassImp(AliFMDAnalysisTaskDensity)
28
29//_____________________________________________________________________
30AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
31: fDebug(0),
8dc7c4c2 32 fOutputList(),
c78bc12b 33 fESD(0x0),
7c3e5162 34 fVertexString(),
35 fVertex(0),
bb8a464f 36 fStandalone(kTRUE),
37 fStatus(kTRUE)
3bb122c7 38{
39 // Default constructor
7c3e5162 40 DefineInput (0, AliESDFMD::Class());
41 DefineInput (1, AliESDVertex::Class());
3bb122c7 42 DefineOutput(0,TList::Class());
43}
44//_____________________________________________________________________
7c3e5162 45AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
3bb122c7 46 AliAnalysisTask(name, "Density"),
47 fDebug(0),
7c3e5162 48 fOutputList(0),
c78bc12b 49 fESD(0x0),
7c3e5162 50 fVertexString(),
51 fVertex(0),
bb8a464f 52 fStandalone(kTRUE),
53 fStatus(kTRUE)
3bb122c7 54{
7c3e5162 55 fStandalone = SE;
56 if(fStandalone) {
57 DefineInput (0, AliESDFMD::Class());
58 DefineInput (1, AliESDVertex::Class());
59 DefineOutput(0, TList::Class());
60 }
f58a4769 61
3bb122c7 62}
63//_____________________________________________________________________
64void AliFMDAnalysisTaskDensity::CreateOutputObjects()
65{
66 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
3bb122c7 67
7c3e5162 68 if(!fOutputList)
69 fOutputList = new TList();
70 fOutputList->SetName("density_list");
3bb122c7 71
bb8a464f 72 fOutputList->Add(&fVertexString);
73
3bb122c7 74 TH2F* hMult = 0;
75
76 Int_t nVtxbins = pars->GetNvtxBins();
77
78 for(Int_t det =1; det<=3;det++)
79 {
3bb122c7 80 Int_t nRings = (det==1 ? 1 : 2);
81 for(Int_t ring = 0;ring<nRings;ring++)
82 {
83 Char_t ringChar = (ring == 0 ? 'I' : 'O');
84 Int_t nSec = (ring == 0 ? 20 : 40);
85
3bb122c7 86 for(Int_t i = 0; i< nVtxbins; i++) {
87 TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
88
89 hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
90 hBg->GetNbinsX(),
91 hBg->GetXaxis()->GetXmin(),
92 hBg->GetXaxis()->GetXmax(),
93 nSec, 0, 2*TMath::Pi());
3bb122c7 94
9f55be54 95 fOutputList->Add(hMult);
3bb122c7 96 }
97 }
98 }
8dc7c4c2 99
bb8a464f 100
7c3e5162 101
3bb122c7 102
103}
104//_____________________________________________________________________
105void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
106{
7c3e5162 107 if(fStandalone) {
108 fESD = (AliESDFMD*)GetInputData(0);
109 fVertex = (AliESDVertex*)GetInputData(1);
110 }
3bb122c7 111}
f58a4769 112
113
114
115
3bb122c7 116//_____________________________________________________________________
117void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
118{
119 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
78f6f750 120 // AliFMDGeometry* geo = AliFMDGeometry::Instance();
3bb122c7 121
7c3e5162 122 //AliESDFMD* fmd = fESD->GetFMDData();
3bb122c7 123
3bb122c7 124 Double_t vertex[3];
7c3e5162 125 fVertex->GetXYZ(vertex);
8dc823cc 126 // Z Vtx cut
bb8a464f 127 if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
128 fStatus = kFALSE;
3bb122c7 129 return;
bb8a464f 130 }
131 else
132 fStatus = kTRUE;
133
3bb122c7 134 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
135 Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
136
137 Int_t vtxbin = (Int_t)vertexBinDouble;
f58a4769 138
8dc7c4c2 139 fVertexString.SetString(Form("%d",vtxbin));
8dc823cc 140 //Reset everything
3bb122c7 141 for(UShort_t det=1;det<=3;det++) {
3bb122c7 142 Int_t nRings = (det==1 ? 1 : 2);
143 for (UShort_t ir = 0; ir < nRings; ir++) {
9f55be54 144 Char_t ring = (ir == 0 ? 'I' : 'O');
145 TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
3bb122c7 146 hMult->Reset();
147 }
148
149 }
150
151
152 for(UShort_t det=1;det<=3;det++) {
3bb122c7 153 Int_t nRings = (det==1 ? 1 : 2);
154 for (UShort_t ir = 0; ir < nRings; ir++) {
bb8a464f 155
3bb122c7 156 Char_t ring = (ir == 0 ? 'I' : 'O');
9f55be54 157 TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
158
3bb122c7 159 UShort_t nsec = (ir == 0 ? 20 : 40);
160 UShort_t nstr = (ir == 0 ? 512 : 256);
bb8a464f 161
3bb122c7 162 for(UShort_t sec =0; sec < nsec; sec++) {
163 for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162 164 Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
9f55be54 165
bb8a464f 166 if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
9f55be54 167
f58a4769 168 Float_t phi = pars->GetPhiFromSector(det,ring,sec);
169 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
41bad769 170
b85ea106 171 Float_t mult_cut = 0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2;
172 if(ring == 'I')
173 mult_cut = 0.10;
174 //Float_t mult_cut = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
bb8a464f 175 Float_t nParticles = 0;
6289b3e8 176 if(fESD->GetUniqueID() == kTRUE) {
177 //proton + proton
cc066cb9 178
179 if(mult > mult_cut) {
6289b3e8 180 nParticles = 1;
cc066cb9 181
182 }
6289b3e8 183 }
184 else {
185
186 //Pb+Pb
187 Float_t mpv = pars->GetMPV(det,ring,eta);
188 Float_t sigma = pars->GetSigma(det,ring,eta);
189 Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
190 Float_t beta = pars->Get3MIPWeight(det,ring,eta);
191
192 Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
193 alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
194 beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
195 Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
196 2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
197 3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
198
199
cc066cb9 200 if(mult > mult_cut) {
6289b3e8 201 if(sumCor) nParticles = weight / sumCor;
202 else nParticles = 1;
cc066cb9 203
6289b3e8 204 }
205 //std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl;
206
207 }
208
209
cc066cb9 210
211
212
1282ce49 213 Float_t correction = GetAcceptanceCorrection(ring,strip);
f58a4769 214
215 //std::cout<<"before "<<correction<<std::endl;
78f6f750 216 if(fESD->GetUniqueID() == kTRUE) {
4fb49e43 217 TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring);
41bad769 218
4fb49e43 219 if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0)
b85ea106 220 correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta));
41bad769 221
78f6f750 222 }
f58a4769 223
6289b3e8 224 if(correction) nParticles = nParticles / correction;
cc066cb9 225 if(nParticles > 0)
226 hMult->Fill(eta,phi,nParticles);
bb8a464f 227
3bb122c7 228
229 }
230 }
bb8a464f 231
3bb122c7 232 }
bb8a464f 233
3bb122c7 234
235
236 }
78f6f750 237
238
7c3e5162 239 if(fStandalone) {
240 PostData(0, fOutputList);
241 }
3bb122c7 242
243}
244//_____________________________________________________________________
1282ce49 245Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
246{
78f6f750 247 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
1282ce49 248
78f6f750 249 //AliFMDRing fmdring(ring);
250 //fmdring.Init();
251 Float_t rad = pars->GetMaxR(ring)-pars->GetMinR(ring);
252 Float_t nstrips = (ring == 'I' ? 512 : 256);
253 Float_t segment = rad / nstrips;
254 Float_t radius = pars->GetMinR(ring) + segment*strip;
255
256 Float_t basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
257 Float_t basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
1282ce49 258 Float_t basearea = basearea1 - basearea2;
259
78f6f750 260 Float_t area1 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
261 Float_t area2 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
1282ce49 262 Float_t area = area1 - area2;
263
264 Float_t correction = area/basearea;
265
266 return correction;
267}
f58a4769 268//_____________________________________________________________________
3bb122c7 269//
270//EOF
271//