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 "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   fArray(),
  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     fArray(),
  51     fESD(0x0),
  52     fVertexString(),
  53     fVertex(0),
  54     fStandalone(kTRUE),
  55     fStatus(kTRUE)
  56 {
  57   fStandalone = SE;
  58   if(fStandalone) {
  59     DefineInput (0, AliESDFMD::Class());
  60     DefineInput (1, AliESDVertex::Class());
  61     DefineOutput(0, TList::Class());
  62   }
  63
  64   fFuncPos = new TF1("funcPos","pol1",0,6);
  65   fFuncPos->SetParameters(0.99925,0.00298301);
  66   fFuncNeg = new TF1("funcNeg","pol1",-6,0);
  67   fFuncNeg->SetParameters(0.987583,-0.0101022);
  68
  69
  70 }
  71 //_____________________________________________________________________
  72 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
  73 {
  74   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  75
  76   fArray.SetName("FMD");
  77   fArray.SetOwner();
  78   if(!fOutputList)
  79     fOutputList = new TList();
  80   fOutputList->SetName("density_list");
  81
  82   fOutputList->Add(&fArray);
  83   fOutputList->Add(&fVertexString);
  84
  85   TH2F* hMult = 0;
  86
  87   Int_t nVtxbins = pars->GetNvtxBins();
  88
  89   for(Int_t det =1; det<=3;det++)
  90     {
  91       TObjArray* detArray = new TObjArray();
  92       detArray->SetName(Form("FMD%d",det));
  93       fArray.AddAtAndExpand(detArray,det);
  94       Int_t nRings = (det==1 ? 1 : 2);
  95       for(Int_t ring = 0;ring<nRings;ring++)
  96         {
  97           Char_t ringChar = (ring == 0 ? 'I' : 'O');
  98           Int_t  nSec     = (ring == 0 ? 20 : 40);
  99
 100           TObjArray* vtxArray = new TObjArray();
 101           vtxArray->SetName(Form("FMD%d%c",det,ringChar));
 102           detArray->AddAtAndExpand(vtxArray,ring);
 103           for(Int_t i = 0; i< nVtxbins; i++) {
 104             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
 105
 106             hMult  = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
 107                               hBg->GetNbinsX(),
 108                               hBg->GetXaxis()->GetXmin(),
 109                               hBg->GetXaxis()->GetXmax(),
 110                               nSec, 0, 2*TMath::Pi());
 111
 112             vtxArray->AddAtAndExpand(hMult,i);
 113           }
 114         }
 115     }
 116
 117
 118
 119
 120 }
 121 //_____________________________________________________________________
 122 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
 123 {
 124   if(fStandalone) {
 125     fESD    = (AliESDFMD*)GetInputData(0);
 126     fVertex = (AliESDVertex*)GetInputData(1);
 127   }
 128 }
 129
 130
 131
 132
 133 //_____________________________________________________________________
 134 void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
 135 {
 136   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 137   AliFMDGeometry* geo       = AliFMDGeometry::Instance();
 138
 139   //AliESDFMD*   fmd = fESD->GetFMDData();
 140
 141   Double_t vertex[3];
 142   fVertex->GetXYZ(vertex);
 143   // Z Vtx cut
 144   if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
 145     fStatus = kFALSE;
 146     return;
 147   }
 148   else
 149     fStatus = kTRUE;
 150
 151   Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 152   Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
 153
 154   Int_t vtxbin = (Int_t)vertexBinDouble;
 155
 156
 157
 158   fVertexString.SetString(Form("%d",vtxbin));
 159   //Reset everything
 160   for(UShort_t det=1;det<=3;det++) {
 161     TObjArray* detArray = (TObjArray*)fArray.At(det);
 162     Int_t nRings = (det==1 ? 1 : 2);
 163     for (UShort_t ir = 0; ir < nRings; ir++) {
 164       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 165
 166       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 167       hMult->Reset();
 168     }
 169
 170   }
 171
 172
 173   for(UShort_t det=1;det<=3;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       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
 178
 179       TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
 180
 181       Char_t   ring = (ir == 0 ? 'I' : 'O');
 182       UShort_t nsec = (ir == 0 ? 20  : 40);
 183       UShort_t nstr = (ir == 0 ? 512 : 256);
 184
 185       for(UShort_t sec =0; sec < nsec;  sec++)  {
 186         for(UShort_t strip = 0; strip < nstr; strip++) {
 187           Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
 188           //Float_t eta = fESD->Eta(det,ring,sec,strip);
 189
 190           if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
 191           //Particle number cut goes here...
 192           //Double_t x,y,z;
 193           //geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 194           // Float_t phi = TMath::ATan2(y,x);
 195           // if(phi<0)
 196           //  phi = phi+2*TMath::Pi();
 197
 198           Float_t phi = pars->GetPhiFromSector(det,ring,sec);
 199           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
 200           //std::cout<<phi<<"     "<<phicalc<<std::endl;
 201           //  Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 202          // Float_t   theta = TMath::ATan2(r,z-vertex[2]);
 203           // Float_t   etacalc   = -1*TMath::Log(TMath::Tan(0.5*theta));
 204
 205            //  std::cout<<eta<<"    "<<etacalc<<std::endl;
 206            //eta = etacalc;
 207
 208           Float_t m   = pars->GetMPV(det,ring,eta);
 209           Float_t s   = pars->GetSigma(det,ring,eta);
 210           //AliFMDParameters* recopars = AliFMDParameters::Instance();
 211
 212           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;
 213           if(ring == 'I')
 214             mult_cut = 0.10;
 215           //mult_cut = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
 216
 217           Float_t nParticles = 0;
 218           if(fESD->GetUniqueID() == kTRUE) {
 219             //proton + proton
 220
 221             if(mult > mult_cut) {
 222               nParticles = 1;
 223
 224             }
 225           }
 226           else {
 227
 228             //Pb+Pb
 229             Float_t mpv   = pars->GetMPV(det,ring,eta);
 230             Float_t sigma = pars->GetSigma(det,ring,eta);
 231             Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
 232             Float_t beta  = pars->Get3MIPWeight(det,ring,eta);
 233
 234             Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
 235               alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 236               beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 237             Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
 238               2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
 239               3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
 240
 241
 242             if(mult > mult_cut) {
 243               if(sumCor) nParticles = weight / sumCor;
 244               else nParticles = 1;
 245
 246             }
 247             //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
 248
 249           }
 250
 251
 252
 253
 254
 255           Float_t correction = GetAcceptanceCorrection(ring,strip);
 256
 257           //std::cout<<"before "<<correction<<std::endl;
 258           if(det == 3)
 259             correction = correction / fFuncNeg->Eval(eta);
 260           else
 261             correction = correction / fFuncPos->Eval(eta);
 262
 263           // std::cout<<correction<<std::endl;
 264           if(correction) nParticles = nParticles / correction;
 265           if(nParticles > 0)
 266             hMult->Fill(eta,phi,nParticles);
 267
 268           //if(det == 1 && ring =='I' && nParticles >0)
 269           //if(nParticles > 0)
 270           //  std::cout<<det<<"    "<<ring<<"    "<<sec<<"    "<<strip<<"   "<<mult<<std::endl;
 271
 272         }
 273       }
 274
 275     }
 276
 277
 278
 279   }
 280   if(fStandalone) {
 281     PostData(0, fOutputList);
 282   }
 283
 284 }
 285 //_____________________________________________________________________
 286 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
 287 {
 288   AliFMDRing fmdring(ring);
 289   fmdring.Init();
 290   Float_t   rad       = fmdring.GetMaxR()-fmdring.GetMinR();
 291   Float_t   segment   = rad / fmdring.GetNStrips();
 292   Float_t   radius    = fmdring.GetMinR() + segment*strip;
 293
 294   Float_t   basearea1 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power(radius,2);
 295   Float_t   basearea2 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power((radius-segment),2);
 296   Float_t   basearea  = basearea1 - basearea2;
 297
 298   Float_t   area1     = 0.5*fmdring.GetStripLength(strip)*TMath::Power(radius,2);
 299   Float_t   area2     = 0.5*fmdring.GetStripLength(strip)*TMath::Power((radius-segment),2);
 300   Float_t   area      = area1 - area2;
 301
 302   Float_t correction = area/basearea;
 303
 304   return correction;
 305 }
 306 //_____________________________________________________________________
 307 Float_t AliFMDAnalysisTaskDensity::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
 308 {
 309   Int_t nsec = (ring == 'I' ? 20 : 40);
 310   Float_t basephi = 0;
 311   if(det == 1)
 312     basephi = 1.72787594;
 313   if(det == 2 && ring == 'I')
 314     basephi = 0.15707963;
 315   if(det == 2 && ring == 'O')
 316     basephi = 0.078539818;
 317   if(det == 3 && ring == 'I')
 318     basephi = 2.984513044;
 319   if(det == 3 && ring == 'O')
 320     basephi = 3.06305289;
 321
 322   Float_t step = 2*TMath::Pi() / nsec;
 323   Float_t phi = 0;
 324   if(det == 3)
 325     phi = basephi - sec*step;
 326   else
 327     phi = basephi + sec*step;
 328
 329   if(phi < 0)
 330     phi = phi +2*TMath::Pi();
 331   if(phi > 2*TMath::Pi() )
 332     phi = phi - 2*TMath::Pi();
 333
 334   return phi;
 335 }
 336
 337
 338 //
 339 //EOF
 340 //