PWG2/FORWARD/analysis2/AliCentralMultiplicityTask.cxx

   1 //====================================================================
   2 //
   3 // Base class for classes that calculate the multiplicity in the
   4 // central region event-by-event
   5 //
   6 // Inputs:
   7 //   - AliESDEvent
   8 //
   9 // Outputs:
  10 //   - AliAODCentralMult
  11 //
  12 // Histograms
  13 //
  14 // Corrections used
  15 #include "AliCentralMultiplicityTask.h"
  16 #include "AliForwardCorrectionManager.h"
  17 #include "AliForwardUtil.h"
  18 #include "AliLog.h"
  19 #include "AliAODHandler.h"
  20 #include "AliInputEventHandler.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliESDEvent.h"
  24 #include "AliMultiplicity.h"
  25 #include <TROOT.h>
  26 #include <TFile.h>
  27 #include <iostream>
  28 #include <iomanip>
  29
  30 //====================================================================
  31 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
  32   : AliAnalysisTaskSE(name),
  33     fData(0),
  34     fList(0),
  35     fAODCentral(kFALSE),
  36     fManager()
  37 {
  38
  39   DefineOutput(1, TList::Class());
  40 }
  41 //____________________________________________________________________
  42 void AliCentralMultiplicityTask::UserCreateOutputObjects() {
  43
  44   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  45   AliAODHandler*      ah =
  46     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  47   if (!ah) AliFatal("No AOD output handler set in analysis manager");
  48
  49
  50   TObject* obj = &fAODCentral;
  51   ah->AddBranch("AliAODCentralMult", &obj);
  52
  53   fList = new TList();
  54   PostData(1,fList);
  55
  56 }
  57 //____________________________________________________________________
  58 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) {
  59
  60   AliESDInputHandler* eventHandler =
  61     dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
  62                                        ->GetInputEventHandler());
  63   if (!eventHandler) {
  64     AliWarning("No inputhandler found for this event!");
  65     return;}
  66
  67   AliESDEvent* esd = eventHandler->GetEvent();
  68
  69   //Selecting only events with |valid vertex| < 10 cm
  70   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
  71   if (!vertex) return;
  72   if(!(vertex->GetStatus())) return;
  73   if(vertex->GetNContributors() <= 0) return ;
  74   if(vertex->GetZRes() > 0.1 ) return;
  75   Double_t vertexXYZ[3]={0,0,0};
  76   vertex->GetXYZ(vertexXYZ);
  77   if(TMath::Abs(vertexXYZ[2]) > 10) return;
  78
  79   Double_t delta = 2 ;
  80   Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
  81   Int_t vtxbin = (Int_t)vertexBinDouble +1 ; //HHD: The vtxbins are 1-10, not 0-9
  82
  83   // Make sure AOD is filled
  84   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  85   AliAODHandler*      ah =
  86     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  87   if (!ah)
  88     AliFatal("No AOD output handler set in analysis manager");
  89
  90   ah->SetFillAOD(kTRUE);
  91
  92   //Doing analysis
  93   fAODCentral.Clear("");
  94   TH2D *aodHist = &(fAODCentral.GetHistogram());
  95
  96   const AliMultiplicity* spdmult = esd->GetMultiplicity();
  97   //Filling clusters in layer 1 used for tracklets...
  98   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
  99     aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 100   }
 101   //...and then the unused ones in layer 1
 102   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 103     aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
 104                    spdmult->GetPhiSingle(j));
 105   }
 106
 107   // Corrections
 108
 109   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
 110   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
 111
 112   aodHist->Divide(hSecMap);
 113
 114   for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
 115     Float_t acccor = hAcceptance->GetBinContent(nx);
 116
 117     Bool_t etabinSeen = kFALSE;
 118     for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
 119       Float_t aodvalue = aodHist->GetBinContent(nx,ny);
 120       Float_t seccor = hSecMap->GetBinContent(nx,ny);
 121       if(seccor > 0.5) etabinSeen = kTRUE;
 122       if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
 123
 124       Float_t aodnew   = aodvalue / acccor ;
 125       aodHist->SetBinContent(nx,ny, aodnew);
 126       Float_t error = aodnew*TMath::Sqrt(TMath::Power(aodHist->GetBinError(nx,ny)/aodvalue,2) +
 127                                          TMath::Power(hAcceptance->GetBinError(nx)/acccor,2) );
 128       aodHist->SetBinError(nx,ny,error);
 129
 130     }
 131     //Filling underflow bin if we eta bin is in range
 132     if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
 133
 134   }
 135
 136   PostData(1,fList);
 137
 138 }
 139 //____________________________________________________________________
 140 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/) {
 141
 142
 143
 144 }
 145 //____________________________________________________________________
 146 void
 147 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
 148 {
 149
 150 }
 151 //____________________________________________________________________
 152 AliCentralMultiplicityTask::Manager::Manager() :
 153   fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
 154   fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
 155   fAcceptance(),
 156   fSecmap(),
 157   fAcceptanceName("centralacceptance"),
 158   fSecMapName("centralsecmap")
 159 {
 160
 161
 162 }
 163 //____________________________________________________________________
 164 const char* AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
 165                                                              UShort_t  sys,
 166                                                              UShort_t  sNN,
 167                                                              Short_t   field) {
 168
 169   TString fname = "";
 170
 171   switch(what) {
 172   case 0:
 173     fname.Append(fSecMapName.Data());
 174     break;
 175   case 1:
 176     fname.Append(fAcceptanceName.Data());
 177     break;
 178   default:
 179     printf("Invalid indentifier for central object, must be 0 or 1!");
 180     break;
 181   }
 182   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
 183                     AliForwardUtil::CollisionSystemString(sys),
 184                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
 185
 186   return fname.Data();
 187
 188 }
 189 //____________________________________________________________________
 190 void AliCentralMultiplicityTask::Manager::Init( UShort_t  sys,
 191                                                 UShort_t  sNN,
 192                                                 Short_t   field) {
 193
 194   TFile fsec(GetFullFileName(0,sys,sNN,field));
 195   fSecmap = dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
 196   if(!fSecmap) {
 197     printf("no central Secondary Map found!") ;
 198     return;
 199   }
 200   TFile facc(GetFullFileName(1,sys,sNN,field));
 201   fAcceptance = dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
 202 if(!fAcceptance) {
 203   printf("no central Acceptance found!") ;
 204   return;
 205  }
 206
 207 }
 208 //
 209 // EOF
 210 //