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