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