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