]> 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
fix codding conventions (Markus)
[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();
3bb122c7 308
da0805e2 309 Int_t delta = fESD->GetOrbitNumber() - fLastOrbit;
310 fLastOrbit = fESD->GetOrbitNumber();
311
312
b64db9b1 313 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
059c7c6b 314 Double_t vertex[3]={0,0,0};
f55d559b 315 Bool_t vtxStatus = pars->GetVertex(fESD,vertex);
bb8a464f 316 fEsdVertex->SetXYZ(vertex);
b64db9b1 317
059c7c6b 318 // std::cout<<vtxStatus<<" "<<vertex[0]<<" "<<vertex[1]<<" "<<vertex[2]<<std::endl;
319
b64db9b1 320 // Process primaries here to get true MC distribution
321 if(pars->GetProcessPrimary())
322 ProcessPrimary();
da0805e2 323 const AliMultiplicity* testmult = fESD->GetMultiplicity();
324 Int_t nTrackLets = testmult->GetNumberOfTracklets();
325 TH2F* hCorrelationClustersTracklets = (TH2F*)fDiagList->FindObject("hCorrelationClustersTracklets");
326 hCorrelationClustersTracklets->Fill(testmult->GetNumberOfSingleClusters(),nTrackLets);
b64db9b1 327
f55d559b 328
6289b3e8 329
059c7c6b 330 Bool_t isTriggered = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
331 Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
332
daedf077 333 Double_t delta2 = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
334 Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta2;
335 Int_t vtxbin = (Int_t)vertexBinDouble;
059c7c6b 336 TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
337 TH1F* hEventsTr = (TH1F*)fDiagList->FindObject("hEventsTr");
338 TH1F* hEventsNSD = (TH1F*)fDiagList->FindObject("hEventsNSD");
339 TH1F* hEventsNSDVtx = (TH1F*)fDiagList->FindObject("hEventsNSDVtx");
340
a82227ad 341 // if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
342 // vtxStatus = kFALSE;
343 // }
059c7c6b 344
345
daedf077 346 if(vtxStatus) hEventsVtx->Fill(vtxbin);
347 if(isTriggered) hEventsTr->Fill(vtxbin);
da0805e2 348
059c7c6b 349 if(vtxStatus && nsd) hEventsNSDVtx->Fill(vtxbin);
350 if(nsd) hEventsNSD->Fill(vtxbin);
351
da0805e2 352 if(!isTriggered || !vtxStatus ) {
bb8a464f 353 fStatus = kFALSE;
354 return;
355 }
356 else
357 fStatus = kTRUE;
b64db9b1 358
da0805e2 359 TH1F* hXvtx = (TH1F*)fDiagList->FindObject("hXvtx");
059c7c6b 360 if(vtxStatus ) hXvtx->Fill(vertex[0]);
da0805e2 361 TH1F* hYvtx = (TH1F*)fDiagList->FindObject("hYvtx");
059c7c6b 362 if(vtxStatus ) hYvtx->Fill(vertex[1]);
f55d559b 363 TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
364 hZvtx->Fill(vertex[2]);
da0805e2 365 //const AliMultiplicity* testmult = fESD->GetMultiplicity();
366 //std::cout<<vertex[2]<<std::endl;
367 //Int_t nTrackLets = testmult->GetNumberOfTracklets();
6289b3e8 368
a82227ad 369 if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
370 fStatus = kFALSE;
371 return;
372 }
da0805e2 373
6289b3e8 374 if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
375 else foutputESDFMD->SetUniqueID(kFALSE);
b64db9b1 376
3bb122c7 377 AliESDFMD* fmd = fESD->GetFMDData();
059c7c6b 378
3bb122c7 379 if (!fmd) return;
da0805e2 380 Int_t nHits[3][2] = {{0,0},{0,0},{0,0}};
381
3bb122c7 382 for(UShort_t det=1;det<=3;det++) {
383 Int_t nRings = (det==1 ? 1 : 2);
384 for (UShort_t ir = 0; ir < nRings; ir++) {
385 Char_t ring = (ir == 0 ? 'I' : 'O');
386 UShort_t nsec = (ir == 0 ? 20 : 40);
387 UShort_t nstr = (ir == 0 ? 512 : 256);
7c3e5162 388
a38d1871 389 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
7c3e5162 390
8dc7c4c2 391 for(UShort_t sec =0; sec < nsec; sec++) {
392 fSharedThis = kFALSE;
393 fSharedPrev = kFALSE;
0210c863 394
3bb122c7 395 for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162 396 foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
3bb122c7 397 Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
bb8a464f 398
8dc7c4c2 399 if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
7c3e5162 400
da0805e2 401
78f6f750 402 //Double_t eta = fmd->Eta(det,ring,sec,strip);
403 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
da0805e2 404
7c3e5162 405 hEdist->Fill(mult);
bb8a464f 406 if(fmd->IsAngleCorrected())
cc066cb9 407 mult = mult/TMath::Cos(Eta2Theta(eta));
0210c863 408 Float_t prevE = 0;
409 Float_t nextE = 0;
8dc823cc 410 if(strip != 0)
bb8a464f 411 if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
0210c863 412 prevE = fmd->Multiplicity(det,ring,sec,strip-1);
bb8a464f 413 if(fmd->IsAngleCorrected())
0210c863 414 prevE = prevE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
bb8a464f 415 }
8dc823cc 416 if(strip != nstr - 1)
bb8a464f 417 if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
0210c863 418 nextE = fmd->Multiplicity(det,ring,sec,strip+1);
bb8a464f 419 if(fmd->IsAngleCorrected())
0210c863 420 nextE = nextE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
bb8a464f 421 }
8dc823cc 422
0210c863 423 Float_t mergedEnergy = GetMultiplicityOfStrip(mult,eta,prevE,nextE,det,ring,sec,strip);
da0805e2 424 //if(mult> (pars->GetMPV(det,ring,eta) - pars->GetSigma(det,ring,eta)))
425 // mergedEnergy = mult;
426 //else mergedEnergy = 0;
427 if(mergedEnergy > 0.3 )
428 nHits[det-1][ir]++;
0210c863 429 foutputESDFMD->SetMultiplicity(det,ring,sec,strip,mergedEnergy);
bb8a464f 430 foutputESDFMD->SetEta(det,ring,sec,strip,eta);
7c3e5162 431
3bb122c7 432 }
433 }
434 }
435 }
da0805e2 436 //cluster cut
437 //if(testmult->GetNumberOfSingleClusters() > 15 + nTrackLets)
438 // {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;}
439
440 TH2F* hCorrelationFMD23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD23");
441 TH2F* hCorrelationFMD12 = (TH2F*)fDiagList->FindObject("hCorrelationFMD12");
442 TH2F* hCorrelationFMD2diff23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD2diff23");
443 TH2F* hCorrelationFMD3diff23 = (TH2F*)fDiagList->FindObject("hCorrelationFMD3diff23");
444 TH2F* hCorrelationFMD1diff13 = (TH2F*)fDiagList->FindObject("hCorrelationFMD1diff13");
445 TH2F* hCorrelationFMD1diff12 = (TH2F*)fDiagList->FindObject("hCorrelationFMD1diff12");
446
447 TH2F* hCorrelationFMDSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMDSPD");
448 TH2F* hCorrelationFMD1SPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD1SPD");
449 TH2F* hCorrelationFMD2ISPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD2ISPD");
450 TH2F* hCorrelationFMD2OSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD2OSPD");
451 TH2F* hCorrelationFMD3ISPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD3ISPD");
452 TH2F* hCorrelationFMD3OSPD = (TH2F*)fDiagList->FindObject("hCorrelationFMD3OSPD");
453
454 TH2F* hCorrelationFMDSPDhits = (TH2F*)fDiagList->FindObject("hCorrelationFMDSPDhits");
455
456 TH2F* hCorrelationSPDTracklets = (TH2F*)fDiagList->FindObject("hCorrelationSPDTracklets");
457 TH2F* hTimeCorrelation = (TH2F*)fDiagList->FindObject("hCorrelationTime");
458 TH2F* hHitsRadius = (TH2F*)fDiagList->FindObject("hCorrelationHitsRadius");
459 TH2F* hHitsX = (TH2F*)fDiagList->FindObject("hCorrelationHitsX");
460 TH2F* hHitsY = (TH2F*)fDiagList->FindObject("hCorrelationHitsY");
461 TH1F* hFMDHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistribution");
462 TH1F* hFMD1HitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD1");
463 TH1F* hFMD2IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2I");
464 TH1F* hFMD2OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2O");
465 TH1F* hFMD3IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3I");
466 TH1F* hFMD3OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3O");
467 TH1F* hCorrelationFMDGoodtracks = (TH1F*)fDiagList->FindObject("hCorrelationFMDGoodtracks");
468 TH1F* hCorrelationFMDBadtracks = (TH1F*)fDiagList->FindObject("hCorrelationFMDBadtracks");
469 TH1F* hCorrelationGoodbadtracks = (TH1F*)fDiagList->FindObject("hCorrelationGoodbadtracks");
470 TH2F* hCorrelationFMDBgCand = (TH2F*)fDiagList->FindObject("hCorrelationFMDBgCand");
471 TH2F* hCorrelationFMDBgCandRelative = (TH2F*)fDiagList->FindObject("hCorrelationFMDBgCandRelative");
472
473 TH2F* hCorrelationFMDFlatTr = (TH2F*)fDiagList->FindObject("hCorrelationFMDFlatTr");
474 TH2F* hCorrelationFMDRatioFlatTr = (TH2F*)fDiagList->FindObject("hCorrelationFMDRatioFlatTr");
475 TH1F* hTrVtxDistribution = (TH1F*)fDiagList->FindObject("hTrVtxDistribution");
476 TH1F* hTrEtaDistribution = (TH1F*)fDiagList->FindObject("hTrEtaDistribution");
477 TH1F* hTrEtaDistribution2 = (TH1F*)fDiagList->FindObject("hTrEtaDistribution2");
478 hCorrelationFMDSPD->Fill(nTrackLets,nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
479 TH1F* hFlatTracks = (TH1F*)fDiagList->FindObject("hFlatTracks");
480 hCorrelationFMD1SPD->Fill(nTrackLets,nHits[0][0]);
481 hCorrelationFMD2ISPD->Fill(nTrackLets,nHits[1][0]);
482 hCorrelationFMD2OSPD->Fill(nTrackLets,nHits[1][1]);
483 hCorrelationFMD3ISPD->Fill(nTrackLets,nHits[2][0]);
484 hCorrelationFMD3OSPD->Fill(nTrackLets,nHits[2][1]);
485 hCorrelationFMDSPDhits->Fill(testmult->GetNumberOfFiredChips(0),nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
486 hCorrelationSPDTracklets->Fill(testmult->GetNumberOfFiredChips(0),nTrackLets);
487
488 hTimeCorrelation->Fill(delta,nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1]);
489 hCorrelationFMD23->Fill(nHits[1][0]+nHits[1][1],nHits[2][0]+nHits[2][1]);
490 hCorrelationFMD12->Fill(nHits[0][0],nHits[1][0]+nHits[1][1]);
491
492 // 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;}
493
494 // if(testmult->GetNumberOfFiredChips(0))
495 // if(testmult->GetNumberOfFiredChips(0) > 15 && ((Float_t)nTrackLets / (Float_t)testmult->GetNumberOfFiredChips(0)) < 0.4)
496 // {fStatus = kFALSE; std::cout<<nTrackLets<<" "<<testmult->GetNumberOfFiredChips(0)<<" "<<nHits[0]<<" "<<nHits[1]<<" "<<nHits[2]<<std::endl; return;}
497
498
499 Float_t diff23 = (Float_t)TMath::Abs(nHits[2][1] + nHits[2][0] - nHits[1][1] - nHits[1][0]);
500
501 Float_t diff13 = TMath::Abs(nHits[2][1] + nHits[2][0] - nHits[1][1] - nHits[1][0] - nHits[0][0]);
502 Float_t diff12 = TMath::Abs(nHits[1][0] - nHits[0][0]);
503
504 hCorrelationFMD1diff12->Fill(nHits[0][0], diff12);
505 hCorrelationFMD1diff13->Fill(nHits[0][0], diff13);
506 hCorrelationFMD2diff23->Fill(nHits[1][1], diff23);
507 hCorrelationFMD3diff23->Fill(nHits[2][1], diff23);
508
509 //
510 Float_t nTotalFMDhits = nHits[0][0]+nHits[1][0]+nHits[1][1]+nHits[2][0]+nHits[2][1] ;
511 Float_t radius = TMath::Sqrt(TMath::Power(vertex[0] + 0.03715,2) + TMath::Power(vertex[1] - 0.1659,2));
512
513 if(vertex[1] !=0 || vertex[1] !=0) {
514 hHitsRadius->Fill(nTotalFMDhits,radius);
515 hHitsX->Fill(nTotalFMDhits,vertex[0]);
516 hHitsY->Fill(nTotalFMDhits,vertex[1]); }
3bb122c7 517
da0805e2 518 hFMDHitDistribution->Fill(nTotalFMDhits);
519 hFMD1HitDistribution->Fill(nHits[0][0]);
520 hFMD2IHitDistribution->Fill(nHits[1][0]);
521 hFMD2OHitDistribution->Fill(nHits[1][1]);
522 hFMD3IHitDistribution->Fill(nHits[2][0]);
523 hFMD3OHitDistribution->Fill(nHits[2][0]);
524
525 // if(radius > 0.5) {fStatus = kFALSE; std::cout<<"FMD : "<<nTotalFMDhits<<std::endl; foutputESDFMD->Clear(); return;}
526 //if(TMath::Abs(vertex[1] - 0.1659) > 0.1 ) {fStatus = kFALSE; std::cout<<"FMD : "<<nTotalFMDhits<<std::endl; foutputESDFMD->Clear(); return;}
527
528 if(nTrackLets < pars->GetLowSPDLimit() || nTrackLets > pars->GetHighSPDLimit())
529 {fStatus = kFALSE; std::cout<<nTrackLets<<" "<<" "<<nHits[0][0]<<" "<<nHits[1][0]<<" "<<nHits[1][1]<<" "<<nHits[2][0]<<" "<<nHits[2][1]<<std::endl; return;}
530
531
532 AliESDtrack* track = 0;
533 Int_t ntracks = fESD->GetNumberOfTracks();
534 Float_t ngood =0, nbad = 0;
535 //std::cout<<" Primary vtx : "<<vertex[0]<<" "<<vertex[1]<<" "<<vertex[2]<<" "<<nTotalFMDhits<<std::endl;
536 Int_t nBgCandidates = 0;
537 Float_t nFlat = 0;
538 for(Int_t i=0;i<ntracks;i++) {
539 track = fESD->GetTrack(i);
540 //std::cout<<track->GetX()-vertex[0]<<" "<<track->GetY()-vertex[1]<<" "<<track->GetZ()-vertex[2]<<std::endl;
541 //std::cout<<track->GetX()<<" "<<track->GetY()<<" "<<track->GetZ()<<std::endl;
542 hTrVtxDistribution->Fill(track->GetZ());
543
544 if(TMath::Abs( track->GetZ()) > 50 && TMath::Abs(track->GetZ()) < 300) { // && TMath::Abs( track->GetY()) < 1)
545 nBgCandidates++;
546 hTrEtaDistribution->Fill(track->Eta());
547 }
548
549
550
551 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) {
552 nbad++;
553 hTrEtaDistribution2->Fill(track->Eta()); }
554 else ngood++;
555
556 if(TMath::Abs(track->Pt()) < 0.1)
557 nFlat++;
558 }
559
560 Float_t ratioFlat = 0;
561 if(fESD->GetNumberOfTracks())
562 ratioFlat = nFlat/(Float_t)fESD->GetNumberOfTracks();
563 hCorrelationFMDFlatTr->Fill(nFlat,nTotalFMDhits);
564 hCorrelationFMDRatioFlatTr->Fill(ratioFlat,nTotalFMDhits);
565 hFlatTracks->Fill(nFlat);
b64db9b1 566
da0805e2 567 // std::cout<<fESD->GetT0zVertex()<<" "<<vertex[2]<<std::endl;
568 Float_t ratioBg = 0;
569 //if(fESD->GetNumberOfTracks() > 0)
570
571 if(fESD->GetNumberOfTracks() > 0)
572 ratioBg = (Float_t)nBgCandidates/(Float_t)fESD->GetNumberOfTracks();
573 hCorrelationFMDBgCand->Fill(nBgCandidates,nTotalFMDhits);
574 hCorrelationFMDBgCandRelative->Fill(ratioBg,nTotalFMDhits);
575
576
aa303f0c 577 // Float_t ratio = (nbad > 0 ? ngood / nbad : 0);
da0805e2 578
579 hCorrelationFMDGoodtracks->Fill(ngood,nTotalFMDhits);
580 hCorrelationFMDBadtracks->Fill(nbad,nTotalFMDhits);
581 hCorrelationGoodbadtracks->Fill(ngood,nbad);
582
7c3e5162 583 if(fStandalone) {
584 PostData(0, foutputESDFMD);
585 PostData(1, fEsdVertex);
586 PostData(2, fESD);
a38d1871 587 PostData(3, fDiagList);
7c3e5162 588 }
3bb122c7 589}
daedf077 590
8dc823cc 591//_____________________________________________________________________
592Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
bb8a464f 593 Float_t eta,
0210c863 594 Float_t prevE,
595 Float_t nextE,
bb8a464f 596 UShort_t det,
597 Char_t ring,
41bad769 598 UShort_t /*sec*/,
599 UShort_t /*strip*/) {
0210c863 600 //analyse and perform sharing on one strip
8dc823cc 601 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
6289b3e8 602
0210c863 603 Float_t mergedEnergy = 0;
cc066cb9 604 //Float_t nParticles = 0;
da0805e2 605 Float_t cutLow = 0.3;//0.15;
606
607 Float_t cutHigh = pars->GetMPV(det,ring,eta) - 2*pars->GetSigma(det,ring,eta);
aa303f0c 608
36f7d459 609 // if(ring == 'I')
610 // cutLow = 0.1;
cc066cb9 611
f58a4769 612 //cutLow = 0;
613 //AliFMDParameters* recopars = AliFMDParameters::Instance();
614 //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
da0805e2 615 //if(foutputESDFMD->GetUniqueID() == kFALSE ) {
616
617 if(mult > 12 || mult < cutLow)
618 {
619 // std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
620 fSharedThis = kFALSE;
621 fSharedPrev = kFALSE;
622 return 0;
623 }
624
cc066cb9 625
626
f58a4769 627
6289b3e8 628 // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
0210c863 629 Float_t totalE = mult;
bb8a464f 630
cc066cb9 631
cc066cb9 632 //std::cout<<det<<ring<<" "<<sec<<" "<<strip<<" "<<cutLow<<std::endl;
c78bc12b 633 if(fSharedThis) {
634 fSharedThis = kFALSE;
635 fSharedPrev = kTRUE;
8dc823cc 636 return 0.;
637 }
638
cc066cb9 639 /* if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
6289b3e8 640 fSharedThis = kFALSE;
641 fSharedPrev = kFALSE;
642 return 0;
cc066cb9 643 }*/
0210c863 644 if(mult<nextE && nextE>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
cc066cb9 645 {
646 fSharedThis = kFALSE;
647 fSharedPrev = kFALSE;
648 return 0;
649 }
da0805e2 650
bb8a464f 651
0210c863 652 if(prevE > cutLow && prevE < cutHigh && !fSharedPrev ) {
653 totalE += prevE;
8dc823cc 654 }
655
0210c863 656 if(nextE > cutLow && nextE < cutHigh ) {
657 totalE += nextE;
c78bc12b 658 fSharedThis = kTRUE;
8dc823cc 659 }
da0805e2 660 totalE = totalE*TMath::Cos(Eta2Theta(eta));
a38d1871 661 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
da0805e2 662 if(totalE > cutLow)
663 hEdist->Fill(totalE);
664
bb8a464f 665
0210c863 666 if(totalE > 0) {
cc066cb9 667
0210c863 668 mergedEnergy = totalE;
c78bc12b 669 fSharedPrev = kTRUE;
cc066cb9 670 // if(det == 1 && ring =='I')
c78bc12b 671 }
da0805e2 672 else{
673 fSharedThis = kFALSE;
674 fSharedPrev = kFALSE;
c78bc12b 675 }
da0805e2 676
0210c863 677 // mergedEnergy = mult;
8dc823cc 678
da0805e2 679
680 /* }
681else {
682 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
683 if(mult > cutLow)
684 fEtotal+=mult;
685 if(mult < cutLow) {
686 mergedEnergy = fEtotal;
687 fEtotal = 0;
688 hEdist->Fill(mergedEnergy);
689
690 }
691
692 }*/
693
0210c863 694 return mergedEnergy;
cc066cb9 695 //}
8dc823cc 696}
da0805e2 697//_____________________________________________________________________
698void AliFMDAnalysisTaskSharing::Terminate(Option_t* /* option*/) {
699
700 TH1F* hFMDHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistribution");
701 TH1F* hFMD1HitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD1");
702 TH1F* hFMD2IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2I");
703 TH1F* hFMD2OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD2O");
704 TH1F* hFMD3IHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3I");
705 TH1F* hFMD3OHitDistribution = (TH1F*)fDiagList->FindObject("hFMDHitDistributionFMD3O");
706
707 TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
708 for(UShort_t det=1;det<=3;det++) {
709 Int_t nRings = (det==1 ? 1 : 2);
710 for (UShort_t ir = 0; ir < nRings; ir++) {
711 Char_t ring = (ir == 0 ? 'I' : 'O');
712
713 TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
714 TH1F* hEdistAfter = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
715 if(hZvtx->GetEntries()) {
716 hEdist->Scale(1./(Float_t)hZvtx->GetEntries());
717 hEdistAfter->Scale(1./(Float_t)hZvtx->GetEntries());
718 }
719
720 }
721
722 }
723 TH1F* hFlatTracks = (TH1F*)fDiagList->FindObject("hFlatTracks");
724 TH1F* hTrVtxDistribution = (TH1F*)fDiagList->FindObject("hTrVtxDistribution");
725 TH1F* hTrEtaDistribution = (TH1F*)fDiagList->FindObject("hTrEtaDistribution");
726 TH1F* hTrEtaDistribution2 = (TH1F*)fDiagList->FindObject("hTrEtaDistribution2");
727 if(hZvtx->GetEntries()) {
728 hFMDHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
729 hFMD1HitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
730 hFMD2IHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
731 hFMD2OHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
732 hFMD3IHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
733 hFMD3OHitDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
734 hFlatTracks->Scale(1./(Float_t)hZvtx->GetEntries());
735 hTrVtxDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
736 hTrEtaDistribution->Scale(1./(Float_t)hZvtx->GetEntries());
737 hTrEtaDistribution2->Scale(1./(Float_t)hZvtx->GetEntries());
738 }
739
059c7c6b 740 //TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
741 //TH1F* hEventsTr = (TH1F*)fDiagList->FindObject("hEventsTr");
daedf077 742 //hEventsVtx->Divide(hEventsTr);
da0805e2 743
744
745}
bb8a464f 746//_____________________________________________________________________
0210c863 747Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) const{
748 //convert the eta of a strip to a theta
bb8a464f 749 Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
750
751 if(eta < 0)
752 theta = theta-TMath::Pi();
753
cc066cb9 754 // std::cout<<"From eta2Theta: "<<theta<<" "<<eta<<std::endl;
bb8a464f 755 return theta;
756
757
758
759}
b64db9b1 760
761
762
bb8a464f 763//_____________________________________________________________________
b64db9b1 764void AliFMDAnalysisTaskSharing::ProcessPrimary() {
aa303f0c 765 //Get the unspoiled MC dN/deta before event cuts
b64db9b1 766 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
767 AliMCEvent* mcEvent = eventHandler->MCEvent();
768 if(!mcEvent)
769 return;
0a2f2742 770 fLastTrackByStrip.Reset(-1);
b64db9b1 771 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
772
773 AliMCParticle* particle = 0;
774
775 AliStack* stack = mcEvent->Stack();
776
777 TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
059c7c6b 778 TH1F* hPrimaryNSD = (TH1F*)fDiagList->FindObject("hMultvsEtaNSDNoCuts");
aa303f0c 779 TH1F* hEnergyOfParticles = (TH1F*)fDiagList->FindObject("hEnergyOfParticles");
b64db9b1 780 AliHeader* header = mcEvent->Header();
781 AliGenEventHeader* genHeader = header->GenEventHeader();
bb8a464f 782
850ad94d 783 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
784
785 if (!pythiaGenHeader) {
786 std::cout<<" no pythia header!"<<std::endl;
787 return;
788 }
bb8a464f 789
da0805e2 790
850ad94d 791 Int_t pythiaType = pythiaGenHeader->ProcessType();
059c7c6b 792 Bool_t nsd = kTRUE;
793 if(pythiaType==92||pythiaType==93)
794 nsd = kFALSE;
795
850ad94d 796
797 /*if(pythiaType==92||pythiaType==93){
798 std::cout<<"single diffractive"<<std::endl;
799 return;
800 }
801 if(pythiaType==94){
802 std::cout<<"double diffractive"<<std::endl;
803 return;
804 }
805 */
daedf077 806 //std::cout<<pythiaType<<std::endl;
850ad94d 807
bb8a464f 808
b64db9b1 809 TArrayF vertex;
810 genHeader->PrimaryVertex(vertex);
0a2f2742 811
b64db9b1 812 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
813 return;
814
815 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
816 Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
817 Int_t vertexBin = (Int_t)vertexBinDouble;
818
819 Bool_t firstTrack = kTRUE;
059c7c6b 820 Bool_t firstTrackNSD = kTRUE;
bb8a464f 821
b64db9b1 822 Int_t nTracks = stack->GetNprimary();
0a2f2742 823 if(pars->GetProcessHits())
824 nTracks = stack->GetNtrack();
b64db9b1 825 TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
059c7c6b 826 TH1F* nMCeventsNSD = (TH1F*)fDiagList->FindObject("nMCEventsNSDNoCuts");
b64db9b1 827 for(Int_t i = 0 ;i<nTracks;i++) {
7aad0c47 828 particle = (AliMCParticle*) mcEvent->GetTrack(i);
b64db9b1 829 if(!particle)
830 continue;
831
832 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
833 hPrimary->Fill(particle->Eta());
059c7c6b 834 if(nsd) {
835 hPrimaryNSD->Fill(particle->Eta());
836
837 if(firstTrackNSD) {
838 nMCeventsNSD->Fill(vertexBin);
839 firstTrackNSD = kFALSE;
840 }
841 }
b64db9b1 842
843 TH1F* hPrimVtxBin = (TH1F*)fDiagList->FindObject(Form("primmult_NoCuts_vtxbin%d",vertexBin));
844 hPrimVtxBin->Fill(particle->Eta());
845
846 if(firstTrack) {
847 nMCevents->Fill(vertexBin);
848 firstTrack = kFALSE;
849 }
850
851 }
0a2f2742 852 if(pars->GetProcessHits()) {
853
854 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
855
856 AliTrackReference* ref = particle->GetTrackReference(j);
857 UShort_t det,sec,strip;
858 Char_t ring;
859 if(ref->DetectorId() != AliTrackReference::kFMD)
860 continue;
aa303f0c 861 if(particle->Charge() != 0)
862 hEnergyOfParticles->Fill(particle->E());
863
0a2f2742 864 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
865 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
866 if(particle->Charge() != 0 && i != thisStripTrack ) {
867 //Double_t x,y,z;
868
869 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
870 TH1F* hHits = (TH1F*)fDiagList->FindObject(Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vertexBin));
871
872
873 hHits->Fill(eta);
874
875 Float_t nstrips = (ring =='O' ? 256 : 512);
876
877 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
878
879 if(strip >0)
880 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
881 if(strip < (nstrips - 1))
882 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
883
884 }
885
886
887 }
888
889
890 }
891
b64db9b1 892 }
893
894}
895
3bb122c7 896//_____________________________________________________________________
897//
898// EOF
899//