PWG2/FORWARD/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 "TF1.h"
  12 #include "AliFMDAnalysisTaskDensity.h"
  13 #include "AliAnalysisManager.h"
  14 #include "AliESDFMD.h"
  15 #include "AliESDEvent.h"
  16 #include "AliAODEvent.h"
  17 #include "AliAODHandler.h"
  18 #include "AliMCEventHandler.h"
  19 #include "AliStack.h"
  20 #include "AliESDVertex.h"
  21 #include "TMath.h"
  22 #include "AliFMDAnaParameters.h"
  23 //#include "AliFMDParameters.h"
  24 //#include "AliFMDGeometry.h"
  25 //#include "AliFMDRing.h"
  26
  27 ClassImp(AliFMDAnalysisTaskDensity)
  28
  29 //_____________________________________________________________________
  30 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
  31 : fDebug(0),
  32   fOutputList(),
  33   fESD(0x0),
  34   fVertexString(),
  35   fVertex(0),
  36   fStandalone(kTRUE),
  37   fStatus(kTRUE)
  38 {
  39   // Default constructor
  40   DefineInput (0, AliESDFMD::Class());
  41   DefineInput (1, AliESDVertex::Class());
  42   DefineOutput(0,TList::Class());
  43 }
  44 //_____________________________________________________________________
  45 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
  46     AliAnalysisTask(name, "Density"),
  47     fDebug(0),
  48     fOutputList(0),
  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 //_____________________________________________________________________
  64 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
  65 {
  66   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  67
  68   if(!fOutputList)
  69     fOutputList = new TList();
  70   fOutputList->SetName("density_list");
  71
  72   fOutputList->Add(&fVertexString);
  73
  74   TH2F* hMult = 0;
  75
  76   Int_t nVtxbins = pars->GetNvtxBins();
  77
  78   for(Int_t det =1; det<=3;det++)
  79     {
  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           for(Int_t i = 0; i< nVtxbins; i++) {
  87             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
  88
  89             hMult  = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
  90                               hBg->GetNbinsX(),
  91                               hBg->GetXaxis()->GetXmin(),
  92                               hBg->GetXaxis()->GetXmax(),
  93                               nSec, 0, 2*TMath::Pi());
  94
  95             fOutputList->Add(hMult);
  96           }
  97         }
  98     }
  99
 100
 101
 102
 103 }
 104 //_____________________________________________________________________
 105 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
 106 {
 107   if(fStandalone) {
 108     fESD    = (AliESDFMD*)GetInputData(0);
 109     fVertex = (AliESDVertex*)GetInputData(1);
 110   }
 111 }
 112
 113
 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     fStatus = kFALSE;
 129     return;
 130   }
 131   else
 132     fStatus = kTRUE;
 133
 134   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 135   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
 136
 137   Int_t vtxbin = (Int_t)vertexBinDouble;
 138
 139   fVertexString.SetString(Form("%d",vtxbin));
 140   //Reset everything
 141   for(UShort_t det=1;det<=3;det++) {
 142     Int_t nRings = (det==1 ? 1 : 2);
 143     for (UShort_t ir = 0; ir < nRings; ir++) {
 144       Char_t   ring = (ir == 0 ? 'I' : 'O');
 145       TH2F* hMult   = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
 146       hMult->Reset();
 147     }
 148
 149   }
 150
 151
 152   for(UShort_t det=1;det<=3;det++) {
 153     Int_t nRings = (det==1 ? 1 : 2);
 154     for (UShort_t ir = 0; ir < nRings; ir++) {
 155
 156       Char_t   ring = (ir == 0 ? 'I' : 'O');
 157       TH2F* hMult   = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
 158
 159       UShort_t nsec = (ir == 0 ? 20  : 40);
 160       UShort_t nstr = (ir == 0 ? 512 : 256);
 161
 162       for(UShort_t sec =0; sec < nsec;  sec++)  {
 163         for(UShort_t strip = 0; strip < nstr; strip++) {
 164           Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
 165
 166           if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
 167
 168           Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 169           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 170
 171           Float_t mult_cut = 0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2;
 172             if(ring == 'I')
 173               mult_cut = 0.10;
 174           //Float_t mult_cut = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
 175           Float_t nParticles = 0;
 176           if(fESD->GetUniqueID() == kTRUE) {
 177             //proton + proton
 178
 179             if(mult > mult_cut) {
 180               nParticles = 1;
 181
 182             }
 183           }
 184           else {
 185
 186             //Pb+Pb
 187             Float_t mpv   = pars->GetMPV(det,ring,eta);
 188             Float_t sigma = pars->GetSigma(det,ring,eta);
 189             Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
 190             Float_t beta  = pars->Get3MIPWeight(det,ring,eta);
 191
 192             Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
 193               alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 194               beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 195             Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
 196               2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 197               3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 198
 199
 200             if(mult > mult_cut) {
 201               if(sumCor) nParticles = weight / sumCor;
 202               else nParticles = 1;
 203
 204             }
 205             //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
 206
 207           }
 208
 209
 210
 211
 212
 213           Float_t correction = GetAcceptanceCorrection(ring,strip);
 214
 215           //std::cout<<"before "<<correction<<std::endl;
 216           if(fESD->GetUniqueID() == kTRUE) {
 217             TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring);
 218
 219             if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0)
 220               correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta));
 221
 222           }
 223
 224           if(correction) nParticles = nParticles / correction;
 225           if(nParticles > 0)
 226             hMult->Fill(eta,phi,nParticles);
 227
 228
 229         }
 230       }
 231
 232     }
 233
 234
 235
 236   }
 237
 238
 239   if(fStandalone) {
 240     PostData(0, fOutputList);
 241   }
 242
 243 }
 244 //_____________________________________________________________________
 245 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
 246 {
 247   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 248
 249   //AliFMDRing fmdring(ring);
 250   //fmdring.Init();
 251   Float_t   rad       = pars->GetMaxR(ring)-pars->GetMinR(ring);
 252   Float_t   nstrips   = (ring == 'I' ? 512 : 256);
 253   Float_t   segment   = rad / nstrips;
 254   Float_t   radius    = pars->GetMinR(ring) + segment*strip;
 255
 256   Float_t   basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
 257   Float_t   basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
 258   Float_t   basearea  = basearea1 - basearea2;
 259
 260   Float_t   area1     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
 261   Float_t   area2     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
 262   Float_t   area      = area1 - area2;
 263
 264   Float_t correction = area/basearea;
 265
 266   return correction;
 267 }
 268 //_____________________________________________________________________
 269 //
 270 //EOF
 271 //