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