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 "AliFMDAnaCalibEnergyDistribution.h"
  21 #include "AliESDVertex.h"
  22 #include "TMath.h"
  23 #include "AliFMDAnaParameters.h"
  24 //#include "AliFMDParameters.h"
  25 //#include "AliFMDGeometry.h"
  26 //#include "AliFMDRing.h"
  27
  28 ClassImp(AliFMDAnalysisTaskDensity)
  29
  30 //_____________________________________________________________________
  31 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
  32 : fDebug(0),
  33   fOutputList(),
  34   fESD(0x0),
  35   fVertexString(),
  36   fVertex(0),
  37   fStandalone(kTRUE),
  38   fStatus(kTRUE)
  39 {
  40   // Default constructor
  41   DefineInput (0, AliESDFMD::Class());
  42   DefineInput (1, AliESDVertex::Class());
  43   DefineOutput(0,TList::Class());
  44 }
  45 //_____________________________________________________________________
  46 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
  47     AliAnalysisTask(name, "Density"),
  48     fDebug(0),
  49     fOutputList(0),
  50     fESD(0x0),
  51     fVertexString(),
  52     fVertex(0),
  53     fStandalone(kTRUE),
  54     fStatus(kTRUE)
  55 {
  56   fStandalone = SE;
  57   if(fStandalone) {
  58     DefineInput (0, AliESDFMD::Class());
  59     DefineInput (1, AliESDVertex::Class());
  60     DefineOutput(0, TList::Class());
  61   }
  62
  63 }
  64 //_____________________________________________________________________
  65 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
  66 {
  67   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  68
  69
  70
  71   if(!fOutputList)
  72     fOutputList = new TList();
  73   fOutputList->SetName("density_list");
  74
  75   fOutputList->Add(&fVertexString);
  76
  77
  78
  79   TH2F* hMult = 0;
  80
  81   Int_t nVtxbins = pars->GetNvtxBins();
  82
  83
  84   for(Int_t det =1; det<=3;det++)
  85     {
  86       Int_t nRings = (det==1 ? 1 : 2);
  87       for(Int_t ring = 0;ring<nRings;ring++)
  88         {
  89           Char_t ringChar = (ring == 0 ? 'I' : 'O');
  90           Int_t  nSec     = (ring == 0 ? 20 : 40);
  91
  92           for(Int_t i = 0; i< nVtxbins; i++) {
  93             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
  94
  95             hMult  = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
  96                               hBg->GetNbinsX(),
  97                               hBg->GetXaxis()->GetXmin(),
  98                               hBg->GetXaxis()->GetXmax(),
  99                               nSec, 0, 2*TMath::Pi());
 100             hMult->Sumw2();
 101
 102             fOutputList->Add(hMult);
 103           }
 104         }
 105     }
 106
 107
 108
 109
 110 }
 111 //_____________________________________________________________________
 112 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
 113 {
 114   if(fStandalone) {
 115     fESD    = (AliESDFMD*)GetInputData(0);
 116     fVertex = (AliESDVertex*)GetInputData(1);
 117   }
 118 }
 119 //_____________________________________________________________________
 120 void AliFMDAnalysisTaskDensity::Init() {
 121   //TFile f("/home/canute/ALICE/FMDanalysis/FirstAnalysis/energydistributions_0_0_1_0_0_0.root");
 122   //fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(f.Get("energydistributions"));
 123
 124 }
 125 //_____________________________________________________________________
 126 void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
 127 {
 128   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 129   // AliFMDGeometry* geo       = AliFMDGeometry::Instance();
 130
 131   //AliESDFMD*   fmd = fESD->GetFMDData();
 132
 133
 134   Double_t vertex[3];
 135   fVertex->GetXYZ(vertex);
 136   // Z Vtx cut
 137
 138   if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
 139     fStatus = kFALSE;
 140     return;
 141   }
 142   else
 143     fStatus = kTRUE;
 144
 145   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 146   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
 147
 148   Int_t vtxbin = (Int_t)vertexBinDouble;
 149
 150   fVertexString.SetString(Form("%d",vtxbin));
 151   //Reset everything
 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       Char_t   ring = (ir == 0 ? 'I' : 'O');
 156       TH2F* hMult   = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
 157       hMult->Reset();
 158     }
 159
 160   }
 161
 162
 163   for(UShort_t det=1;det<=3;det++) {
 164     Int_t nRings = (det==1 ? 1 : 2);
 165     for (UShort_t ir = 0; ir < nRings; ir++) {
 166
 167       Char_t   ring = (ir == 0 ? 'I' : 'O');
 168       TH2F* hMult   = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
 169
 170       UShort_t nsec = (ir == 0 ? 20  : 40);
 171       UShort_t nstr = (ir == 0 ? 512 : 256);
 172       /*
 173      TH1F* hEnergyDist       = pars->GetEmptyEnergyDistribution(det,ring);
 174       TF1*  fitFunc           = hEnergyDist->GetFunction("FMDfitFunc");
 175       TH1F* hSignalDist    = pars->GetRingEnergyDistribution(det, ring);
 176       TF1*  fitFuncSignal  = hSignalDist->GetFunction("FMDfitFunc");
 177       */
 178       for(UShort_t sec =0; sec < nsec;  sec++)  {
 179         for(UShort_t strip = 0; strip < nstr; strip++) {
 180           Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
 181
 182           if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
 183
 184           Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 185           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 186
 187           Float_t mult_cut = 0.3;//pars->GetMPV(det,ring,eta);// - pars->GetSigma(det,ring,eta);//0.25;//0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2;
 188           // if(ring == 'I')
 189           //  mult_cut = 0.10;
 190           //Float_t mult_cut = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
 191           Float_t nParticles = 0;
 192           if(fESD->GetUniqueID() == kTRUE) {
 193             //proton + proton
 194
 195             if(mult > mult_cut) {
 196               nParticles = 1.;
 197             }
 198           }
 199           else {
 200
 201             //Pb+Pb
 202             Float_t mpv   = pars->GetMPV(det,ring,eta);
 203             Float_t sigma = pars->GetSigma(det,ring,eta);
 204             Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
 205             Float_t beta  = pars->Get3MIPWeight(det,ring,eta);
 206
 207             Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
 208               alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 209               beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 210             Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
 211               2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 212               3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 213
 214
 215             if(mult > mult_cut) {
 216               if(sumCor) nParticles = weight / sumCor;
 217               else nParticles = 1;
 218
 219             }
 220             //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
 221
 222           }
 223
 224
 225
 226
 227
 228           Float_t correction = GetAcceptanceCorrection(ring,strip);
 229
 230           //std::cout<<"before "<<correction<<std::endl;
 231           if(fESD->GetUniqueID() == kTRUE) {
 232             TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring);
 233
 234             if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0)
 235               correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta));
 236
 237           }
 238
 239           //Dead channel correction:
 240           TH1F* hFMDDeadCorrection = pars->GetFMDDeadCorrection(vtxbin);
 241           if(hFMDDeadCorrection)
 242             if(hFMDDeadCorrection->GetBinContent(hFMDDeadCorrection->FindBin(eta)) != 0)
 243               correction = correction*hFMDDeadCorrection->GetBinContent(hFMDDeadCorrection->FindBin(eta));
 244
 245
 246           //std::cout<<det<<"   "<<ring<<"    "<<sec<<"    "<<strip<<"    "<<vertex[2]<<"   "<<eta<<std::endl;
 247
 248           //Float_t signal = hSignalDist->Integral(hSignalDist->FindBin(0.5),hSignalDist->FindBin(2)) ;//pars->GetConstant(det,ring,eta);
 249           //Float_t bkg =  hEnergyDist->GetBinContent(hEnergyDist->FindBin(mult));
 250           //Float_t signal =  hSignalDist->GetBinContent(hSignalDist->FindBin(mult));
 251           /*
 252           if(fitFunc && fitFuncSignal && pars->IsRealData()) {
 253             Float_t bkg = fitFunc->Eval(mult);
 254             Float_t signal = fitFuncSignal->Eval(mult);
 255
 256             if(signal > 0) {
 257               correction = correction/(1-(bkg/signal));
 258               //test
 259               // if(TMath::Abs(eta)<3) correction = 1.15*correction;
 260               // if(det == 2 && ring == 'I') correction = correction/(1-bkg/signal);
 261               //if(det == 2 && ring == 'O') correction = correction/(1-bkg/signal);
 262               //if(det == 3 && ring == 'I') correction = correction/(1-bkg/signal);
 263               //if(det == 3 && ring == 'O') correction = correction/(1-bkg/signal);
 264
 265             }
 266             }*/
 267           if(correction) nParticles = nParticles / correction;
 268           if(nParticles > 0)
 269             hMult->Fill(eta,phi,nParticles);
 270
 271
 272         }
 273       }
 274
 275     }
 276
 277
 278
 279   }
 280
 281
 282
 283   if(fStandalone) {
 284     PostData(0, fOutputList);
 285   }
 286
 287 }
 288 //_____________________________________________________________________
 289 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
 290 {
 291   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 292
 293   //AliFMDRing fmdring(ring);
 294   //fmdring.Init();
 295   Float_t   rad       = pars->GetMaxR(ring)-pars->GetMinR(ring);
 296   Float_t   nstrips   = (ring == 'I' ? 512 : 256);
 297   Float_t   segment   = rad / nstrips;
 298   Float_t   radius    = pars->GetMinR(ring) + segment*strip;
 299
 300   Float_t   basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
 301   Float_t   basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
 302   Float_t   basearea  = basearea1 - basearea2;
 303
 304   Float_t   area1     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
 305   Float_t   area2     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
 306   Float_t   area      = area1 - area2;
 307
 308   Float_t correction = area/basearea;
 309
 310   return correction;
 311 }
 312 //_____________________________________________________________________
 313 //
 314 //EOF
 315 //