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   fStatus(kTRUE)
  37 {
  38   // Default constructor
  39   DefineInput (0, AliESDFMD::Class());
  40   DefineInput (1, AliESDVertex::Class());
  41   DefineOutput(0,TList::Class());
  42 }
  43 //_____________________________________________________________________
  44 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
  45     AliAnalysisTask(name, "Density"),
  46     fDebug(0),
  47     fOutputList(0),
  48     fArray(),
  49     fESD(0x0),
  50     fVertexString(),
  51     fVertex(0),
  52     fStandalone(kTRUE),
  53     fStatus(kTRUE)
  54 {
  55   fStandalone = SE;
  56   if(fStandalone) {
  57     DefineInput (0, AliESDFMD::Class());
  58     DefineInput (1, AliESDVertex::Class());
  59     DefineOutput(0, TList::Class());
  60   }
  61 }
  62 //_____________________________________________________________________
  63 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
  64 {
  65   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  66
  67   fArray.SetName("FMD");
  68   fArray.SetOwner();
  69   if(!fOutputList)
  70     fOutputList = new TList();
  71   fOutputList->SetName("density_list");
  72
  73   fOutputList->Add(&fArray);
  74   fOutputList->Add(&fVertexString);
  75
  76   TH2F* hMult = 0;
  77
  78   Int_t nVtxbins = pars->GetNvtxBins();
  79
  80   for(Int_t det =1; det<=3;det++)
  81     {
  82       TObjArray* detArray = new TObjArray();
  83       detArray->SetName(Form("FMD%d",det));
  84       fArray.AddAtAndExpand(detArray,det);
  85       Int_t nRings = (det==1 ? 1 : 2);
  86       for(Int_t ring = 0;ring<nRings;ring++)
  87         {
  88           Char_t ringChar = (ring == 0 ? 'I' : 'O');
  89           Int_t  nSec     = (ring == 0 ? 20 : 40);
  90
  91           TObjArray* vtxArray = new TObjArray();
  92           vtxArray->SetName(Form("FMD%d%c",det,ringChar));
  93           detArray->AddAtAndExpand(vtxArray,ring);
  94           for(Int_t i = 0; i< nVtxbins; i++) {
  95             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
  96
  97             hMult  = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
  98                               hBg->GetNbinsX(),
  99                               hBg->GetXaxis()->GetXmin(),
 100                               hBg->GetXaxis()->GetXmax(),
 101                               nSec, 0, 2*TMath::Pi());
 102
 103             vtxArray->AddAtAndExpand(hMult,i);
 104           }
 105         }
 106     }
 107
 108
 109
 110
 111 }
 112 //_____________________________________________________________________
 113 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
 114 {
 115   if(fStandalone) {
 116     fESD    = (AliESDFMD*)GetInputData(0);
 117     fVertex = (AliESDVertex*)GetInputData(1);
 118   }
 119 }
 120 //_____________________________________________________________________
 121 void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
 122 {
 123   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 124   AliFMDGeometry* geo       = AliFMDGeometry::Instance();
 125
 126   //AliESDFMD*   fmd = fESD->GetFMDData();
 127
 128   Double_t vertex[3];
 129   fVertex->GetXYZ(vertex);
 130   // Z Vtx cut
 131   if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
 132     fStatus = kFALSE;
 133     return;
 134   }
 135   else
 136     fStatus = kTRUE;
 137
 138   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 139   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
 140
 141   Int_t vtxbin = (Int_t)vertexBinDouble;
 142
 143   fVertexString.SetString(Form("%d",vtxbin));
 144   //Reset everything
 145   for(UShort_t det=1;det<=3;det++) {
 146     TObjArray* detArray = (TObjArray*)fArray.At(det);
 147     Int_t nRings = (det==1 ? 1 : 2);
 148     for (UShort_t ir = 0; ir < nRings; ir++) {
 149       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 150
 151       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 152       hMult->Reset();
 153     }
 154
 155   }
 156
 157
 158   for(UShort_t det=1;det<=3;det++) {
 159     TObjArray* detArray = (TObjArray*)fArray.At(det);
 160     Int_t nRings = (det==1 ? 1 : 2);
 161     for (UShort_t ir = 0; ir < nRings; ir++) {
 162       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 163
 164       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 165
 166       Char_t   ring = (ir == 0 ? 'I' : 'O');
 167       UShort_t nsec = (ir == 0 ? 20  : 40);
 168       UShort_t nstr = (ir == 0 ? 512 : 256);
 169
 170       for(UShort_t sec =0; sec < nsec;  sec++)  {
 171         for(UShort_t strip = 0; strip < nstr; strip++) {
 172           Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
 173           Float_t eta = fESD->Eta(det,ring,sec,strip);
 174           Float_t mult_cut = 0.2;
 175           if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
 176           //Particle number cut goes here...
 177           Float_t nParticles = 0;
 178           if(fESD->GetUniqueID() == kTRUE) {
 179             //proton + proton
 180             if(mult > mult_cut)
 181               nParticles = 1;
 182           }
 183           else {
 184
 185             //Pb+Pb
 186             Float_t mpv   = pars->GetMPV(det,ring,eta);
 187             Float_t sigma = pars->GetSigma(det,ring,eta);
 188             Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
 189             Float_t beta  = pars->Get3MIPWeight(det,ring,eta);
 190
 191             Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
 192               alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 193               beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 194             Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
 195               2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 196               3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 197
 198
 199             if(mult > 0){//mult_cut) {
 200               if(sumCor) nParticles = weight / sumCor;
 201               else nParticles = 1;
 202             }
 203             //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
 204
 205           }
 206
 207
 208           Double_t x,y,z;
 209           geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 210           Float_t phi = TMath::ATan2(y,x);
 211           if(phi<0)
 212             phi = phi+2*TMath::Pi();
 213           Float_t correction = GetAcceptanceCorrection(ring,strip);
 214           if(correction) nParticles = nParticles / correction;
 215           hMult->Fill(eta,phi,nParticles);
 216
 217
 218         }
 219       }
 220
 221     }
 222
 223
 224
 225   }
 226   if(fStandalone) {
 227     PostData(0, fOutputList);
 228   }
 229
 230 }
 231 //_____________________________________________________________________
 232 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
 233 {
 234   AliFMDRing fmdring(ring);
 235   fmdring.Init();
 236   Float_t   rad       = fmdring.GetMaxR()-fmdring.GetMinR();
 237   Float_t   segment   = rad / fmdring.GetNStrips();
 238   Float_t   radius    = fmdring.GetMinR() + segment*strip;
 239
 240   Float_t   basearea1 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power(radius,2);
 241   Float_t   basearea2 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power((radius-segment),2);
 242   Float_t   basearea  = basearea1 - basearea2;
 243
 244   Float_t   area1     = 0.5*fmdring.GetStripLength(strip)*TMath::Power(radius,2);
 245   Float_t   area2     = 0.5*fmdring.GetStripLength(strip)*TMath::Power((radius-segment),2);
 246   Float_t   area      = area1 - area2;
 247
 248   Float_t correction = area/basearea;
 249
 250   return correction;
 251 }
 252
 253 //
 254 //EOF
 255 //