edd93a102f9c80e2b233ad6ecec0950b85c89ae3
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoCutMonitorParticlePtPDG.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 //                                                                            //
3 // AliFemtoCutMonitorParticlePtPDG - the cut monitor for particles to study    //
4 // the difference between reconstructed and true momentum                     //
5 //                                                                            //
6 ////////////////////////////////////////////////////////////////////////////////
7 #include "AliFemtoCutMonitorParticlePtPDG.h"
8 #include "AliFemtoModelHiddenInfo.h"
9 #include <TH1D.h>
10 #include <TH2D.h>
11 #include <TList.h>
12 #include <TMath.h>
13
14 AliFemtoCutMonitorParticlePtPDG::AliFemtoCutMonitorParticlePtPDG():
15   fPtPDG(0),ftpcHist(0),fPtGoodPi(0),fPtFakePi(0),fPtGoodK(0),fPtFakeK(0),
16   fPtGoodP(0),fPtFakeP(0),fPtRPi(0),fPtRK(0),fPtRP(0),
17   fPtContP(0),
18   fPtContPi(0),
19   fPtContMup(0),
20   fPtContElp(0),
21   fMass(0.13957)
22 {
23   // Default constructor
24   fPtPDG = new TH2D("PtPDG", "PDG vs Pt",  5, .0, 5.0, 100, 0.1, 2.0);
25   ftpcHist=new TH2D("tpcHist","TPC dE/dX vs momentum",100,0.1,2.7,100,0.,6.);
26   fPtGoodPi = new TH1D("PtGoodPi", "good pions Pt",              100, 0.1, 2.0);
27   fPtRPi = new TH1D("PtRPi", "right pdg pions Pt",               100, 0.1, 2.0);
28   fPtFakePi = new TH1D("PtFakePi", "fake pions Pt",              100, 0.1, 2.0);
29   fPtGoodK = new TH1D("PtGoodK", "good kaons Pt",                100, 0.1, 2.0);
30   fPtRK = new TH1D("PtRK", "right pdg kaons Pt",                 100, 0.1, 2.0);
31   fPtFakeK = new TH1D("PtFakeK", "fake kaons Pt",                100, 0.1, 2.0); 
32   fPtGoodP = new TH1D("PtGoodP", "good protons Pt",              100, 0.1, 2.0);
33   fPtFakeP = new TH1D("PtFakeP", "fake protons Pt",              100, 0.1, 2.0); 
34   fPtRP = new TH1D("PtRP", "right pdg protons Pt",               100, 0.1, 2.0);   
35   
36   fPtContP = new TH1D("PtContP", "contamination Pt",               100, 0.1, 2.0);   
37   fPtContPi = new TH1D("PtContPi", "contamination Pt",               100, 0.1, 2.0);   
38   fPtContMup = new TH1D("PtContMup", "contamination Pt",               100, 0.1, 2.0);   
39   fPtContElp = new TH1D("PtContElp", "contamination Pt",               100, 0.1, 2.0);   
40
41 }
42
43 AliFemtoCutMonitorParticlePtPDG::AliFemtoCutMonitorParticlePtPDG(const char *aName, float aMass):
44   AliFemtoCutMonitor(),
45   fPtPDG(0),ftpcHist(0),fPtGoodPi(0),fPtFakePi(0),fPtGoodK(0),fPtFakeK(0),
46   fPtGoodP(0),fPtFakeP(0),fPtRPi(0),fPtRK(0),fPtRP(0),
47   fPtContP(0),
48   fPtContPi(0),
49   fPtContMup(0),
50   fPtContElp(0),
51   fMass(aMass)
52 {
53   // Normal constructor
54   char name[200];
55   snprintf(name, 200, "PtPDG%s", aName);
56   fPtPDG = new TH2D(name, "PDG vs Pt", 10, 0.0, 5.0, 100, 0.1, 2.0);
57   snprintf(name, 200, "tpcHist%s", aName);
58   ftpcHist=new TH2D(name,"TPC dE/dX vs momentum",100,0.1,2.7,100,0.,6.);
59   snprintf(name, 200, "PtGoodPi%s", aName);
60   fPtGoodPi = new TH1D(name, "good pions Pt",                    100, 0.1, 2.0);
61   snprintf(name, 200, "PtFakePi%s", aName);
62   fPtFakePi = new TH1D(name, "fake pions Pt",              100, 0.1, 2.0);
63   snprintf(name, 200, "PtRPi%s", aName);
64   fPtRPi = new TH1D(name, "right pdg pions Pt",               100, 0.1, 2.0);
65   snprintf(name, 200, "PtGoodK%s", aName);
66   fPtGoodK = new TH1D(name, "good kaons Pt",                     100, 0.1, 2.0);
67   snprintf(name, 200, "PtFakeK%s", aName);
68   fPtFakeK = new TH1D(name, "fake kaons Pt",                100, 0.1, 2.0);
69   snprintf(name, 200, "PtRK%s", aName);
70   fPtRK = new TH1D(name, "right pdg kaons Pt",                 100, 0.1, 2.0);  
71    snprintf(name, 200, "PtGoodP%s", aName);
72   fPtGoodP = new TH1D(name, "good protons Pt",                     100, 0.1, 2.0);
73   snprintf(name, 200, "PtFakeP%s", aName);
74   fPtFakeP = new TH1D(name, "fake protons Pt",                100, 0.1, 2.0);
75   snprintf(name, 200, "PtRP%s", aName);
76   fPtRP = new TH1D(name, "right pdg protons Pt",                 100, 0.1, 2.0);   
77
78   snprintf(name, 200, "PtContP%s", aName);
79   fPtContP = new TH1D(name, "contamination",                 100, 0.1, 2.0);   
80   snprintf(name, 200, "PtContPi%s", aName);
81   fPtContPi = new TH1D(name, "contamination",                 100, 0.1, 2.0);   
82   snprintf(name, 200, "PtContMup%s", aName);
83   fPtContMup = new TH1D(name, "contamination",                 100, 0.1, 2.0);   
84   snprintf(name, 200, "PtContElp%s", aName);
85   fPtContElp = new TH1D(name, "contamination",                 100, 0.1, 2.0);   
86  
87 }
88
89 AliFemtoCutMonitorParticlePtPDG::AliFemtoCutMonitorParticlePtPDG(const AliFemtoCutMonitorParticlePtPDG &aCut):
90   AliFemtoCutMonitor(),
91   fPtPDG(0),ftpcHist(0),fPtGoodPi(0),fPtFakePi(0),fPtGoodK(0),fPtFakeK(0),
92   fPtGoodP(0),fPtFakeP(0),fPtRPi(0),fPtRK(0),fPtRP(0),
93   fPtContP(0),
94   fPtContPi(0),
95   fPtContMup(0),
96   fPtContElp(0), 
97   fMass(0.13957)
98 {
99   // copy constructor
100   if (fPtPDG) delete fPtPDG;
101   fPtPDG = new TH2D(*aCut.fPtPDG);
102   ftpcHist= new TH2D(*aCut.ftpcHist);
103   fPtGoodPi= new TH1D(*aCut.fPtGoodPi);
104   fPtFakePi= new TH1D(*aCut.fPtFakePi);
105   fPtGoodK= new TH1D(*aCut.fPtGoodK);
106   fPtFakeK= new TH1D(*aCut.fPtFakePi);
107   fPtGoodP= new TH1D(*aCut.fPtGoodP);
108   fPtFakeP= new TH1D(*aCut.fPtFakePi);
109   fPtRPi= new TH1D(*aCut.fPtRPi);
110   fPtRK= new TH1D(*aCut.fPtRK);
111   fPtRP= new TH1D(*aCut.fPtRP);  
112   
113   fPtContP= new TH1D(*aCut.fPtContP);
114   fPtContPi= new TH1D(*aCut.fPtContPi);
115   fPtContMup= new TH1D(*aCut.fPtContMup);
116   fPtContElp= new TH1D(*aCut.fPtContElp);
117   
118   fMass = aCut.fMass; 
119 }
120
121 AliFemtoCutMonitorParticlePtPDG::~AliFemtoCutMonitorParticlePtPDG()
122 {
123   // Destructor
124   delete fPtPDG;
125   delete ftpcHist;
126   delete fPtGoodPi;
127   delete fPtFakePi;
128   delete fPtGoodK;
129   delete fPtFakeK;
130   delete fPtGoodP;
131   delete fPtFakeP; 
132   delete fPtRPi;
133   delete fPtRK;
134   delete fPtRP;
135
136   delete fPtContP;
137   delete fPtContPi;
138   delete fPtContMup;
139   delete fPtContElp;
140
141 }
142
143 AliFemtoCutMonitorParticlePtPDG& AliFemtoCutMonitorParticlePtPDG::operator=(const AliFemtoCutMonitorParticlePtPDG& aCut)
144 {
145   // assignment operator
146   if (this == &aCut) 
147     return *this;
148
149   if (fPtPDG) delete fPtPDG;
150   fPtPDG = new TH2D(*aCut.fPtPDG);
151   
152   if (ftpcHist) delete ftpcHist;
153   ftpcHist = new TH2D(*aCut.ftpcHist);
154   
155    if (fPtGoodPi) delete fPtGoodPi;
156   fPtGoodPi = new TH1D(*aCut.fPtGoodPi);
157   
158    if (fPtFakePi) delete fPtFakePi;
159   fPtFakePi = new TH1D(*aCut.fPtFakePi);
160   
161    if (fPtRPi) delete fPtRPi;
162   fPtRPi = new TH1D(*aCut.fPtRPi);
163   
164   if (fPtGoodK) delete fPtGoodK;
165   fPtGoodK = new TH1D(*aCut.fPtGoodK);
166   
167    if (fPtFakeK) delete fPtFakeK;
168   fPtFakeK = new TH1D(*aCut.fPtFakeK);
169   
170    if (fPtRK) delete fPtRK;
171   fPtRK = new TH1D(*aCut.fPtRK);  
172    
173   if (fPtGoodP) delete fPtGoodP;
174   fPtGoodP = new TH1D(*aCut.fPtGoodP);
175   
176    if (fPtFakeP) delete fPtFakeP;
177   fPtFakeP = new TH1D(*aCut.fPtFakeP);
178   
179    if (fPtRP) delete fPtRP;
180   fPtRP = new TH1D(*aCut.fPtRP);
181  
182    if (fPtContP) delete fPtContP;
183   fPtContP = new TH1D(*aCut.fPtContP);
184
185    if (fPtContPi) delete fPtContPi;
186   fPtContPi = new TH1D(*aCut.fPtContPi);
187   
188      if (fPtContMup) delete fPtContMup;
189   fPtContMup = new TH1D(*aCut.fPtContMup);
190
191      if (fPtContElp) delete fPtContElp;
192   fPtContElp = new TH1D(*aCut.fPtContElp);
193   
194   return *this;
195 }
196
197 AliFemtoString AliFemtoCutMonitorParticlePtPDG::Report(){ 
198   // Prepare report from the execution
199   string stemp = "*** AliFemtoCutMonitorParticlePtPDG report"; 
200   AliFemtoString returnThis = stemp;
201   return returnThis; 
202 }
203
204 void AliFemtoCutMonitorParticlePtPDG::Fill(const AliFemtoTrack* aTrack)
205 {
206   // Fill in the monitor histograms with the values from the current track
207   //  float tEnergy = ::sqrt(aTrack->P().mag2()+fMass*fMass);
208   //  float tRapidity = 0.5*::log((tEnergy+aTrack->P().z())/(tEnergy-aTrack->P().z()));
209   float tPt = ::sqrt((aTrack->P().x())*(aTrack->P().x())+(aTrack->P().y())*(aTrack->P().y()));
210   //  float tEta = -TMath::Log(TMath::Tan(aTrack->P().theta()/2.0));
211   //  float tPhi = aTrack->P().phi();
212   float tP = ::sqrt((aTrack->P().z())*(aTrack->P().z())+(aTrack->P().x())*(aTrack->P().x())+(aTrack->P().y())*(aTrack->P().y()));;
213   float dedx = aTrack->TPCsignalN();
214   float w[10];
215   w[0] = aTrack->PidProbElectron();
216   w[1] = aTrack->PidProbMuon();
217   w[2] = aTrack->PidProbPion();
218   w[3] = aTrack->PidProbKaon(); 
219   w[4] = aTrack->PidProbProton();
220   
221    Int_t pdg1=0;
222   AliFemtoModelHiddenInfo *info = ( AliFemtoModelHiddenInfo *) aTrack->GetHiddenInfo();
223   if(info)pdg1 = info->GetPDGPid();
224
225
226 //most probable particle  
227   fPtGoodPi->Fill(tPt);
228   fPtGoodP->Fill(tPt);
229   fPtGoodK->Fill(tPt);
230
231 //contaminations 
232   if (abs(pdg1)!=321)fPtFakeK->Fill(tPt);
233   if (abs(pdg1)!=211)fPtFakePi->Fill(tPt);
234   if (abs(pdg1)!=2212)fPtFakeP->Fill(tPt);
235
236                
237 //contaminations for kaons 
238   if (abs(pdg1)==2212)fPtContP->Fill(tPt);
239   if (abs(pdg1)==211)fPtContPi->Fill(tPt);
240   if (abs(pdg1)==13)fPtContMup->Fill(tPt);
241   if (abs(pdg1)==11)fPtContElp->Fill(tPt);
242              
243   Float_t pdg=-1.0;
244   if(abs(pdg1)==211){
245    pdg=2.0;
246    fPtRPi->Fill(tPt);
247   }
248
249   if(abs(pdg1)==321){
250    pdg=3.0;
251    fPtRK->Fill(tPt);
252    }
253   
254   if(abs(pdg1)==2212){
255    pdg=4.0;
256    fPtRP->Fill(tPt);
257    }
258   
259   if(abs(pdg1)==11)pdg=0.0; //+-electron
260   if(abs(pdg1)==13)pdg=1.0;  //+-muon
261   
262   //cout<<"pdg from CutMonitor.."<<pdg1<<"pdg"<<pdg<<endl;
263  
264    
265   fPtPDG->Fill(pdg, tPt);
266   ftpcHist->Fill(tP,dedx);
267   
268  
269 }
270
271 void AliFemtoCutMonitorParticlePtPDG::Write()
272 {
273   // Write out the relevant histograms
274   
275   fPtPDG->Write();
276   ftpcHist->Write();
277   fPtGoodPi->Write();
278   fPtFakePi->Write();
279   fPtGoodK->Write();
280   fPtFakeK->Write();
281   fPtGoodP->Write();
282   fPtFakeP->Write(); 
283   fPtRPi->Write();
284   fPtRK->Write();
285   fPtRP->Write();
286   fPtContP->Write();
287   fPtContPi->Write();
288   fPtContMup->Write();
289   fPtContElp->Write();
290 }
291
292 TList *AliFemtoCutMonitorParticlePtPDG::GetOutputList()
293 {
294   TList *tOutputList = new TList();
295   tOutputList->Add(fPtPDG);
296   tOutputList->Add(ftpcHist);
297   tOutputList->Add(fPtGoodPi);
298   tOutputList->Add(fPtFakePi);
299   tOutputList->Add(fPtGoodK);
300   tOutputList->Add(fPtFakeK);
301   tOutputList->Add(fPtGoodP);
302   tOutputList->Add(fPtFakeP); 
303   tOutputList->Add(fPtRPi);
304   tOutputList->Add(fPtRK);
305   tOutputList->Add(fPtRP);
306  
307   tOutputList->Add(fPtContP);
308   tOutputList->Add(fPtContPi);
309   tOutputList->Add(fPtContMup);
310   tOutputList->Add(fPtContElp); 
311
312   return tOutputList;
313 }