Various upgrades and bug fixes etc
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskSharing.cxx
1  
2 #include <TROOT.h>
3 #include <TSystem.h>
4 #include <TInterpreter.h>
5 #include <TChain.h>
6 #include <TFile.h>
7 #include <TList.h>
8 #include <iostream>
9 #include <TMath.h>
10 //#include "AliFMDDebug.h"
11 #include "AliFMDAnalysisTaskSharing.h"
12 #include "AliAnalysisManager.h"
13 #include "AliESDFMD.h"
14 //#include "AliFMDGeometry.h"
15 #include "AliMCEventHandler.h"
16 #include "AliMCEvent.h"
17 #include "AliStack.h"
18 #include "AliESDVertex.h"
19 #include "AliMultiplicity.h"
20 #include "AliFMDAnaParameters.h"
21 #include "TH1F.h"
22 #include "TObjString.h"
23 //#include "/home/canute/ALICE/AliRoot/PWG0/AliPWG0Helper.h"
24 //#include "AliFMDParameters.h"
25 #include "AliGenEventHeader.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliHeader.h"
28 #include "AliStack.h"
29 #include "AliMCParticle.h"
30 #include "AliFMDStripIndex.h"
31 #include "AliESDVZERO.h"
32 #include "AliESDtrack.h"
33 #include "AliGenDPMjetEventHeader.h"
34
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
42 ClassImp(AliFMDAnalysisTaskSharing)
43
44 //_____________________________________________________________________
45 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
46 : fDebug(0),
47   fESD(0x0),
48   foutputESDFMD(),
49   fSharedThis(kFALSE),
50   fSharedPrev(kFALSE),
51   fDiagList(0),
52   fStandalone(kTRUE),
53   fEsdVertex(0),
54   fStatus(kTRUE),
55   fLastTrackByStrip(0),
56   fLastOrbit(0)
57 {
58   // Default constructor
59   DefineInput (0, AliESDEvent::Class());
60   DefineOutput(0, AliESDFMD::Class());
61   DefineOutput(1, AliESDVertex::Class());
62   DefineOutput(2, AliESDEvent::Class());
63   DefineOutput(3, TList::Class());
64 }
65 //_____________________________________________________________________
66 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
67     AliAnalysisTask(name, "AnalysisTaskFMD"),
68     fDebug(0),
69     fESD(0x0),
70     foutputESDFMD(),
71     fSharedThis(kFALSE),
72     fSharedPrev(kFALSE),
73     fDiagList(0),
74     fStandalone(kTRUE),
75     fEsdVertex(0),
76     fStatus(kTRUE),
77     fLastTrackByStrip(0),
78     fLastOrbit(0)
79 {
80   // named constructor
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   }
89 }
90 //_____________________________________________________________________
91 Float_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 }
101 //_____________________________________________________________________
102 Float_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
114 //_____________________________________________________________________
115 void AliFMDAnalysisTaskSharing::CreateOutputObjects()
116 {
117   // Create the output objects
118   if(!foutputESDFMD)
119     foutputESDFMD = new AliESDFMD();
120   
121   if(!fEsdVertex)
122     fEsdVertex    = new AliESDVertex();
123   //Diagnostics
124   if(!fDiagList)
125     fDiagList = new TList();
126   
127   fDiagList->SetName("Sharing diagnostics");
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);
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   
145   TH1F* hXvtx = new TH1F("hXvtx","x vertex distribution",100,-2,2);
146   TH1F* hYvtx = new TH1F("hYvtx","y vertex distribution",100,-2,2);
147   TH1F* hZvtx = new TH1F("hZvtx","z vertex distribution",4*pars->GetNvtxBins(),-4*pars->GetVtxCutZ(),4*pars->GetVtxCutZ());
148   
149   fDiagList->Add(hXvtx);
150   fDiagList->Add(hYvtx);
151   fDiagList->Add(hZvtx);
152   
153   TH1F* hPrimVertexBin = 0;
154   TH1F* hHits = 0;
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   
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),
174                                      1000,0,25);
175       TH1F* hEdistAfter  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
176                                      Form("Edist_after_sharing_FMD%d%c", det, ringChar),
177                                      1000,0,25);
178       
179       
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);
183       fDiagList->Add(hEdist);
184       fDiagList->Add(hEdistAfter);
185       //fDiagList->Add(hNstripsHit);
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);
194
195       }
196       
197     }
198   }
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);
226 fDiagList->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);
264  
265  TH1F* hEnergyOfParticles = new TH1F("hEnergyOfParticles","hEnergyOfParticles",1000000,0,10);
266  fDiagList->Add(hEnergyOfParticles);
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);
277   TH1F*  nMCevents = new TH1F("nMCEventsNoCuts","nMCEventsNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
278   TH1F*  nMCeventsNSD= new TH1F("nMCEventsNSDNoCuts","nMCEventsNSDNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
279   fDiagList->Add(nMCevents);
280   fDiagList->Add(nMCeventsNSD);
281   
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());
286   fDiagList->Add(hEventsVtx);
287   fDiagList->Add(hEventsTr);
288   fDiagList->Add(hEventsNSD);
289   fDiagList->Add(hEventsNSDVtx);
290   
291 }
292 //_____________________________________________________________________
293 void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
294 {
295   // connect the input data
296   if(fStandalone)
297     fESD = (AliESDEvent*)GetInputData(0);
298 }
299 //_____________________________________________________________________
300 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
301 {
302   //perform analysis on one event
303   AliESD* old = fESD->GetAliESDOld();
304   if (old) {
305     fESD->CopyFromOldESD();
306   }
307   
308   foutputESDFMD->Clear();
309
310     
311   Int_t delta = fESD->GetOrbitNumber() - fLastOrbit;
312   fLastOrbit = fESD->GetOrbitNumber();
313   
314   
315   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
316   Double_t vertex[3]={0,0,0};
317   Bool_t vtxStatus = pars->GetVertex(fESD,vertex);
318   fEsdVertex->SetXYZ(vertex);
319   
320   // std::cout<<vtxStatus<<"  "<<vertex[0]<<"   "<<vertex[1]<<"   "<<vertex[2]<<std::endl;
321   
322   // Process primaries here to get true MC distribution
323   if(pars->GetProcessPrimary())
324     ProcessPrimary();
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);  
329   
330   
331   
332   Bool_t isTriggered = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
333   Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
334
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;
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   
343   // if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
344   //  vtxStatus = kFALSE;
345   // }
346   
347   
348   if(vtxStatus) hEventsVtx->Fill(vtxbin);
349   if(isTriggered) hEventsTr->Fill(vtxbin);
350   
351   if(vtxStatus && nsd) hEventsNSDVtx->Fill(vtxbin);
352   if(nsd) hEventsNSD->Fill(vtxbin);
353   
354   if(!isTriggered || !vtxStatus ) {
355     fStatus = kFALSE;
356     return;
357   }
358   else
359     fStatus = kTRUE;
360   
361   TH1F* hXvtx = (TH1F*)fDiagList->FindObject("hXvtx");
362   if(vtxStatus ) hXvtx->Fill(vertex[0]);
363   TH1F* hYvtx = (TH1F*)fDiagList->FindObject("hYvtx");
364   if(vtxStatus ) hYvtx->Fill(vertex[1]);
365   TH1F* hZvtx = (TH1F*)fDiagList->FindObject("hZvtx");
366   hZvtx->Fill(vertex[2]);
367   //const AliMultiplicity* testmult = fESD->GetMultiplicity();
368   //std::cout<<vertex[2]<<std::endl;
369   //Int_t nTrackLets = testmult->GetNumberOfTracklets();
370   
371   if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
372     fStatus = kFALSE;
373     return;
374   }
375     
376   if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
377   else foutputESDFMD->SetUniqueID(kFALSE);
378   
379   AliESDFMD* fmd = fESD->GetFMDData();
380    
381   if (!fmd) return;
382   Int_t nHits[3][2] = {{0,0},{0,0},{0,0}};
383
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);
390       
391       TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
392       
393       for(UShort_t sec =0; sec < nsec;  sec++) {
394         fSharedThis      = kFALSE;
395         fSharedPrev      = kFALSE;
396         
397         for(UShort_t strip = 0; strip < nstr; strip++) {
398           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
399           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
400           
401           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
402           
403          
404           //Double_t eta = fmd->Eta(det,ring,sec,strip);
405           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
406                   
407           hEdist->Fill(mult);
408           if(fmd->IsAngleCorrected())
409             mult = mult/TMath::Cos(Eta2Theta(eta));
410           Float_t prevE = 0;
411           Float_t nextE = 0;
412           if(strip != 0)
413             if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
414               prevE = fmd->Multiplicity(det,ring,sec,strip-1);
415               if(fmd->IsAngleCorrected())
416                 prevE = prevE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
417             }
418           if(strip != nstr - 1)
419             if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
420               nextE = fmd->Multiplicity(det,ring,sec,strip+1);
421               if(fmd->IsAngleCorrected())
422                 nextE = nextE/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
423             }
424           
425           Float_t mergedEnergy = GetMultiplicityOfStrip(mult,eta,prevE,nextE,det,ring,sec,strip);
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]++;
431           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,mergedEnergy);
432           foutputESDFMD->SetEta(det,ring,sec,strip,eta);
433           
434         }
435       }
436     }
437   }
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]); }
519   
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);
568    
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   
579   // Float_t ratio =  (nbad > 0 ? ngood / nbad  : 0);
580   
581   hCorrelationFMDGoodtracks->Fill(ngood,nTotalFMDhits);
582   hCorrelationFMDBadtracks->Fill(nbad,nTotalFMDhits);
583   hCorrelationGoodbadtracks->Fill(ngood,nbad);
584     
585   if(fStandalone) {
586     PostData(0, foutputESDFMD); 
587     PostData(1, fEsdVertex); 
588     PostData(2, fESD); 
589     PostData(3, fDiagList); 
590   }
591 }
592
593 //_____________________________________________________________________
594 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
595                                                           Float_t eta,
596                                                           Float_t prevE,
597                                                           Float_t nextE,
598                                                           UShort_t   det,
599                                                           Char_t  ring,
600                                                           UShort_t /*sec*/,
601                                                           UShort_t /*strip*/) {
602   //analyse and perform sharing on one strip
603   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
604  
605   Float_t mergedEnergy = 0;
606   //Float_t nParticles = 0;
607   Float_t cutLow  = 0.3;//0.15;
608   
609   Float_t cutHigh = pars->GetMPV(det,ring,eta) - 2*pars->GetSigma(det,ring,eta);
610
611   // if(ring == 'I')
612   //  cutLow = 0.1;
613   
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());
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   
627   
628   
629   
630   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
631   Float_t totalE  = mult;
632   
633
634     //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
635   if(fSharedThis) {
636     fSharedThis      = kFALSE;
637     fSharedPrev      = kTRUE;
638     return 0.;
639   }
640   
641   /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
642     fSharedThis      = kFALSE;
643     fSharedPrev      = kFALSE;
644     return 0;
645     }*/
646   if(mult<nextE && nextE>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
647     {
648       fSharedThis      = kFALSE;
649       fSharedPrev      = kFALSE;
650       return 0;
651     }
652   
653   
654   if(prevE > cutLow && prevE < cutHigh && !fSharedPrev ) {
655     totalE += prevE;
656   }
657   
658   if(nextE > cutLow && nextE < cutHigh ) {
659     totalE += nextE;
660     fSharedThis      = kTRUE;
661   }
662   totalE = totalE*TMath::Cos(Eta2Theta(eta));
663   TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
664   if(totalE > cutLow)
665     hEdist->Fill(totalE);
666   
667   
668   if(totalE > 0) {
669     
670     mergedEnergy = totalE;
671     fSharedPrev      = kTRUE;
672     // if(det == 1 && ring =='I')
673   }
674   else{
675       fSharedThis      = kFALSE;
676       fSharedPrev      = kFALSE;
677   }
678
679   // mergedEnergy = mult;
680   
681
682   /*   }
683 else {
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   
696   return mergedEnergy; 
697   //}  
698 }
699 //_____________________________________________________________________
700 void 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   
742   //TH1F* hEventsVtx = (TH1F*)fDiagList->FindObject("hEventsVtx");
743   //TH1F* hEventsTr  = (TH1F*)fDiagList->FindObject("hEventsTr");
744   //hEventsVtx->Divide(hEventsTr);
745   
746   
747 }
748 //_____________________________________________________________________
749 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) const{
750   //convert the eta of a strip to a theta
751   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
752   
753   if(eta < 0)
754     theta = theta-TMath::Pi();
755   
756   //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
757   return theta;
758   
759
760
761 }
762
763
764
765 //_____________________________________________________________________
766 void AliFMDAnalysisTaskSharing::ProcessPrimary() {
767   //Get the unspoiled MC dN/deta before event cuts
768   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
769   AliMCEvent* mcEvent = eventHandler->MCEvent();
770   if(!mcEvent)
771     return;
772   fLastTrackByStrip.Reset(-1);
773   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
774   
775   AliMCParticle* particle = 0;
776   
777   AliStack* stack = mcEvent->Stack();
778   
779   TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
780   TH1F* hPrimaryNSD = (TH1F*)fDiagList->FindObject("hMultvsEtaNSDNoCuts");
781   TH1F* hEnergyOfParticles = (TH1F*)fDiagList->FindObject("hEnergyOfParticles");
782   AliHeader* header            = mcEvent->Header();
783   AliGenEventHeader* genHeader = header->GenEventHeader();
784   
785   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
786   AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(header->GenEventHeader());
787   
788  
789   Bool_t nsd = kTRUE;
790   
791   
792   if (!pythiaGenHeader && !dpmHeader) {
793     std::cout<<" no pythia or dpm header!"<<std::endl;
794   }
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     }
809   }
810
811   TArrayF vertex;
812   genHeader->PrimaryVertex(vertex);
813   
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;
822   Bool_t firstTrackNSD = kTRUE;
823   
824   Int_t nTracks = stack->GetNprimary();
825   if(pars->GetProcessHits())
826     nTracks = stack->GetNtrack();
827   TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
828   TH1F* nMCeventsNSD = (TH1F*)fDiagList->FindObject("nMCEventsNSDNoCuts");
829   for(Int_t i = 0 ;i<nTracks;i++) {
830     particle = (AliMCParticle*) mcEvent->GetTrack(i);
831     if(!particle)
832       continue;
833     
834     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
835       hPrimary->Fill(particle->Eta());
836       if(nsd) {
837         hPrimaryNSD->Fill(particle->Eta());
838         
839         if(firstTrackNSD) {
840           nMCeventsNSD->Fill(vertexBin);
841           firstTrackNSD = kFALSE;
842         }
843       }
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     }
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;
863         if(particle->Charge() != 0)
864           hEnergyOfParticles->Fill(particle->E());
865         
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     
894   }
895
896 }
897
898 //_____________________________________________________________________
899 //
900 // EOF
901 //