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