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