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