The revisited SPD cluster analysis in the new AOD framework and the corrections for...
[u/mrichter/AliRoot.git] / 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 //