FMD/analysis/AliFMDAnalysisTaskDensity.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 "TAxis.h"
  10 #include "TH2F.h"
  11 #include "AliFMDAnalysisTaskDensity.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 "AliESDVertex.h"
  20 #include "TMath.h"
  21 #include "AliFMDAnaParameters.h"
  22 #include "AliFMDGeometry.h"
  23 #include "AliFMDRing.h"
  24
  25 ClassImp(AliFMDAnalysisTaskDensity)
  26
  27 //_____________________________________________________________________
  28 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
  29 : fDebug(0),
  30   fOutputList(),
  31   fArray(),
  32   fESD(0x0),
  33   fVertexString(),
  34   fVertex(0),
  35   fStandalone(kTRUE)
  36 {
  37   // Default constructor
  38   DefineInput (0, AliESDFMD::Class());
  39   DefineInput (1, AliESDVertex::Class());
  40   DefineOutput(0,TList::Class());
  41 }
  42 //_____________________________________________________________________
  43 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
  44     AliAnalysisTask(name, "Density"),
  45     fDebug(0),
  46     fOutputList(0),
  47     fArray(),
  48     fESD(0x0),
  49     fVertexString(),
  50     fVertex(0),
  51     fStandalone(kTRUE)
  52 {
  53   fStandalone = SE;
  54   if(fStandalone) {
  55     DefineInput (0, AliESDFMD::Class());
  56     DefineInput (1, AliESDVertex::Class());
  57     DefineOutput(0, TList::Class());
  58   }
  59 }
  60 //_____________________________________________________________________
  61 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
  62 {
  63   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  64
  65   fArray.SetName("FMD");
  66   fArray.SetOwner();
  67   if(!fOutputList)
  68     fOutputList = new TList();
  69   fOutputList->SetName("density_list");
  70
  71   TH2F* hMult = 0;
  72
  73   Int_t nVtxbins = pars->GetNvtxBins();
  74
  75   for(Int_t det =1; det<=3;det++)
  76     {
  77       TObjArray* detArray = new TObjArray();
  78       detArray->SetName(Form("FMD%d",det));
  79       fArray.AddAtAndExpand(detArray,det);
  80       Int_t nRings = (det==1 ? 1 : 2);
  81       for(Int_t ring = 0;ring<nRings;ring++)
  82         {
  83           Char_t ringChar = (ring == 0 ? 'I' : 'O');
  84           Int_t  nSec     = (ring == 0 ? 20 : 40);
  85
  86           TObjArray* vtxArray = new TObjArray();
  87           vtxArray->SetName(Form("FMD%d%c",det,ringChar));
  88           detArray->AddAtAndExpand(vtxArray,ring);
  89           for(Int_t i = 0; i< nVtxbins; i++) {
  90             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
  91
  92             hMult  = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
  93                               hBg->GetNbinsX(),
  94                               hBg->GetXaxis()->GetXmin(),
  95                               hBg->GetXaxis()->GetXmax(),
  96                               nSec, 0, 2*TMath::Pi());
  97             vtxArray->AddAtAndExpand(hMult,i);
  98
  99           }
 100         }
 101     }
 102
 103   fOutputList->Add(&fArray);
 104   fOutputList->Add(&fVertexString);
 105
 106
 107 }
 108 //_____________________________________________________________________
 109 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
 110 {
 111   if(fStandalone) {
 112     fESD    = (AliESDFMD*)GetInputData(0);
 113     fVertex = (AliESDVertex*)GetInputData(1);
 114   }
 115 }
 116 //_____________________________________________________________________
 117 void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
 118 {
 119   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 120   AliFMDGeometry* geo       = AliFMDGeometry::Instance();
 121
 122   //AliESDFMD*   fmd = fESD->GetFMDData();
 123
 124   Double_t vertex[3];
 125   fVertex->GetXYZ(vertex);
 126   // Z Vtx cut
 127   if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
 128     return;
 129   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 130   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
 131
 132   Int_t vtxbin = (Int_t)vertexBinDouble;
 133
 134   fVertexString.SetString(Form("%d",vtxbin));
 135   //Reset everything
 136   for(UShort_t det=1;det<=3;det++) {
 137     TObjArray* detArray = (TObjArray*)fArray.At(det);
 138     Int_t nRings = (det==1 ? 1 : 2);
 139     for (UShort_t ir = 0; ir < nRings; ir++) {
 140       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 141
 142       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 143       hMult->Reset();
 144     }
 145
 146   }
 147
 148
 149   for(UShort_t det=1;det<=3;det++) {
 150     TObjArray* detArray = (TObjArray*)fArray.At(det);
 151     Int_t nRings = (det==1 ? 1 : 2);
 152     for (UShort_t ir = 0; ir < nRings; ir++) {
 153       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 154
 155       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 156       Char_t   ring = (ir == 0 ? 'I' : 'O');
 157       UShort_t nsec = (ir == 0 ? 20  : 40);
 158       UShort_t nstr = (ir == 0 ? 512 : 256);
 159       for(UShort_t sec =0; sec < nsec;  sec++)  {
 160         for(UShort_t strip = 0; strip < nstr; strip++) {
 161           Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
 162           if(mult < 1 || mult == AliESDFMD::kInvalidMult) continue;
 163           Float_t eta = fESD->Eta(det,ring,sec,strip);
 164           Double_t x,y,z;
 165           geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 166           Float_t phi = TMath::ATan2(y,x);
 167           if(phi<0)
 168             phi = phi+2*TMath::Pi();
 169           Float_t correction = GetAcceptanceCorrection(ring,strip);
 170           if(correction) mult = mult / correction;
 171           hMult->Fill(eta,phi,mult);
 172
 173         }
 174       }
 175     }
 176
 177
 178   }
 179   if(fStandalone) {
 180     PostData(0, fOutputList);
 181   }
 182
 183 }
 184 //_____________________________________________________________________
 185 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
 186 {
 187   AliFMDRing fmdring(ring);
 188   fmdring.Init();
 189   Float_t   rad       = fmdring.GetMaxR()-fmdring.GetMinR();
 190   Float_t   segment   = rad / fmdring.GetNStrips();
 191   Float_t   radius    = fmdring.GetMinR() + segment*strip;
 192
 193   Float_t   basearea1 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power(radius,2);
 194   Float_t   basearea2 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power((radius-segment),2);
 195   Float_t   basearea  = basearea1 - basearea2;
 196
 197   Float_t   area1     = 0.5*fmdring.GetStripLength(strip)*TMath::Power(radius,2);
 198   Float_t   area2     = 0.5*fmdring.GetStripLength(strip)*TMath::Power((radius-segment),2);
 199   Float_t   area      = area1 - area2;
 200
 201   Float_t correction = area/basearea;
 202
 203   return correction;
 204 }
 205
 206 //
 207 //EOF
 208 //