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