]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis/AliFMDAnalysisTaskSharing.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Compilation warnings
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskSharing.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>
bb8a464f 9#include <TMath.h>
78f6f750 10//#include "AliFMDDebug.h"
3bb122c7 11#include "AliFMDAnalysisTaskSharing.h"
12#include "AliAnalysisManager.h"
13#include "AliESDFMD.h"
78f6f750 14//#include "AliFMDGeometry.h"
3bb122c7 15#include "AliMCEventHandler.h"
b64db9b1 16#include "AliMCEvent.h"
3bb122c7 17#include "AliStack.h"
18#include "AliESDVertex.h"
6289b3e8 19#include "AliMultiplicity.h"
3bb122c7 20#include "AliFMDAnaParameters.h"
0210c863 21#include "TH1F.h"
22#include "TObjString.h"
bb0a45c3 23//#include "/home/canute/ALICE/AliRoot/PWG0/AliPWG0Helper.h"
78f6f750 24//#include "AliFMDParameters.h"
b64db9b1 25#include "AliGenEventHeader.h"
da0805e2 26#include "AliGenPythiaEventHeader.h"
b64db9b1 27#include "AliHeader.h"
28#include "AliStack.h"
29#include "AliMCParticle.h"
0a2f2742 30#include "AliFMDStripIndex.h"
da0805e2 31#include "AliESDVZERO.h"
32#include "AliESDtrack.h"
3bb122c7 33
0210c863 34// This is the task to do the FMD sharing or hit merging.
35// It reads the input ESDFMD data and posts an ESDFMD object to
36// the tasks that must be performed after this task ie.
37// Density, BackgroundCorrection and Dndeta.
38// Author: Hans Hjersing Dalsgaard, hans.dalsgaard@cern.ch
39
40
3bb122c7 41ClassImp(AliFMDAnalysisTaskSharing)
42
43//_____________________________________________________________________
44AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
45: fDebug(0),
46 fESD(0x0),
8dc7c4c2 47 foutputESDFMD(),
c78bc12b 48 fSharedThis(kFALSE),
7c3e5162 49 fSharedPrev(kFALSE),
a38d1871 50 fDiagList(0),
7c3e5162 51 fStandalone(kTRUE),
bb8a464f 52 fEsdVertex(0),
0a2f2742 53 fStatus(kTRUE),
aa303f0c 54 fLastTrackByStrip(0),
55 fLastOrbit(0)
3bb122c7 56{
57 // Default constructor
58 DefineInput (0, AliESDEvent::Class());
7c3e5162 59 DefineOutput(0, AliESDFMD::Class());
60 DefineOutput(1, AliESDVertex::Class());
61 DefineOutput(2, AliESDEvent::Class());
62 DefineOutput(3, TList::Class());
3bb122c7 63}
64//_____________________________________________________________________
7c3e5162 65AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
3bb122c7 66 AliAnalysisTask(name, "AnalysisTaskFMD"),
67 fDebug(0),
68 fESD(0x0),
8dc7c4c2 69 foutputESDFMD(),
c78bc12b 70 fSharedThis(kFALSE),
7c3e5162 71 fSharedPrev(kFALSE),
a38d1871 72 fDiagList(0),
7c3e5162 73 fStandalone(kTRUE),
bb8a464f 74 fEsdVertex(0),
0a2f2742 75 fStatus(kTRUE),
da0805e2 76 fLastTrackByStrip(0),
aa303f0c 77 fLastOrbit(0)
3bb122c7 78{
0210c863 79 // named constructor
7c3e5162 80 fStandalone = SE;
81 if(fStandalone) {
82 DefineInput (0, AliESDEvent::Class());
83 DefineOutput(0, AliESDFMD::Class());
84 DefineOutput(1, AliESDVertex::Class());
85 DefineOutput(2, AliESDEvent::Class());
86 DefineOutput(3, TList::Class());
87 }
3bb122c7 88}
89//_____________________________________________________________________
daedf077 90Float_t AliFMDAnalysisTaskSharing::GetVtxEfficiencyFromData(){
91
92 TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
93 TH1F* hEventsTr = (TH1F*)fDiagList->FindObject("hEventsTr");
94
95 if(hEventsTr->GetEntries() != 0 && hEventsVtx->GetEntries() !=0 && hEventsTr->GetEntries() != hEventsVtx->GetEntries())
96 return hEventsVtx->GetEntries() / hEventsTr->GetEntries();
97 else return -1;
98
99}
059c7c6b 100//_____________________________________________________________________
101Float_t AliFMDAnalysisTaskSharing::GetNSDVtxEfficiencyFromData(){
102
103 TH1F* hEventsNSDVtx = (TH1F*)fDiagList->FindObject("hEventsNSDVtx");
104 TH1F* hEventsNSD = (TH1F*)fDiagList->FindObject("hEventsNSD");
105
106 if(hEventsNSD->GetEntries() != 0 && hEventsNSDVtx->GetEntries() !=0 && hEventsNSD->GetEntries() != hEventsNSDVtx->GetEntries())
107 return hEventsNSDVtx->GetEntries() / hEventsNSD->GetEntries();
108 else return -1;
109
110}
111
112
daedf077 113//_____________________________________________________________________
3bb122c7 114void AliFMDAnalysisTaskSharing::CreateOutputObjects()
115{
0210c863 116 // Create the output objects
7c3e5162 117 if(!foutputESDFMD)
118 foutputESDFMD = new AliESDFMD();
119
120 if(!fEsdVertex)
121 fEsdVertex = new AliESDVertex();
122 //Diagnostics
a38d1871 123 if(!fDiagList)
124 fDiagList = new TList();
125
126 fDiagList->SetName("Sharing diagnostics");
b64db9b1 127
128 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
129 TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', 0);
130 TH1F* hPrimary = new TH1F("hMultvsEtaNoCuts","hMultvsEtaNoCuts",
131 hBg->GetNbinsX(),
132 hBg->GetXaxis()->GetXmin(),
133 hBg->GetXaxis()->GetXmax());
134 hPrimary->Sumw2();
135 fDiagList->Add(hPrimary);
059c7c6b 136
137 TH1F* hPrimaryNSD = new TH1F("hMultvsEtaNSDNoCuts","hMultvsEtaNSDNoCuts",
138 hBg->GetNbinsX(),
139 hBg->GetXaxis()->GetXmin(),
140 hBg->GetXaxis()->GetXmax());
141 hPrimaryNSD->Sumw2();
142 fDiagList->Add(hPrimaryNSD);
143
da0805e2 144 TH1F* hXvtx = new TH1F("hXvtx","x vertex distribution",100,-2,2);
145 TH1F* hYvtx = new TH1F("hYvtx","y vertex distribution",100,-2,2);
850ad94d 146 TH1F* hZvtx = new TH1F("hZvtx","z vertex distribution",4*pars->GetNvtxBins(),-4*pars->GetVtxCutZ(),4*pars->GetVtxCutZ());
f55d559b 147
da0805e2 148 fDiagList->Add(hXvtx);
149 fDiagList->Add(hYvtx);
f55d559b 150 fDiagList->Add(hZvtx);
151
b64db9b1 152 TH1F* hPrimVertexBin = 0;
0a2f2742 153 TH1F* hHits = 0;
b64db9b1 154 for(Int_t i = 0; i< pars->GetNvtxBins(); i++) {
155
156 hPrimVertexBin = new TH1F(Form("primmult_NoCuts_vtxbin%d",i),
157 Form("primmult_NoCuts_vtxbin%d",i),
158 hBg->GetNbinsX(),
159 hBg->GetXaxis()->GetXmin(),
160 hBg->GetXaxis()->GetXmax());
161 hPrimVertexBin->Sumw2();
162 fDiagList->Add(hPrimVertexBin);
163
164 }
165
7c3e5162 166 for(Int_t det = 1; det<=3; det++) {
167 Int_t nRings = (det==1 ? 1 : 2);
168
169 for(Int_t iring = 0;iring<nRings; iring++) {
170 Char_t ringChar = (iring == 0 ? 'I' : 'O');
171 TH1F* hEdist = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
172 Form("Edist_before_sharing_FMD%d%c", det, ringChar),
bb8a464f 173 1000,0,25);
0210c863 174 TH1F* hEdistAfter = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
7c3e5162 175 Form("Edist_after_sharing_FMD%d%c", det, ringChar),
bb8a464f 176 1000,0,25);
177
178
b64db9b1 179 //TH1F* hNstripsHit = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
180 // Form("N_strips_hit_FMD%d%c",det,ringChar),
181 // 25,0,25);
a38d1871 182 fDiagList->Add(hEdist);
0210c863 183 fDiagList->Add(hEdistAfter);
b64db9b1 184 //fDiagList->Add(hNstripsHit);
0a2f2742 185
186 for(Int_t i = 0; i< pars->GetNvtxBins(); i++) {
187 hHits = new TH1F(Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
188 hBg->GetNbinsX(),
189 hBg->GetXaxis()->GetXmin(),
190 hBg->GetXaxis()->GetXmax());
191 hHits->Sumw2();
192 fDiagList->Add(hHits);
7c3e5162 193
0a2f2742 194 }
195
7c3e5162 196 }
197 }
da0805e2 198 TH2F* hCorrelationFMDSPDhits = new TH2F("hCorrelationFMDSPDhits","hCorrelationFMDSPDhits;SPD;FMD ",100,0,200,100,0,500);
199 TH2F* hCorrelationFMDSPD = new TH2F("hCorrelationFMDSPD","hCorrelationFMDSPD;SPD;FMD ",100,0,200,100,0,500);
200 TH2F* hCorrelationFMD1SPD = new TH2F("hCorrelationFMD1SPD","hCorrelationFMD1SPD;SPD;FMD1 ",100,0,200,100,0,200);
201 TH2F* hCorrelationFMD2ISPD = new TH2F("hCorrelationFMD2ISPD","hCorrelationFMD2ISPD;SPD;FMD2I ",100,0,200,100,0,200);
202 TH2F* hCorrelationFMD2OSPD = new TH2F("hCorrelationFMD2OSPD","hCorrelationFMD2OSPD;SPD;FMD2O ",100,0,200,100,0,200);
203 TH2F* hCorrelationFMD3ISPD = new TH2F("hCorrelationFMD3ISPD","hCorrelationFMD3ISPD;SPD;FMD3I ",100,0,200,100,0,200);
204 TH2F* hCorrelationFMD3OSPD = new TH2F("hCorrelationFMD3OSPD","hCorrelationFMD3OSPD;SPD;FMD3O ",100,0,200,100,0,200);
205 TH2F* hCorrelationFMDGoodtracks = new TH2F("hCorrelationFMDGoodtracks","hCorrelationGoodtracks;good tracks;FMD ",100,0,200,100,0,200);
206 TH2F* hCorrelationFMDBadtracks = new TH2F("hCorrelationFMDBadtracks","hCorrelationBadtracks;bad tracks;FMD ",100,0,200,100,0,200);
207 TH2F* hCorrelationGoodbadtracks = new TH2F("hCorrelationGoodbadtracks","hCorrelationGoodbadtracks;good tracks;bad tracks ",100,0,200,100,0,200);
208 TH2F* hCorrelationSPDTracklets = new TH2F("hCorrelationSPDTracklets","hCorrelationSPDTracklets;hits ; tracklets ",100,0,500,100,0,200);
209 TH2F* hCorrelationClustersTracklets = new TH2F("hCorrelationClustersTracklets","hCorrelationClustersTracklets;clusters ; tracklets ",500,0,500,100,0,100);
210 TH2F* hCorrelationHitsRadius = new TH2F("hCorrelationHitsRadius","hCorrelationHitsRadius;hits ; radius ",100,0,500,100,0,10);
211 TH2F* hCorrelationHitsX = new TH2F("hCorrelationHitsX","hCorrelationHitsX;hits ; X ",100,0,500,100,-5,5);
212 TH2F* hCorrelationHitsY = new TH2F("hCorrelationHitsY","hCorrelationHitsY;hits ; Y ",100,0,500,100,-5,5);
213 fDiagList->Add(hCorrelationHitsRadius);
214 fDiagList->Add(hCorrelationHitsX);
215 fDiagList->Add(hCorrelationHitsY);
216 fDiagList->Add(hCorrelationFMDSPD);
217 fDiagList->Add(hCorrelationFMD1SPD);
218 fDiagList->Add(hCorrelationFMD2ISPD);
219 fDiagList->Add(hCorrelationFMD2OSPD);
220 fDiagList->Add(hCorrelationFMD3ISPD);
221 fDiagList->Add(hCorrelationFMD3OSPD);
222 fDiagList->Add(hCorrelationFMDGoodtracks);
223 fDiagList->Add(hCorrelationFMDBadtracks);
224 fDiagList->Add(hCorrelationGoodbadtracks);
225fDiagList->Add(hCorrelationFMDSPDhits);
226 fDiagList->Add(hCorrelationClustersTracklets);
227 fDiagList->Add(hCorrelationSPDTracklets);
228 TH2F* hCorrelationFMD23 = new TH2F("hCorrelationFMD23","hCorrelationFMD23;FMD2 ;FMD3 ",100,0,500,100,0,500);
229 TH2F* hCorrelationFMD2diff23 = new TH2F("hCorrelationFMD2diff23","hCorrelationFMD2diff23;FMD2 ;diff FMD23 ",100,0,100,100,0,100);
230 TH2F* hCorrelationFMD3diff23 = new TH2F("hCorrelationFMD3diff23","hCorrelationFMD3diff23;FMD3 ;diff FMD23 ",100,0,100,100,0,100);
231 TH2F* hCorrelationFMD1diff13 = new TH2F("hCorrelationFMD1diff13","hCorrelationFMD1diff13;FMD1 ;diff FMD13 ",100,0,100,100,0,100);
232 TH2F* hCorrelationFMD1diff12 = new TH2F("hCorrelationFMD1diff12","hCorrelationFMD1diff12;FMD1 ;diff FMD12 ",100,0,100,100,0,100);
233 TH2F* hCorrelationFMD12 = new TH2F("hCorrelationFMD12","hCorrelationFMD12;FMD1 ;FMD2 ",100,0,500,100,0,500);
234 TH2F* hCorrelationFMDBgCand = new TH2F("hCorrelationFMDBgCand","hCorrelationFMDBgCand;Bg Tr ;FMD ",100,0,100,500,0,500);
235
236 TH2F* hCorrelationFMDFlatTr = new TH2F("hCorrelationFMDFlatTr","hCorrelationFMDFlatTr;Bg Tr ;FMD ",100,0,100,500,0,500);
237 TH2F* hCorrelationFMDRatioFlatTr = new TH2F("hCorrelationFMDRatioFlatTr","hCorrelationFMDRatioFlatTr;Bg Tr ;FMD ",100,0,1,500,0,500);
238 fDiagList->Add(hCorrelationFMDBgCand);
239 fDiagList->Add(hCorrelationFMDFlatTr);
240 fDiagList->Add(hCorrelationFMDRatioFlatTr);
241 TH2F* hCorrelationFMDBgCandRelative = new TH2F("hCorrelationFMDBgCandRelative","hCorrelationFMDBgCandRelative;Bg Tr ;FMD ",100,0,2,500,0,500);
242 fDiagList->Add(hCorrelationFMDBgCandRelative);
243 fDiagList->Add(hCorrelationFMD2diff23);
244 fDiagList->Add(hCorrelationFMD3diff23);
245 fDiagList->Add(hCorrelationFMD1diff13);
246 fDiagList->Add(hCorrelationFMD1diff12);
247 fDiagList->Add(hCorrelationFMD23);
248 fDiagList->Add(hCorrelationFMD12);
249 TH2F* hTimeCorrelation = new TH2F("hCorrelationTime","hCorrelationTime ; time ; FMD hits",500,0,500,100,0,200);
250 fDiagList->Add(hTimeCorrelation);
251 TH1F* hHitDistribution = new TH1F("hFMDHitDistribution","hFMDHitDistribution ; FMD hits",500,0,500);
252
253 TH1F* hHitDistributionFMD1 = new TH1F("hFMDHitDistributionFMD1","hFMDHitDistributionFMD1 ; FMD1 hits",500,0,500);
254 TH1F* hHitDistributionFMD2I = new TH1F("hFMDHitDistributionFMD2I","hFMDHitDistributionFMD2I ; FMD2I hits",500,0,500);
255 TH1F* hHitDistributionFMD2O = new TH1F("hFMDHitDistributionFMD2O","hFMDHitDistributionFMD2O ; FMD2O hits",500,0,500);
256 TH1F* hHitDistributionFMD3I = new TH1F("hFMDHitDistributionFMD3I","hFMDHitDistributionFMD3I ; FMD3I hits",500,0,500);
257 TH1F* hHitDistributionFMD3O = new TH1F("hFMDHitDistributionFMD3O","hFMDHitDistributionFMD3O ; FMD3O hits",500,0,500);
258 TH1F* hTrVtxDistribution = new TH1F("hTrVtxDistribution","hTrVtxDistribution ; TrVtx",200,-500,500);
259 TH1F* hTrEtaDistribution = new TH1F("hTrEtaDistribution","hTrEtaDistribution ; TrEta",200,-9,9);
260 TH1F* hTrEtaDistribution2 = new TH1F("hTrEtaDistribution2","hTrEtaDistribution2 ; TrEta",200,-9,9);
261
262 TH1F* hFlatTracks = new TH1F("hFlatTracks","hFlatTracks ; Horizontal tracks",100,0,100);
aa303f0c 263
264 TH1F* hEnergyOfParticles = new TH1F("hEnergyOfParticles","hEnergyOfParticles",1000000,0,10);
265 fDiagList->Add(hEnergyOfParticles);
da0805e2 266 fDiagList->Add(hTrVtxDistribution);
267 fDiagList->Add(hTrEtaDistribution);
268 fDiagList->Add(hTrEtaDistribution2);
269 fDiagList->Add(hFlatTracks);
270 fDiagList->Add(hHitDistribution);
271 fDiagList->Add(hHitDistributionFMD1);
272 fDiagList->Add(hHitDistributionFMD2I);
273 fDiagList->Add(hHitDistributionFMD2O);
274 fDiagList->Add(hHitDistributionFMD3I);
275 fDiagList->Add(hHitDistributionFMD3O);
b64db9b1 276 TH1F* nMCevents = new TH1F("nMCEventsNoCuts","nMCEventsNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
059c7c6b 277 TH1F* nMCeventsNSD= new TH1F("nMCEventsNSDNoCuts","nMCEventsNSDNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
b64db9b1 278 fDiagList->Add(nMCevents);
059c7c6b 279 fDiagList->Add(nMCeventsNSD);
daedf077 280
059c7c6b 281 TH1F* hEventsVtx = new TH1F("hEventsVtx","hEventsVtx",pars->GetNvtxBins(),0,pars->GetNvtxBins());
282 TH1F* hEventsTr = new TH1F("hEventsTr","hEventsTr",pars->GetNvtxBins(),0,pars->GetNvtxBins());
283 TH1F* hEventsNSD = new TH1F("hEventsNSD","hEventsNSD",pars->GetNvtxBins(),0,pars->GetNvtxBins());
284 TH1F* hEventsNSDVtx = new TH1F("hEventsNSDVtx","hEventsNSDVtx",pars->GetNvtxBins(),0,pars->GetNvtxBins());
daedf077 285 fDiagList->Add(hEventsVtx);
286 fDiagList->Add(hEventsTr);
059c7c6b 287 fDiagList->Add(hEventsNSD);
288 fDiagList->Add(hEventsNSDVtx);
daedf077 289
3bb122c7 290}
291//_____________________________________________________________________
292void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
293{
0210c863 294 // connect the input data
7c3e5162 295 if(fStandalone)
296 fESD = (AliESDEvent*)GetInputData(0);
3bb122c7 297}
298//_____________________________________________________________________
299void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
300{
0210c863 301 //perform analysis on one event
3bb122c7 302 AliESD* old = fESD->GetAliESDOld();
303 if (old) {
304 fESD->CopyFromOldESD();
305 }
306
7c3e5162 307 foutputESDFMD->Clear();
541c19ed 308
309
da0805e2 310 Int_t delta = fESD->GetOrbitNumber() - fLastOrbit;
311 fLastOrbit = fESD->GetOrbitNumber();
312
313
b64db9b1 314 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
059c7c6b 315 Double_t vertex[3]={0,0,0};
f55d559b 316 Bool_t vtxStatus = pars->GetVertex(fESD,vertex);
bb8a464f 317 fEsdVertex->SetXYZ(vertex);
b64db9b1 318
059c7c6b 319 // std::cout<<vtxStatus<<" "<<vertex[0]<<" "<<vertex[1]<<" "<<vertex[2]<<std::endl;
320
b64db9b1 321 // Process primaries here to get true MC distribution
322 if(pars->GetProcessPrimary())
323 ProcessPrimary();
da0805e2 324 const AliMultiplicity* testmult = fESD->GetMultiplicity();
325 Int_t nTrackLets = testmult->GetNumberOfTracklets();
326 TH2F* hCorrelationClustersTracklets = (TH2F*)fDiagList->FindObject("hCorrelationClustersTracklets");
327 hCorrelationClustersTracklets->Fill(testmult->GetNumberOfSingleClusters(),nTrackLets);
b64db9b1 328
f55d559b 329
6289b3e8 330
059c7c6b 331 Bool_t isTriggered = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
332 Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
333
daedf077 334 Double_t delta2 = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
335 Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta2;
336 Int_t vtxbin = (Int_t)vertexBinDouble;
059c7c6b 337 TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
338 TH1F* hEventsTr = (TH1F*)fDiagList->FindObject("hEventsTr");
339 TH1F* hEventsNSD = (TH1F*)fDiagList->FindObject("hEventsNSD");
340 TH1F* hEventsNSDVtx = (TH1F*)fDiagList->FindObject("hEventsNSDVtx");
341
a82227ad 342 // if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
343 // vtxStatus = kFALSE;
344 // }
059c7c6b 345
346
daedf077 347 if(vtxStatus) hEventsVtx->Fill(vtxbin);
348 if(isTriggered) hEventsTr->Fill(vtxbin);
da0805e2 349
059c7c6b 350 if(vtxStatus && nsd) hEventsNSDVtx->Fill(vtxbin);
351 if(nsd) hEventsNSD->Fill(vtxbin);
352
da0805e2 353 if(!isTriggered || !vtxStatus ) {
bb8a464f 354 fStatus = kFALSE;
355 return;
356 }
357 else
358 fStatus = kTRUE;
b64db9b1 359
da0805e2 360 TH1F* hXvtx = (TH1F*)fDiagList->FindObject("hXvtx");
059c7c6b 361 if(vtxStatus ) hXvtx->Fill(vertex[0]);
da0805e2 362 TH1F* hYvtx = (TH1F*)fDiagList->FindObject("hYvtx");
059c7c6b 363 if(vtxStatus ) hYvtx->Fill(vertex[1]);
f55d559b 364 TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
365 hZvtx->Fill(vertex[2]);
da0805e2 366 //const AliMultiplicity* testmult = fESD->GetMultiplicity();
367 //std::cout<<vertex[2]<<std::endl;
368 //Int_t nTrackLets = testmult->GetNumberOfTracklets();
6289b3e8 369
a82227ad 370 if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
371 fStatus = kFALSE;
372 return;
373 }
da0805e2 374
6289b3e8 375 if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
376 else foutputESDFMD->SetUniqueID(kFALSE);
b64db9b1 377
3bb122c7 378 AliESDFMD* fmd = fESD->GetFMDData();
059c7c6b 379
3bb122c7 380 if (!fmd) return;
da0805e2 381 Int_t nHits[3][2] = {{0,0},{0,0},{0,0}};
382
3bb122c7 383 for(UShort_t det=1;det<=3;det++) {
384 Int_t nRings = (det==1 ? 1 : 2);
385 for (UShort_t ir = 0; ir < nRings; ir++) {
386 Char_t ring = (ir == 0 ? 'I' : 'O');
387 UShort_t nsec = (ir == 0 ? 20 : 40);
388 UShort_t nstr = (ir == 0 ? 512 : 256);
7c3e5162 389
a38d1871 390 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
7c3e5162 391
8dc7c4c2 392 for(UShort_t sec =0; sec < nsec; sec++) {
393 fSharedThis = kFALSE;
394 fSharedPrev = kFALSE;
0210c863 395
3bb122c7 396 for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162 397 foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
3bb122c7 398 Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
bb8a464f 399
8dc7c4c2 400 if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
7c3e5162 401
da0805e2 402
78f6f750 403 //Double_t eta = fmd->Eta(det,ring,sec,strip);
404 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
da0805e2 405
7c3e5162 406 hEdist->Fill(mult);
bb8a464f 407 if(fmd->IsAngleCorrected())
cc066cb9 408 mult = mult/TMath::Cos(Eta2Theta(eta));
0210c863 409 Float_t prevE = 0;
410 Float_t nextE = 0;
8dc823cc 411 if(strip != 0)
bb8a464f 412 if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
0210c863 413 prevE = fmd->Multiplicity(det,ring,sec,strip-1);
bb8a464f 414 if(fmd->IsAngleCorrected())
0210c863 415 prevE = prevE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
bb8a464f 416 }
8dc823cc 417 if(strip != nstr - 1)
bb8a464f 418 if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
0210c863 419 nextE = fmd->Multiplicity(det,ring,sec,strip+1);
bb8a464f 420 if(fmd->IsAngleCorrected())
0210c863 421 nextE = nextE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
bb8a464f 422 }
8dc823cc 423
0210c863 424 Float_t mergedEnergy = GetMultiplicityOfStrip(mult,eta,prevE,nextE,det,ring,sec,strip);
da0805e2 425 //if(mult> (pars->GetMPV(det,ring,eta) - pars->GetSigma(det,ring,eta)))
426 // mergedEnergy = mult;
427 //else mergedEnergy = 0;
428 if(mergedEnergy > 0.3 )
429 nHits[det-1][ir]++;
0210c863 430 foutputESDFMD->SetMultiplicity(det,ring,sec,strip,mergedEnergy);
bb8a464f 431 foutputESDFMD->SetEta(det,ring,sec,strip,eta);
7c3e5162 432
3bb122c7 433 }
434 }
435 }
436 }
da0805e2 437 //cluster cut
438 //if(testmult->GetNumberOfSingleClusters() > 15 + nTrackLets)
439 // {fStatus = kFALSE; std::cout<<"FMD : "<<nHits[0][0]<<" "<<nHits[1][0]<<" "<<nHits[1][1]<<" "<<nHits[2][0]<<" "<<nHits[2][1]<<" tracks "<<testmult->GetNumberOfSingleClusters()<<" "<<nTrackLets<<std::endl; return;}
440
441 TH2F* hCorrelationFMD23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD23");
442 TH2F* hCorrelationFMD12 = (TH2F*)fDiagList->FindObject("hCorrelationFMD12");
443 TH2F* hCorrelationFMD2diff23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD2diff23");
444 TH2F* hCorrelationFMD3diff23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD3diff23");
445 TH2F* hCorrelationFMD1diff13 = (TH2F*)fDiagList->FindObject("hCorrelationFMD1diff13");
446 TH2F* hCorrelationFMD1diff12 = (TH2F*)fDiagList->FindObject("hCorrelationFMD1diff12");
447
448 TH2F* hCorrelationFMDSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMDSPD");
449 TH2F* hCorrelationFMD1SPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD1SPD");
450 TH2F* hCorrelationFMD2ISPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD2ISPD");
451 TH2F* hCorrelationFMD2OSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD2OSPD");
452 TH2F* hCorrelationFMD3ISPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD3ISPD");
453 TH2F* hCorrelationFMD3OSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD3OSPD");
454
455 TH2F* hCorrelationFMDSPDhits = (TH2F*)fDiagList->FindObject("hCorrelationFMDSPDhits");
456
457 TH2F* hCorrelationSPDTracklets = (TH2F*)fDiagList->FindObject("hCorrelationSPDTracklets");
458 TH2F* hTimeCorrelation = (TH2F*)fDiagList->FindObject("hCorrelationTime");
459 TH2F* hHitsRadius = (TH2F*)fDiagList->FindObject("hCorrelationHitsRadius");
460 TH2F* hHitsX = (TH2F*)fDiagList->FindObject("hCorrelationHitsX");
461 TH2F* hHitsY = (TH2F*)fDiagList->FindObject("hCorrelationHitsY");
462 TH1F* hFMDHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistribution");
463 TH1F* hFMD1HitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD1");
464 TH1F* hFMD2IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2I");
465 TH1F* hFMD2OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2O");
466 TH1F* hFMD3IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3I");
467 TH1F* hFMD3OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3O");
468 TH1F* hCorrelationFMDGoodtracks = (TH1F*)fDiagList->FindObject("hCorrelationFMDGoodtracks");
469 TH1F* hCorrelationFMDBadtracks = (TH1F*)fDiagList->FindObject("hCorrelationFMDBadtracks");
470 TH1F* hCorrelationGoodbadtracks = (TH1F*)fDiagList->FindObject("hCorrelationGoodbadtracks");
471 TH2F* hCorrelationFMDBgCand = (TH2F*)fDiagList->FindObject("hCorrelationFMDBgCand");
472 TH2F* hCorrelationFMDBgCandRelative = (TH2F*)fDiagList->FindObject("hCorrelationFMDBgCandRelative");
473
474 TH2F* hCorrelationFMDFlatTr = (TH2F*)fDiagList->FindObject("hCorrelationFMDFlatTr");
475 TH2F* hCorrelationFMDRatioFlatTr = (TH2F*)fDiagList->FindObject("hCorrelationFMDRatioFlatTr");
476 TH1F* hTrVtxDistribution = (TH1F*)fDiagList->FindObject("hTrVtxDistribution");
477 TH1F* hTrEtaDistribution = (TH1F*)fDiagList->FindObject("hTrEtaDistribution");
478 TH1F* hTrEtaDistribution2 = (TH1F*)fDiagList->FindObject("hTrEtaDistribution2");
479 hCorrelationFMDSPD->Fill(nTrackLets,nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
480 TH1F* hFlatTracks = (TH1F*)fDiagList->FindObject("hFlatTracks");
481 hCorrelationFMD1SPD->Fill(nTrackLets,nHits[0][0]);
482 hCorrelationFMD2ISPD->Fill(nTrackLets,nHits[1][0]);
483 hCorrelationFMD2OSPD->Fill(nTrackLets,nHits[1][1]);
484 hCorrelationFMD3ISPD->Fill(nTrackLets,nHits[2][0]);
485 hCorrelationFMD3OSPD->Fill(nTrackLets,nHits[2][1]);
486 hCorrelationFMDSPDhits->Fill(testmult->GetNumberOfFiredChips(0),nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
487 hCorrelationSPDTracklets->Fill(testmult->GetNumberOfFiredChips(0),nTrackLets);
488
489 hTimeCorrelation->Fill(delta,nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
490 hCorrelationFMD23->Fill(nHits[1][0]+nHits[1][1],nHits[2][0]+nHits[2][1]);
491 hCorrelationFMD12->Fill(nHits[0][0],nHits[1][0]+nHits[1][1]);
492
493 // if(TMath::Abs(nHits[1]-nHits[2]) > 15 && (nHits[1]+nHits[2]) > 35) {fStatus = kFALSE; std::cout<<"difference : "<<TMath::Abs(nHits[1]-nHits[2])<<std::endl; return;}
494
495 // if(testmult->GetNumberOfFiredChips(0))
496 // if(testmult->GetNumberOfFiredChips(0) > 15 && ((Float_t)nTrackLets / (Float_t)testmult->GetNumberOfFiredChips(0)) < 0.4)
497 // {fStatus = kFALSE; std::cout<<nTrackLets<<" "<<testmult->GetNumberOfFiredChips(0)<<" "<<nHits[0]<<" "<<nHits[1]<<" "<<nHits[2]<<std::endl; return;}
498
499
500 Float_t diff23 = (Float_t)TMath::Abs(nHits[2][1] + nHits[2][0] - nHits[1][1] - nHits[1][0]);
501
502 Float_t diff13 = TMath::Abs(nHits[2][1] + nHits[2][0] - nHits[1][1] - nHits[1][0] - nHits[0][0]);
503 Float_t diff12 = TMath::Abs(nHits[1][0] - nHits[0][0]);
504
505 hCorrelationFMD1diff12->Fill(nHits[0][0], diff12);
506 hCorrelationFMD1diff13->Fill(nHits[0][0], diff13);
507 hCorrelationFMD2diff23->Fill(nHits[1][1], diff23);
508 hCorrelationFMD3diff23->Fill(nHits[2][1], diff23);
509
510 //
511 Float_t nTotalFMDhits = nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1] ;
512 Float_t radius = TMath::Sqrt(TMath::Power(vertex[0] + 0.03715,2) + TMath::Power(vertex[1] - 0.1659,2));
513
514 if(vertex[1] !=0 || vertex[1] !=0) {
515 hHitsRadius->Fill(nTotalFMDhits,radius);
516 hHitsX->Fill(nTotalFMDhits,vertex[0]);
517 hHitsY->Fill(nTotalFMDhits,vertex[1]); }
3bb122c7 518
da0805e2 519 hFMDHitDistribution->Fill(nTotalFMDhits);
520 hFMD1HitDistribution->Fill(nHits[0][0]);
521 hFMD2IHitDistribution->Fill(nHits[1][0]);
522 hFMD2OHitDistribution->Fill(nHits[1][1]);
523 hFMD3IHitDistribution->Fill(nHits[2][0]);
524 hFMD3OHitDistribution->Fill(nHits[2][0]);
525
526 // if(radius > 0.5) {fStatus = kFALSE; std::cout<<"FMD : "<<nTotalFMDhits<<std::endl; foutputESDFMD->Clear(); return;}
527 //if(TMath::Abs(vertex[1] - 0.1659) > 0.1 ) {fStatus = kFALSE; std::cout<<"FMD : "<<nTotalFMDhits<<std::endl; foutputESDFMD->Clear(); return;}
528
529 if(nTrackLets < pars->GetLowSPDLimit() || nTrackLets > pars->GetHighSPDLimit())
530 {fStatus = kFALSE; std::cout<<nTrackLets<<" "<<" "<<nHits[0][0]<<" "<<nHits[1][0]<<" "<<nHits[1][1]<<" "<<nHits[2][0]<<" "<<nHits[2][1]<<std::endl; return;}
531
532
533 AliESDtrack* track = 0;
534 Int_t ntracks = fESD->GetNumberOfTracks();
535 Float_t ngood =0, nbad = 0;
536 //std::cout<<" Primary vtx : "<<vertex[0]<<" "<<vertex[1]<<" "<<vertex[2]<<" "<<nTotalFMDhits<<std::endl;
537 Int_t nBgCandidates = 0;
538 Float_t nFlat = 0;
539 for(Int_t i=0;i<ntracks;i++) {
540 track = fESD->GetTrack(i);
541 //std::cout<<track->GetX()-vertex[0]<<" "<<track->GetY()-vertex[1]<<" "<<track->GetZ()-vertex[2]<<std::endl;
542 //std::cout<<track->GetX()<<" "<<track->GetY()<<" "<<track->GetZ()<<std::endl;
543 hTrVtxDistribution->Fill(track->GetZ());
544
545 if(TMath::Abs( track->GetZ()) > 50 && TMath::Abs(track->GetZ()) < 300) { // && TMath::Abs( track->GetY()) < 1)
546 nBgCandidates++;
547 hTrEtaDistribution->Fill(track->Eta());
548 }
549
550
551
552 if(TMath::Abs(track->GetX()-vertex[0]) > 0.3 || TMath::Abs(track->GetY()-vertex[1]) > 0.3 || TMath::Abs(track->GetZ()-vertex[2]) > 0.3) {
553 nbad++;
554 hTrEtaDistribution2->Fill(track->Eta()); }
555 else ngood++;
556
557 if(TMath::Abs(track->Pt()) < 0.1)
558 nFlat++;
559 }
560
561 Float_t ratioFlat = 0;
562 if(fESD->GetNumberOfTracks())
563 ratioFlat = nFlat/(Float_t)fESD->GetNumberOfTracks();
564 hCorrelationFMDFlatTr->Fill(nFlat,nTotalFMDhits);
565 hCorrelationFMDRatioFlatTr->Fill(ratioFlat,nTotalFMDhits);
566 hFlatTracks->Fill(nFlat);
b64db9b1 567
da0805e2 568 // std::cout<<fESD->GetT0zVertex()<<" "<<vertex[2]<<std::endl;
569 Float_t ratioBg = 0;
570 //if(fESD->GetNumberOfTracks() > 0)
571
572 if(fESD->GetNumberOfTracks() > 0)
573 ratioBg = (Float_t)nBgCandidates/(Float_t)fESD->GetNumberOfTracks();
574 hCorrelationFMDBgCand->Fill(nBgCandidates,nTotalFMDhits);
575 hCorrelationFMDBgCandRelative->Fill(ratioBg,nTotalFMDhits);
576
577
aa303f0c 578 // Float_t ratio = (nbad > 0 ? ngood / nbad : 0);
da0805e2 579
580 hCorrelationFMDGoodtracks->Fill(ngood,nTotalFMDhits);
581 hCorrelationFMDBadtracks->Fill(nbad,nTotalFMDhits);
582 hCorrelationGoodbadtracks->Fill(ngood,nbad);
583
7c3e5162 584 if(fStandalone) {
585 PostData(0, foutputESDFMD);
586 PostData(1, fEsdVertex);
587 PostData(2, fESD);
a38d1871 588 PostData(3, fDiagList);
7c3e5162 589 }
3bb122c7 590}
daedf077 591
8dc823cc 592//_____________________________________________________________________
593Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
bb8a464f 594 Float_t eta,
0210c863 595 Float_t prevE,
596 Float_t nextE,
bb8a464f 597 UShort_t det,
598 Char_t ring,
41bad769 599 UShort_t /*sec*/,
600 UShort_t /*strip*/) {
0210c863 601 //analyse and perform sharing on one strip
8dc823cc 602 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
6289b3e8 603
0210c863 604 Float_t mergedEnergy = 0;
cc066cb9 605 //Float_t nParticles = 0;
da0805e2 606 Float_t cutLow = 0.3;//0.15;
607
608 Float_t cutHigh = pars->GetMPV(det,ring,eta) - 2*pars->GetSigma(det,ring,eta);
aa303f0c 609
36f7d459 610 // if(ring == 'I')
611 // cutLow = 0.1;
cc066cb9 612
f58a4769 613 //cutLow = 0;
614 //AliFMDParameters* recopars = AliFMDParameters::Instance();
615 //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
da0805e2 616 //if(foutputESDFMD->GetUniqueID() == kFALSE ) {
617
618 if(mult > 12 || mult < cutLow)
619 {
620 // std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
621 fSharedThis = kFALSE;
622 fSharedPrev = kFALSE;
623 return 0;
624 }
625
cc066cb9 626
627
f58a4769 628
6289b3e8 629 // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
0210c863 630 Float_t totalE = mult;
bb8a464f 631
cc066cb9 632
cc066cb9 633 //std::cout<<det<<ring<<" "<<sec<<" "<<strip<<" "<<cutLow<<std::endl;
c78bc12b 634 if(fSharedThis) {
635 fSharedThis = kFALSE;
636 fSharedPrev = kTRUE;
8dc823cc 637 return 0.;
638 }
639
cc066cb9 640 /* if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
6289b3e8 641 fSharedThis = kFALSE;
642 fSharedPrev = kFALSE;
643 return 0;
cc066cb9 644 }*/
0210c863 645 if(mult<nextE && nextE>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
cc066cb9 646 {
647 fSharedThis = kFALSE;
648 fSharedPrev = kFALSE;
649 return 0;
650 }
da0805e2 651
bb8a464f 652
0210c863 653 if(prevE > cutLow && prevE < cutHigh && !fSharedPrev ) {
654 totalE += prevE;
8dc823cc 655 }
656
0210c863 657 if(nextE > cutLow && nextE < cutHigh ) {
658 totalE += nextE;
c78bc12b 659 fSharedThis = kTRUE;
8dc823cc 660 }
da0805e2 661 totalE = totalE*TMath::Cos(Eta2Theta(eta));
a38d1871 662 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
da0805e2 663 if(totalE > cutLow)
664 hEdist->Fill(totalE);
665
bb8a464f 666
0210c863 667 if(totalE > 0) {
cc066cb9 668
0210c863 669 mergedEnergy = totalE;
c78bc12b 670 fSharedPrev = kTRUE;
cc066cb9 671 // if(det == 1 && ring =='I')
c78bc12b 672 }
da0805e2 673 else{
674 fSharedThis = kFALSE;
675 fSharedPrev = kFALSE;
c78bc12b 676 }
da0805e2 677
0210c863 678 // mergedEnergy = mult;
8dc823cc 679
da0805e2 680
681 /* }
682else {
683 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
684 if(mult > cutLow)
685 fEtotal+=mult;
686 if(mult < cutLow) {
687 mergedEnergy = fEtotal;
688 fEtotal = 0;
689 hEdist->Fill(mergedEnergy);
690
691 }
692
693 }*/
694
0210c863 695 return mergedEnergy;
cc066cb9 696 //}
8dc823cc 697}
da0805e2 698//_____________________________________________________________________
699void AliFMDAnalysisTaskSharing::Terminate(Option_t* /* option*/) {
700
701 TH1F* hFMDHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistribution");
702 TH1F* hFMD1HitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD1");
703 TH1F* hFMD2IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2I");
704 TH1F* hFMD2OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2O");
705 TH1F* hFMD3IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3I");
706 TH1F* hFMD3OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3O");
707
708 TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
709 for(UShort_t det=1;det<=3;det++) {
710 Int_t nRings = (det==1 ? 1 : 2);
711 for (UShort_t ir = 0; ir < nRings; ir++) {
712 Char_t ring = (ir == 0 ? 'I' : 'O');
713
714 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
715 TH1F* hEdistAfter = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
716 if(hZvtx->GetEntries()) {
717 hEdist->Scale(1./(Float_t)hZvtx->GetEntries());
718 hEdistAfter->Scale(1./(Float_t)hZvtx->GetEntries());
719 }
720
721 }
722
723 }
724 TH1F* hFlatTracks = (TH1F*)fDiagList->FindObject("hFlatTracks");
725 TH1F* hTrVtxDistribution = (TH1F*)fDiagList->FindObject("hTrVtxDistribution");
726 TH1F* hTrEtaDistribution = (TH1F*)fDiagList->FindObject("hTrEtaDistribution");
727 TH1F* hTrEtaDistribution2 = (TH1F*)fDiagList->FindObject("hTrEtaDistribution2");
728 if(hZvtx->GetEntries()) {
729 hFMDHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
730 hFMD1HitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
731 hFMD2IHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
732 hFMD2OHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
733 hFMD3IHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
734 hFMD3OHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
735 hFlatTracks->Scale(1./(Float_t)hZvtx->GetEntries());
736 hTrVtxDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
737 hTrEtaDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
738 hTrEtaDistribution2->Scale(1./(Float_t)hZvtx->GetEntries());
739 }
740
059c7c6b 741 //TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
742 //TH1F* hEventsTr = (TH1F*)fDiagList->FindObject("hEventsTr");
daedf077 743 //hEventsVtx->Divide(hEventsTr);
da0805e2 744
745
746}
bb8a464f 747//_____________________________________________________________________
0210c863 748Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) const{
749 //convert the eta of a strip to a theta
bb8a464f 750 Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
751
752 if(eta < 0)
753 theta = theta-TMath::Pi();
754
cc066cb9 755 // std::cout<<"From eta2Theta: "<<theta<<" "<<eta<<std::endl;
bb8a464f 756 return theta;
757
758
759
760}
b64db9b1 761
762
763
bb8a464f 764//_____________________________________________________________________
b64db9b1 765void AliFMDAnalysisTaskSharing::ProcessPrimary() {
aa303f0c 766 //Get the unspoiled MC dN/deta before event cuts
b64db9b1 767 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
768 AliMCEvent* mcEvent = eventHandler->MCEvent();
769 if(!mcEvent)
770 return;
0a2f2742 771 fLastTrackByStrip.Reset(-1);
b64db9b1 772 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
773
774 AliMCParticle* particle = 0;
775
776 AliStack* stack = mcEvent->Stack();
777
778 TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
059c7c6b 779 TH1F* hPrimaryNSD = (TH1F*)fDiagList->FindObject("hMultvsEtaNSDNoCuts");
aa303f0c 780 TH1F* hEnergyOfParticles = (TH1F*)fDiagList->FindObject("hEnergyOfParticles");
b64db9b1 781 AliHeader* header = mcEvent->Header();
782 AliGenEventHeader* genHeader = header->GenEventHeader();
bb8a464f 783
850ad94d 784 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
0f6643af 785 Bool_t nsd = kTRUE;
850ad94d 786 if (!pythiaGenHeader) {
787 std::cout<<" no pythia header!"<<std::endl;
0f6643af 788 // return;
850ad94d 789 }
0f6643af 790 else {
bb8a464f 791
da0805e2 792
850ad94d 793 Int_t pythiaType = pythiaGenHeader->ProcessType();
0f6643af 794
059c7c6b 795 if(pythiaType==92||pythiaType==93)
796 nsd = kFALSE;
0f6643af 797 }
850ad94d 798
799 /*if(pythiaType==92||pythiaType==93){
800 std::cout<<"single diffractive"<<std::endl;
801 return;
802 }
803 if(pythiaType==94){
804 std::cout<<"double diffractive"<<std::endl;
805 return;
806 }
807 */
daedf077 808 //std::cout<<pythiaType<<std::endl;
850ad94d 809
bb8a464f 810
b64db9b1 811 TArrayF vertex;
812 genHeader->PrimaryVertex(vertex);
0a2f2742 813
b64db9b1 814 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
815 return;
816
817 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
818 Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
819 Int_t vertexBin = (Int_t)vertexBinDouble;
820
821 Bool_t firstTrack = kTRUE;
059c7c6b 822 Bool_t firstTrackNSD = kTRUE;
bb8a464f 823
b64db9b1 824 Int_t nTracks = stack->GetNprimary();
0a2f2742 825 if(pars->GetProcessHits())
826 nTracks = stack->GetNtrack();
b64db9b1 827 TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
059c7c6b 828 TH1F* nMCeventsNSD = (TH1F*)fDiagList->FindObject("nMCEventsNSDNoCuts");
b64db9b1 829 for(Int_t i = 0 ;i<nTracks;i++) {
7aad0c47 830 particle = (AliMCParticle*) mcEvent->GetTrack(i);
b64db9b1 831 if(!particle)
832 continue;
833
834 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
835 hPrimary->Fill(particle->Eta());
059c7c6b 836 if(nsd) {
837 hPrimaryNSD->Fill(particle->Eta());
838
839 if(firstTrackNSD) {
840 nMCeventsNSD->Fill(vertexBin);
841 firstTrackNSD = kFALSE;
842 }
843 }
b64db9b1 844
845 TH1F* hPrimVtxBin = (TH1F*)fDiagList->FindObject(Form("primmult_NoCuts_vtxbin%d",vertexBin));
846 hPrimVtxBin->Fill(particle->Eta());
847
848 if(firstTrack) {
849 nMCevents->Fill(vertexBin);
850 firstTrack = kFALSE;
851 }
852
853 }
0a2f2742 854 if(pars->GetProcessHits()) {
855
856 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
857
858 AliTrackReference* ref = particle->GetTrackReference(j);
859 UShort_t det,sec,strip;
860 Char_t ring;
861 if(ref->DetectorId() != AliTrackReference::kFMD)
862 continue;
aa303f0c 863 if(particle->Charge() != 0)
864 hEnergyOfParticles->Fill(particle->E());
865
0a2f2742 866 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
867 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
868 if(particle->Charge() != 0 && i != thisStripTrack ) {
869 //Double_t x,y,z;
870
871 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
872 TH1F* hHits = (TH1F*)fDiagList->FindObject(Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vertexBin));
873
874
875 hHits->Fill(eta);
876
877 Float_t nstrips = (ring =='O' ? 256 : 512);
878
879 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
880
881 if(strip >0)
882 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
883 if(strip < (nstrips - 1))
884 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
885
886 }
887
888
889 }
890
891
892 }
893
b64db9b1 894 }
895
896}
897
3bb122c7 898//_____________________________________________________________________
899//
900// EOF
901//