Next round of upgrades and cleanups of code. The code is now more streamlined and...
[u/mrichter/AliRoot.git] / 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   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           
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 //