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