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 "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   TH2F* hBg = 0;
  90   for(Int_t i = 0; i< nVtxbins; i++) {
  91
  92     for(Int_t det =1; det<=3;det++)
  93       {
  94         Int_t nRings = (det==1 ? 1 : 2);
  95         for(Int_t ring = 0;ring<nRings;ring++)
  96           {
  97             Char_t ringChar = (ring == 0 ? 'I' : 'O');
  98             Int_t  nSec     = (ring == 0 ? 20 : 40);
  99
 100
 101
 102             hBg = pars->GetBackgroundCorrection(det, ringChar, i);
 103             hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
 104                               hBg->GetNbinsX(),
 105                               hBg->GetXaxis()->GetXmin(),
 106                               hBg->GetXaxis()->GetXmax(),
 107                               nSec, 0, 2*TMath::Pi());
 108
 109             hHits  = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
 110                               hBg->GetNbinsX(),
 111                               hBg->GetXaxis()->GetXmin(),
 112                               hBg->GetXaxis()->GetXmax());
 113             hHits->Sumw2();
 114             fOutputList->Add(hHits);
 115
 116             hMult->Sumw2();
 117             fOutputList->Add(hMult);
 118
 119           }
 120       }
 121   }
 122
 123   for(Int_t i = 0; i< nVtxbins; i++) {
 124
 125     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
 126                               Form("primmult_vtxbin%d",i),
 127                               hBgTmp->GetNbinsX(),
 128                               hBgTmp->GetXaxis()->GetXmin(),
 129                               hBgTmp->GetXaxis()->GetXmax());
 130     hPrimVertexBin->Sumw2();
 131     fOutputList->Add(hPrimVertexBin);
 132
 133   }
 134
 135   fNevents.SetBins(nVtxbins,0,nVtxbins);
 136   fNevents.SetName("nEvents");
 137   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
 138   fNMCevents.SetName("nMCEvents");
 139   fOutputList->Add(&fNevents);
 140   fOutputList->Add(&fNMCevents);
 141
 142 }
 143 //_____________________________________________________________________
 144 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
 145 {
 146   if(fStandalone) {
 147     fInputList   = (TList*)GetInputData(0);
 148   }
 149 }
 150 //_____________________________________________________________________
 151 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
 152 {
 153   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 154   fVertexString = (TObjString*)fInputList->At(0);
 155
 156   Int_t vtxbin   = fVertexString->GetString().Atoi();
 157
 158   fNevents.Fill(vtxbin);
 159
 160   for(UShort_t det=1;det<=3;det++) {
 161     Int_t nRings = (det==1 ? 1 : 2);
 162     for (UShort_t ir = 0; ir < nRings; ir++) {
 163       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 164
 165       TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 166
 167
 168       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
 169
 170       hMultTotal->Add(hMultInput);
 171
 172     }
 173   }
 174
 175   if(pars->GetProcessPrimary())
 176     ProcessPrimary();
 177
 178   if(fStandalone) {
 179     PostData(0, fOutputList);
 180   }
 181
 182 }
 183 //_____________________________________________________________________
 184 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
 185
 186
 187   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 188
 189   Int_t nVtxbins = pars->GetNvtxBins();
 190
 191   for(UShort_t det=1;det<=3;det++) {
 192     Int_t nRings = (det==1 ? 1 : 2);
 193     for (UShort_t ir = 0; ir < nRings; ir++) {
 194       Char_t ringChar = (ir == 0 ? 'I' : 'O');
 195       for(Int_t i =0; i<nVtxbins; i++) {
 196
 197         TH1F* hSharingEff = pars->GetSharingEfficiency(det,ringChar,i);
 198         TH1F* hSharingEffTrVtx = pars->GetSharingEfficiencyTrVtx(det,ringChar,i);
 199         TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
 200         TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
 201
 202         for(Int_t nx=1; nx<hMultTotal->GetNbinsX(); nx++) {
 203           Float_t correction = hSharingEff->GetBinContent(nx);
 204           Float_t correctionTrVtx = hSharingEffTrVtx->GetBinContent(nx);
 205           for(Int_t ny=1; ny<hMultTotal->GetNbinsY(); ny++) {
 206
 207             if(correction != 0){
 208               hMultTotal->SetBinContent(nx,ny,hMultTotal->GetBinContent(nx,ny)/correction);
 209               Float_t error = TMath::Sqrt(TMath::Power(hMultTotal->GetBinError(nx,ny),2) + TMath::Power(hMultTotal->GetBinContent(nx,ny)*hSharingEff->GetBinError(nx),2)) / correction;
 210               hMultTotal->SetBinError(nx,ny,error);
 211             }
 212             if(correctionTrVtx != 0){
 213               hMultTrVtx->SetBinContent(nx,ny,hMultTrVtx->GetBinContent(nx,ny)/correctionTrVtx);
 214               Float_t error = TMath::Sqrt(TMath::Power(hMultTrVtx->GetBinError(nx,ny),2) + TMath::Power(hMultTrVtx->GetBinContent(nx,ny)*hSharingEffTrVtx->GetBinError(nx),2)) / correctionTrVtx;
 215               hMultTrVtx->SetBinError(nx,ny,error);
 216             }
 217           }
 218
 219         }
 220
 221         //hMultTotal->Divide(hSharingEff);
 222
 223         hMultTotal->Scale(1/pars->GetEventSelectionEfficiency(i));
 224
 225         //hMultTrVtx->Divide(hSharingEffTrVtx);
 226
 227         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 228         TH1D* hMultProjTrVtx   = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
 229         fOutputList->Add(hMultTrVtx);
 230         fOutputList->Add(hMultProj);
 231         fOutputList->Add(hMultProjTrVtx);
 232       }
 233     }
 234   }
 235
 236 }
 237 //_____________________________________________________________________
 238 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
 239
 240   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 241   AliMCEvent* mcEvent = eventHandler->MCEvent();
 242   if(!mcEvent)
 243     return;
 244
 245   fLastTrackByStrip.Reset(-1);
 246
 247   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 248
 249   AliMCParticle* particle = 0;
 250   AliStack* stack = mcEvent->Stack();
 251
 252   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
 253   AliHeader* header            = mcEvent->Header();
 254   AliGenEventHeader* genHeader = header->GenEventHeader();
 255
 256   TArrayF vertex;
 257   genHeader->PrimaryVertex(vertex);
 258   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 259     return;
 260   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 261   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 262   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 263
 264   Bool_t firstTrack = kTRUE;
 265
 266   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
 267   Int_t nTracks = stack->GetNprimary();
 268   if(pars->GetProcessHits())
 269     nTracks = stack->GetNtrack();
 270
 271   for(Int_t i = 0 ;i<nTracks;i++) {
 272     particle = (AliMCParticle*) mcEvent->GetTrack(i);
 273     if(!particle)
 274       continue;
 275
 276     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
 277       hPrimary->Fill(particle->Eta());
 278
 279
 280       TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
 281       hPrimVtxBin->Fill(particle->Eta());
 282       if(firstTrack) {
 283         fNMCevents.Fill(vertexBin);
 284         firstTrack = kFALSE;
 285       }
 286
 287     }
 288     if(pars->GetProcessHits()) {
 289
 290       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
 291
 292         AliTrackReference* ref = particle->GetTrackReference(j);
 293         UShort_t det,sec,strip;
 294         Char_t   ring;
 295         if(ref->DetectorId() != AliTrackReference::kFMD)
 296           continue;
 297         AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
 298         Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
 299         if(particle->Charge() != 0 && i != thisStripTrack ) {
 300           //Double_t x,y,z;
 301
 302           Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
 303           TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
 304
 305
 306           hHits->Fill(eta);
 307
 308           Float_t nstrips = (ring =='O' ? 256 : 512);
 309
 310           fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
 311
 312           if(strip >0)
 313             fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
 314           if(strip < (nstrips - 1))
 315             fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
 316
 317         }
 318
 319
 320       }
 321
 322
 323     }
 324
 325
 326
 327   }
 328 }
 329 //_____________________________________________________________________
 330 //
 331 //
 332 // EOF