1317a71969d3179d06a31e89fd9ce1c0512838db
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskCollector.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 <TMath.h>
9 #include <iostream>
10
11 #include "AliFMDAnalysisTaskCollector.h"
12 #include "AliAnalysisManager.h"
13 #include "AliESDFMD.h"
14 #include "AliESDEvent.h"
15 #include "AliAODEvent.h"
16 #include "AliAODHandler.h"
17 #include "AliMCEventHandler.h"
18 #include "AliStack.h"
19 #include "AliESDVertex.h"
20 #include "AliFMDAnaParameters.h"
21 //#include "AliFMDGeometry.h"
22
23
24 ClassImp(AliFMDAnalysisTaskCollector)
25
26 //____________________________________________________________________
27 Double_t  AliFMDAnalysisTaskCollector::TripleLandau(Double_t *x, Double_t *par) {
28   
29   Double_t energy        = x[0];
30   Double_t constant = par[0];
31   Double_t mpv      = par[1];
32   Double_t sigma    = par[2];
33   Double_t alpha    = par[3];
34   Double_t beta     = par[4];
35  
36   Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
37                          alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
38                          beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
39   
40   return f;
41 }
42 //____________________________________________________________________
43
44 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
45 : fDebug(0),
46   fOutputList(0),
47   fArray(0),
48   fZvtxDist(0)
49 {
50   // Default constructor
51   
52  
53 }
54 //____________________________________________________________________
55 AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
56     AliAnalysisTaskSE(name),
57     fDebug(0),
58     fOutputList(0),
59     fArray(0),
60     fZvtxDist(0)
61 {
62   // Default constructor
63   
64   DefineOutput(1, TList::Class());
65 }
66 //____________________________________________________________________
67 void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
68 {
69   // Create the output container
70   printf("AnalysisTaskFMD::CreateOutPutData() \n");
71   
72   fOutputList = new TList();//(TList*)GetOutputData(0);
73   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
74   
75   fArray     = new TObjArray();
76   fArray->SetName("FMD");
77   fArray->SetOwner();
78   TH1F* hEdist = 0;
79  
80   for(Int_t nEta = 0; nEta <= pars->GetNetaBins()+1; nEta++) {
81     TObjArray* etaArray = new TObjArray();
82     fArray->AddAtAndExpand(etaArray,nEta);
83     for(Int_t det =1; det<=3;det++)
84       {
85         TObjArray* detArray = new TObjArray();
86         detArray->SetName(Form("FMD%d",det));
87         etaArray->AddAtAndExpand(detArray,det);
88         Int_t nRings = (det==1 ? 1 : 2);
89         for(Int_t ring = 0;ring<nRings;ring++)
90           {
91             Char_t ringChar = (ring == 0 ? 'I' : 'O');
92             hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
93             hEdist->SetXTitle("#Delta E / E_{MIP}");
94             fOutputList->Add(hEdist);
95             detArray->AddAtAndExpand(hEdist,ring);
96           } 
97       }
98     
99   }
100   
101   fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
102   fZvtxDist->SetXTitle("z vertex");
103   fOutputList->Add(fZvtxDist);
104 }
105
106 //____________________________________________________________________
107 void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
108 {
109   
110   AliESDEvent* esd = (AliESDEvent*)InputEvent();
111   AliESD* old = esd->GetAliESDOld();
112   if (old) {
113     esd->CopyFromOldESD();
114   }
115   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
116   
117   Bool_t trigger = pars->IsEventTriggered(esd);
118   if(!trigger)
119     return;
120   Double_t vertex[3];
121   
122   pars->GetVertex(esd,vertex);
123   if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
124     return;
125   
126   fZvtxDist->Fill(vertex[2]);
127   
128   if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
129     return;
130   
131   AliESDFMD* fmd = esd->GetFMDData();
132   if (!fmd) return;
133   
134   for(UShort_t det=1;det<=3;det++) {
135       
136     Int_t nRings = (det==1 ? 1 : 2);
137     for (UShort_t ir = 0; ir < nRings; ir++) {
138   
139       Char_t   ring = (ir == 0 ? 'I' : 'O');
140       UShort_t nsec = (ir == 0 ? 20  : 40);
141       UShort_t nstr = (ir == 0 ? 512 : 256);
142       TH2F* hBg = pars->GetBackgroundCorrection(det,ring,0);
143       
144       for(UShort_t sec =0; sec < nsec;  sec++)  {
145         for(UShort_t strip = 0; strip < nstr; strip++) {
146           
147           
148           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
149           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
150           
151           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
152           
153           Int_t nEta = hBg->GetXaxis()->FindBin(eta);
154           
155           TObjArray* etaArray = (TObjArray*)fArray->At(nEta);
156           TObjArray* detArray = (TObjArray*)etaArray->At(det);
157           TH1F* Edist = (TH1F*)detArray->At(ir);
158           
159           Edist->Fill(mult);
160           
161         }
162       }
163     }
164   }
165   
166   PostData(1, fOutputList); 
167   
168 }
169
170 //____________________________________________________________________
171 void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {
172
173   //speciesOption:
174   //0: p+p Landau fit
175   //1: Pb+Pb triple landau convolution fit
176   
177   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
178   pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
179   
180   TFile fin(filename,"UPDATE");
181   
182   TList* list = (TList*)fin.Get("energyDist");
183   
184   AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
185   
186   EnergyDist->SetNetaBins(pars->GetNetaBins());
187   EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
188     
189   for(Int_t nEta = 1; nEta <= pars->GetNetaBins(); nEta++) {
190   
191     for(Int_t det = 1; det<=3; det++) {
192       Int_t nRings  =  (det == 1 ? 1 : 2);
193       for(Int_t ring = 0;ring<nRings; ring++) {
194         Char_t ringChar = (ring == 0 ? 'I' : 'O');
195         
196         TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
197         TF1* fitFunc = 0 ;
198         
199         if(hEdist->GetEntries() != 0) {
200           
201           hEdist->GetXaxis()->SetRangeUser(0.2,hEdist->GetXaxis()->GetXmax());
202           
203           if(speciesOption == 0)
204             fitFunc =  new TF1("FMDfitFunc","landau",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
205           if(speciesOption == 1) {
206             fitFunc = new TF1("FMDfitFunc",TripleLandau,hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,5,5);
207             fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
208             fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
209             fitFunc->SetParLimits(1,0.6,1.2);
210             fitFunc->SetParLimits(3,0,1);
211             fitFunc->SetParLimits(4,0,1);
212             
213           }
214           
215           
216           hEdist->Fit(fitFunc,"","",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
217           fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
218           
219         }
220         
221         TH2F* hBg = pars->GetBackgroundCorrection(det,ringChar,0);
222         EnergyDist->SetEnergyDistribution(det,ringChar,hBg->GetXaxis()->GetBinCenter(nEta),hEdist);
223       }
224       
225     }
226     
227   }
228   
229   fin.Write();
230   fin.Close();
231   
232   if(store) {
233     TFile fcalib("$ALICE_ROOT/FMD/Correction/EnergyDistribution/energydistributions.root","RECREATE");
234     EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
235     fcalib.Close();
236   }
237
238
239
240 }
241
242 //____________________________________________________________________
243 //
244 // EOF
245 //