Upgraded analysis tasks to conply with the requirements of the Analysis trains
[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   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(fESD->GetUniqueID() == kTRUE) {
259             if(det == 3) 
260               correction = correction / fFuncNeg->Eval(eta);
261             else
262               correction = correction / fFuncPos->Eval(eta);
263           }
264           
265           // std::cout<<correction<<std::endl;
266           if(correction) nParticles = nParticles / correction;
267           if(nParticles > 0)
268             hMult->Fill(eta,phi,nParticles);
269           
270           //if(det == 1 && ring =='I' && nParticles >0)
271           //if(nParticles > 0)
272           //  std::cout<<det<<"    "<<ring<<"    "<<sec<<"    "<<strip<<"   "<<mult<<std::endl;
273           
274         }
275       }
276       
277     }
278     
279         
280   
281   }
282   
283
284   if(fStandalone) {
285     PostData(0, fOutputList); 
286   }
287   
288 }
289 //_____________________________________________________________________
290 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
291 {
292   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
293   
294   //AliFMDRing fmdring(ring);
295   //fmdring.Init();
296   Float_t   rad       = pars->GetMaxR(ring)-pars->GetMinR(ring);
297   Float_t   nstrips   = (ring == 'I' ? 512 : 256);
298   Float_t   segment   = rad / nstrips;
299   Float_t   radius    = pars->GetMinR(ring) + segment*strip;
300   
301   Float_t   basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
302   Float_t   basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
303   Float_t   basearea  = basearea1 - basearea2;
304   
305   Float_t   area1     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
306   Float_t   area2     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
307   Float_t   area      = area1 - area2;
308   
309   Float_t correction = area/basearea;
310   
311   return correction;
312 }
313 //_____________________________________________________________________
314 /*Float_t AliFMDAnalysisTaskDensity::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
315 {
316   Int_t nsec = (ring == 'I' ? 20 : 40);
317   Float_t basephi = 0;
318   if(det == 1) 
319     basephi = 1.72787594; 
320   if(det == 2 && ring == 'I')
321     basephi = 0.15707963;
322   if(det == 2 && ring == 'O')
323     basephi = 0.078539818;
324   if(det == 3 && ring == 'I')
325     basephi = 2.984513044;
326   if(det == 3 && ring == 'O')
327     basephi = 3.06305289;
328   
329   Float_t step = 2*TMath::Pi() / nsec;
330   Float_t phi = 0;
331   if(det == 3)
332     phi = basephi - sec*step;
333   else
334     phi = basephi + sec*step;
335   
336   if(phi < 0) 
337     phi = phi +2*TMath::Pi();
338   if(phi > 2*TMath::Pi() )
339     phi = phi - 2*TMath::Pi();
340   
341   return phi;
342 }
343
344 */
345 //
346 //EOF
347 //