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