double counting for Nsigma
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 //-----------------------------------------------------------------\r
17 //         AliAnalysisTaskSpectraAllChAOD class\r
18 //-----------------------------------------------------------------\r
19 \r
20 #include "TChain.h"\r
21 #include "TTree.h"\r
22 #include "TLegend.h"\r
23 #include "TH1F.h"\r
24 #include "TH2F.h"\r
25 #include "THnSparse.h"\r
26 #include "TCanvas.h"\r
27 #include "AliAnalysisTask.h"\r
28 #include "AliAODTrack.h"\r
29 #include "AliAODMCParticle.h"\r
30 #include "AliVParticle.h"\r
31 #include "AliAODEvent.h"\r
32 #include "AliAODInputHandler.h"\r
33 #include "AliAnalysisTaskSpectraAllChAOD.h"\r
34 #include "AliAnalysisTaskESDfilter.h"\r
35 #include "AliAnalysisDataContainer.h"\r
36 #include "AliSpectraAODTrackCuts.h"\r
37 #include "AliSpectraAODEventCuts.h"\r
38 #include "AliHelperPID.h"\r
39 #include "AliPIDCombined.h"\r
40 #include "AliCentrality.h"\r
41 #include "TProof.h"\r
42 #include "AliVEvent.h"\r
43 #include "AliStack.h"\r
44 #include <TMCProcess.h>\r
45 \r
46 #include <iostream>\r
47 \r
48 using namespace AliHelperPIDNameSpace;\r
49 using namespace std;\r
50 \r
51 ClassImp(AliAnalysisTaskSpectraAllChAOD)\r
52 \r
53 //________________________________________________________________________\r
54 AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  \r
55   fAOD(0x0),\r
56   fTrackCuts(0x0),\r
57   fEventCuts(0x0),\r
58   fHelperPID(0x0),\r
59   fIsMC(0),\r
60   fCharge(0),\r
61   fVZEROside(0),\r
62   fOutput(0x0),\r
63   fnCentBins(20),\r
64   fnQvecBins(40)\r
65 {\r
66   // Default constructor\r
67   DefineInput(0, TChain::Class());\r
68   DefineOutput(1, TList::Class());\r
69   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
70   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
71   DefineOutput(4, AliHelperPID::Class());\r
72 }\r
73 \r
74 //________________________________________________________________________\r
75 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()\r
76 {\r
77   // create output objects\r
78   fOutput=new TList();\r
79   fOutput->SetOwner();\r
80   fOutput->SetName("fOutput");\r
81   \r
82   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
83   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
84   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");\r
85   \r
86   // binning common to all the THn\r
87   const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.};\r
88   const Int_t nptBins=26;\r
89   \r
90   //dimensions of THnSparse for tracks\r
91   const Int_t nvartrk=8;\r
92   //                                             pt          cent        Q vec     IDrec     IDgen       isph           iswd      y\r
93   Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,   fnQvecBins,        4,        3,         2,          2,       2};\r
94   Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,            0.,     -1.5,      -0.5,      -0.5,        -0.5,   -0.5};\r
95   Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,            8.,      2.5,      2.5,       1.5,         1.5,     0.5};    \r
96   THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);\r
97   NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");\r
98   NSparseHistTrk->GetAxis(0)->SetName("pT_rec");\r
99   NSparseHistTrk->SetBinEdges(0,ptBins);\r
100   NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
101   NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
102   NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");\r
103   NSparseHistTrk->GetAxis(2)->SetName("Q_vec");\r
104   NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");\r
105   NSparseHistTrk->GetAxis(3)->SetName("ID_rec");\r
106   NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");\r
107   NSparseHistTrk->GetAxis(4)->SetName("ID_gen");\r
108   NSparseHistTrk->GetAxis(5)->SetTitle("isph");\r
109   NSparseHistTrk->GetAxis(5)->SetName("isph");\r
110   NSparseHistTrk->GetAxis(6)->SetTitle("iswd");\r
111   NSparseHistTrk->GetAxis(6)->SetName("iswd");\r
112   NSparseHistTrk->GetAxis(7)->SetTitle("y");\r
113   NSparseHistTrk->GetAxis(7)->SetName("y");\r
114   fOutput->Add(NSparseHistTrk);\r
115   \r
116   //dimensions of THnSparse for stack\r
117   const Int_t nvarst=5;\r
118   //                                             pt          cent    IDgen        isph        y\r
119   Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2};\r
120   Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5};\r
121   Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5};\r
122   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);\r
123   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");\r
124   NSparseHistSt->SetBinEdges(0,ptBins);\r
125   NSparseHistSt->GetAxis(0)->SetName("pT_rec");\r
126   NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
127   NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
128   NSparseHistSt->GetAxis(2)->SetTitle("ID gen");\r
129   NSparseHistSt->GetAxis(2)->SetName("ID_gen");\r
130   NSparseHistSt->GetAxis(3)->SetTitle("isph");\r
131   NSparseHistSt->GetAxis(3)->SetName("isph");\r
132   NSparseHistSt->GetAxis(4)->SetTitle("y");\r
133   NSparseHistSt->GetAxis(4)->SetName("y");\r
134   fOutput->Add(NSparseHistSt);\r
135   \r
136   //dimensions of THnSparse for the normalization\r
137   const Int_t nvarev=2;\r
138   //                                             cent             Q vec   \r
139   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins};\r
140   Double_t xminHistRealEv[nvarev] = {           0.,               0.};\r
141   Double_t xmaxHistRealEv[nvarev] = {       100.,              10.};\r
142   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);\r
143   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
144   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
145   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");\r
146   NSparseHistEv->GetAxis(1)->SetName("Q_vec");\r
147   fOutput->Add(NSparseHistEv);\r
148   \r
149   PostData(1, fOutput  );\r
150   PostData(2, fEventCuts);\r
151   PostData(3, fTrackCuts);\r
152   PostData(4, fHelperPID);\r
153 }\r
154 \r
155 //________________________________________________________________________\r
156 \r
157 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)\r
158 {\r
159   //Printf("An event");\r
160   // main event loop\r
161   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
162   if (!fAOD) {\r
163     AliWarning("ERROR: AliAODEvent not available \n");\r
164     return;\r
165   }\r
166   \r
167   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
168     {\r
169       AliFatal("Not processing AODs");\r
170     }\r
171   \r
172   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
173 \r
174   //Default TPC priors\r
175   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?\r
176   \r
177   Double_t Qvec=0.;//in case of MC we save space in the memory\r
178   if(!fIsMC){\r
179     if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();\r
180     else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();\r
181   }\r
182 \r
183   Double_t Cent=fEventCuts->GetCent();\r
184     \r
185   // First do MC to fill up the MC particle array\r
186   TClonesArray *arrayMC = 0;\r
187   if (fIsMC)\r
188     {\r
189       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
190       if (!arrayMC) {\r
191         AliFatal("Error: MC particles branch not found!\n");\r
192       }\r
193       Int_t nMC = arrayMC->GetEntries();\r
194       for (Int_t iMC = 0; iMC < nMC; iMC++)\r
195         {\r
196           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
197           if(!partMC->Charge()) continue;//Skip neutrals\r
198           if(fCharge != 0 && partMC->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge\r
199           \r
200           if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!\r
201           \r
202           //pt     cent    IDgen        isph        y\r
203           Double_t varSt[5];\r
204           varSt[0]=partMC->Pt();\r
205           varSt[1]=Cent;\r
206           varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);\r
207           varSt[3]=(Double_t)partMC->IsPhysicalPrimary();\r
208           varSt[4]=partMC->Y();\r
209           ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop\r
210           \r
211           //Printf("a particle");\r
212           \r
213         }\r
214     }\r
215   \r
216   //main loop on tracks\r
217   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {\r
218     AliAODTrack* track = fAOD->GetTrack(iTracks);\r
219     if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge \r
220     if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)\r
221     Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector\r
222     Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));\r
223     Int_t IDgen=kSpUndefined;//set if MC\r
224     Int_t isph=-999;\r
225     Int_t iswd=-999;\r
226     \r
227     if (arrayMC) {\r
228       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
229       if (!partMC) { \r
230         AliError("Cannot get MC particle");\r
231         continue; \r
232       }\r
233       IDgen=fHelperPID->GetParticleSpecies(partMC);\r
234       isph=partMC->IsPhysicalPrimary();\r
235       iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions\r
236     }\r
237     \r
238     //pt     cent    Q vec     IDrec     IDgen       isph           iswd      y\r
239     Double_t varTrk[8];\r
240     varTrk[0]=track->Pt();\r
241     varTrk[1]=Cent;\r
242     varTrk[2]=Qvec;\r
243     varTrk[3]=(Double_t)IDrec;\r
244     varTrk[4]=(Double_t)IDgen;\r
245     varTrk[5]=(Double_t)isph;\r
246     varTrk[6]=(Double_t)iswd;\r
247     varTrk[7]=y;\r
248     ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
249     \r
250     //for nsigma PID fill double counting of ID\r
251     if(fHelperPID->GetPIDType()<kBayes){//only nsigma\r
252       Bool_t *HasDC;\r
253       HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting\r
254       for(Int_t ipart=0;ipart<kNSpecies;ipart++){\r
255         if(HasDC[ipart]==kTRUE){\r
256           varTrk[3]=(Double_t)ipart;\r
257           ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
258         }\r
259       }\r
260     }\r
261     \r
262     //fill all charged (-1)\r
263     varTrk[3]=-1.;\r
264     ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
265     \r
266     //Printf("a track");\r
267     \r
268   } // end loop on tracks\r
269   \r
270   Double_t varEv[2];\r
271   varEv[0]=Cent;\r
272   varEv[1]=Qvec;\r
273   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop\r
274   \r
275   PostData(1, fOutput  );\r
276   PostData(2, fEventCuts);\r
277   PostData(3, fTrackCuts);\r
278   PostData(4, fHelperPID);\r
279 }\r
280 \r
281 //_________________________________________________________________\r
282 void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)\r
283 {\r
284   // Terminate\r
285 }\r