PWG2/FORWARD/analysis/AliFMDAnalysisTaskDndeta.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 "TH1F.h"
  10 #include "TH2F.h"
  11 #include "AliFMDAnalysisTaskDndeta.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliESDFMD.h"
  14 #include "AliESDEvent.h"
  15 #include "AliAODEvent.h"
  16 #include "AliAODHandler.h"
  17 #include "AliMCEventHandler.h"
  18 #include "AliStack.h"
  19 #include "AliLog.h"
  20 #include "AliESDVertex.h"
  21 #include "TMath.h"
  22 #include "AliFMDAnaParameters.h"
  23 //#include "AliFMDGeometry.h"
  24 #include "AliGenEventHeader.h"
  25 #include "AliGenPythiaEventHeader.h"
  26 #include "AliHeader.h"
  27 //#include "TDatabasePDG.h"
  28 //#include "TParticlePDG.h"
  29 #include "AliFMDStripIndex.h"
  30 #include "AliESDInputHandler.h"
  31 ClassImp(AliFMDAnalysisTaskDndeta)
  32
  33
  34 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
  35 : fDebug(0),
  36   fOutputList(0),
  37   fInputList(0),
  38   fVertexString(0x0),
  39   fNevents(),
  40   fNMCevents(),
  41   fStandalone(kTRUE),
  42   fLastTrackByStrip(0)
  43 {
  44   // Default constructor
  45   DefineInput (0, TList::Class());
  46   DefineOutput(0, TList::Class());
  47 }
  48 //_____________________________________________________________________
  49 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
  50     AliAnalysisTask(name, "Density"),
  51     fDebug(0),
  52     fOutputList(0),
  53     fInputList(0),
  54     fVertexString(0x0),
  55     fNevents(),
  56     fNMCevents(),
  57     fStandalone(kTRUE),
  58     fLastTrackByStrip(0)
  59 {
  60   fStandalone = SE;
  61   if(fStandalone) {
  62     DefineInput (0, TList::Class());
  63     DefineInput(1, TObjString::Class());
  64     DefineOutput(0, TList::Class());
  65
  66   }
  67 }
  68 //_____________________________________________________________________
  69 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
  70 {
  71   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  72
  73   if(!fOutputList)
  74     fOutputList = new TList();
  75   fOutputList->SetName("BackgroundCorrected");
  76
  77
  78   TH2F* hMult = 0;
  79   TH1F* hHits = 0;
  80   TH2F* hMultTrVtx = 0;
  81   TH1F* hPrimVertexBin = 0;
  82
  83
  84   TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
  85   TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
  86                             hBgTmp->GetNbinsX(),
  87                             hBgTmp->GetXaxis()->GetXmin(),
  88                             hBgTmp->GetXaxis()->GetXmax());
  89   hPrimary->Sumw2();
  90   fOutputList->Add(hPrimary);
  91   Int_t nVtxbins = pars->GetNvtxBins();
  92   TH2F* hBg = 0;
  93   for(Int_t i = 0; i< nVtxbins; i++) {
  94
  95     for(Int_t det =1; det<=3;det++)
  96       {
  97         Int_t nRings = (det==1 ? 1 : 2);
  98         for(Int_t ring = 0;ring<nRings;ring++)
  99           {
 100             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 101             Int_t  nSec     = (ring == 0 ? 20 : 40);
 102
 103
 104
 105             hBg = pars->GetBackgroundCorrection(det, ringChar, i);
 106             hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
 107                               hBg->GetNbinsX(),
 108                               hBg->GetXaxis()->GetXmin(),
 109                               hBg->GetXaxis()->GetXmax(),
 110                               nSec, 0, 2*TMath::Pi());
 111             hMultTrVtx  = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
 112                                    hBg->GetNbinsX(),
 113                                    hBg->GetXaxis()->GetXmin(),
 114                                    hBg->GetXaxis()->GetXmax(),
 115                                    nSec, 0, 2*TMath::Pi());
 116             hHits  = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
 117                               hBg->GetNbinsX(),
 118                               hBg->GetXaxis()->GetXmin(),
 119                               hBg->GetXaxis()->GetXmax());
 120             hHits->Sumw2();
 121             fOutputList->Add(hHits);
 122
 123             hMult->Sumw2();
 124             fOutputList->Add(hMult);
 125
 126             hMultTrVtx->Sumw2();
 127             fOutputList->Add(hMultTrVtx);
 128
 129           }
 130       }
 131   }
 132
 133   for(Int_t i = 0; i< nVtxbins; i++) {
 134
 135     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
 136                               Form("primmult_vtxbin%d",i),
 137                               hBgTmp->GetNbinsX(),
 138                               hBgTmp->GetXaxis()->GetXmin(),
 139                               hBgTmp->GetXaxis()->GetXmax());
 140     hPrimVertexBin->Sumw2();
 141     fOutputList->Add(hPrimVertexBin);
 142     //SPD part
 143     TH2F* hSPDMult  = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i),
 144                          hBgTmp->GetNbinsX(),
 145                          hBgTmp->GetXaxis()->GetXmin(),
 146                          hBgTmp->GetXaxis()->GetXmax(),
 147                          20, 0, 2*TMath::Pi());
 148     hSPDMult->Sumw2();
 149     fOutputList->Add(hSPDMult);
 150     TH2F* hSPDMultTrVtx  = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i),
 151                          hBgTmp->GetNbinsX(),
 152                          hBgTmp->GetXaxis()->GetXmin(),
 153                          hBgTmp->GetXaxis()->GetXmax(),
 154                          20, 0, 2*TMath::Pi());
 155     hSPDMultTrVtx->Sumw2();
 156     fOutputList->Add(hSPDMultTrVtx);
 157   }
 158
 159   fNevents.SetBins(nVtxbins,0,nVtxbins);
 160   fNevents.SetName("nEvents");
 161   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
 162   fNMCevents.SetName("nMCEvents");
 163   fOutputList->Add(&fNevents);
 164   fOutputList->Add(&fNMCevents);
 165
 166 }
 167 //_____________________________________________________________________
 168 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
 169 {
 170   if(fStandalone) {
 171     fInputList   = (TList*)GetInputData(0);
 172   }
 173 }
 174 //_____________________________________________________________________
 175 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
 176 {
 177   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 178   fVertexString = (TObjString*)fInputList->At(0);
 179
 180   Int_t vtxbin   = fVertexString->GetString().Atoi();
 181
 182   fNevents.Fill(vtxbin);
 183
 184   for(UShort_t det=1;det<=3;det++) {
 185     Int_t nRings = (det==1 ? 1 : 2);
 186     for (UShort_t ir = 0; ir < nRings; ir++) {
 187       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 188
 189       TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 190       TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 191
 192       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 193       TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 194
 195       hMultTotal->Add(hMultInput);
 196       hMultTotalTrVtx->Add(hMultInputTrVtx);
 197
 198     }
 199   }
 200
 201   //SPD code
 202   TH2F* hMultSPDTotal      = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin));
 203   TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin));
 204   TH2F* hMultSPDInput      = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
 205   TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
 206   hMultSPDTotal->Add(hMultSPDInput);
 207   hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx);
 208
 209
 210   if(pars->GetProcessPrimary())
 211     ProcessPrimary();
 212
 213   if(fStandalone) {
 214     PostData(0, fOutputList);
 215   }
 216
 217 }
 218 //_____________________________________________________________________
 219 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
 220
 221   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 222
 223   Int_t nVtxbins = pars->GetNvtxBins();
 224
 225   for(UShort_t det=1;det<=3;det++) {
 226     Int_t nRings = (det==1 ? 1 : 2);
 227     for (UShort_t ir = 0; ir < nRings; ir++) {
 228       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 229       for(Int_t i =0; i<nVtxbins; i++) {
 230
 231         TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
 232         //TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
 233         TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i));
 234
 235         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 236         TH1D* hMultProjTrVtx   = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 237         //fOutputList->Add(hMultTrVtx);
 238         fOutputList->Add(hMultProj);
 239         fOutputList->Add(hMultProjTrVtx);
 240       }
 241     }
 242   }
 243
 244   for(Int_t i =0; i<nVtxbins; i++) {
 245
 246     TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i));
 247     TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i));
 248
 249     TH1D* hMultProj   = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY());
 250     TH1D* hMultProjTrVtx   = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY());
 251
 252     fOutputList->Add(hMultProj);
 253     fOutputList->Add(hMultProjTrVtx);
 254
 255   }
 256
 257   std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" events"<<std::endl;
 258 }
 259 //_____________________________________________________________________
 260 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
 261
 262   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 263   AliMCEvent* mcEvent = eventHandler->MCEvent();
 264   if(!mcEvent)
 265     return;
 266
 267   fLastTrackByStrip.Reset(-1);
 268
 269   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 270
 271   AliMCParticle* particle = 0;
 272   AliStack* stack = mcEvent->Stack();
 273
 274   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
 275   AliHeader* header            = mcEvent->Header();
 276   AliGenEventHeader* genHeader = header->GenEventHeader();
 277
 278   TArrayF vertex;
 279   genHeader->PrimaryVertex(vertex);
 280   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 281     return;
 282   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 283   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 284   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 285
 286   Bool_t firstTrack = kTRUE;
 287
 288   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
 289   Int_t nTracks = stack->GetNprimary();
 290   if(pars->GetProcessHits())
 291     nTracks = stack->GetNtrack();
 292
 293   for(Int_t i = 0 ;i<nTracks;i++) {
 294     particle = (AliMCParticle*) mcEvent->GetTrack(i);
 295     if(!particle)
 296       continue;
 297
 298     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 299       hPrimary->Fill(particle->Eta());
 300
 301
 302       TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
 303       hPrimVtxBin->Fill(particle->Eta());
 304       if(firstTrack) {
 305         fNMCevents.Fill(vertexBin);
 306         firstTrack = kFALSE;
 307       }
 308
 309     }
 310     if(pars->GetProcessHits()) {
 311
 312       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 313
 314         AliTrackReference* ref = particle->GetTrackReference(j);
 315         UShort_t det,sec,strip;
 316         Char_t   ring;
 317         if(ref->DetectorId() != AliTrackReference::kFMD)
 318           continue;
 319         AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 320         Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 321         if(particle->Charge() != 0 && i != thisStripTrack ) {
 322           //Double_t x,y,z;
 323
 324           Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
 325           TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
 326
 327
 328           hHits->Fill(eta);
 329
 330           Float_t nstrips = (ring =='O' ? 256 : 512);
 331
 332           fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 333
 334           if(strip >0)
 335             fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 336           if(strip < (nstrips - 1))
 337             fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 338
 339         }
 340
 341
 342       }
 343
 344
 345     }
 346
 347
 348
 349   }
 350 }
 351 //_____________________________________________________________________
 352 //
 353 //
 354 // EOF