Port of changes from v4-07-Release and additional rule conformance
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoShareQualityCorrFctn.cxx
CommitLineData
0215f606 1////////////////////////////////////////////////////////////////////////////////
2/// ///
3/// AliFemtoShareQualityCorrFctn - A correlation function that saves the ///
4/// amount of sharing and splitting hits per pair as a function of qinv ///
5/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
6/// ///
7////////////////////////////////////////////////////////////////////////////////
67427ff7 8
b2f60a91 9#include "AliFemtoShareQualityCorrFctn.h"
65423af9 10//#include "AliFemtoHisto.hh"
67427ff7 11#include <cstdio>
12
13#ifdef __ROOT__
14ClassImp(AliFemtoShareQualityCorrFctn)
15#endif
16
17//____________________________
0215f606 18AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
19 fShareNumerator(0),
20 fShareDenominator(0),
21 fQualityNumerator(0),
22 fQualityDenominator(0)
23{
67427ff7 24 // set up numerator
25 // title = "Num Qinv (MeV/c)";
d92ed900 26 char tTitNum[100] = "NumShare";
27 strcat(tTitNum,title);
cc5faabc 28 fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
67427ff7 29 // set up denominator
30 //title = "Den Qinv (MeV/c)";
d92ed900 31 char tTitDen[100] = "DenShare";
32 strcat(tTitDen,title);
cc5faabc 33 fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
67427ff7 34
d92ed900 35 char tTit2Num[100] = "NumQuality";
36 strcat(tTit2Num,title);
cc5faabc 37 fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
67427ff7 38 // set up denominator
39 //title = "Den Qinv (MeV/c)";
d92ed900 40 char tTit2Den[100] = "DenQuality";
41 strcat(tTit2Den,title);
cc5faabc 42 fQualityDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
67427ff7 43 // set up ratio
44 //title = "Ratio Qinv (MeV/c)";
45 // this next bit is unfortunately needed so that we can have many histos of same "title"
46 // it is neccessary if we typedef TH2D to TH1d (which we do)
47 //mShareNumerator->SetDirectory(0);
48 //mShareDenominator->SetDirectory(0);
49 //mRatio->SetDirectory(0);
50
51 // to enable error bar calculation...
52 fShareNumerator->Sumw2();
53 fShareDenominator->Sumw2();
54
55 fQualityNumerator->Sumw2();
56 fQualityDenominator->Sumw2();
57}
58
59//____________________________
0215f606 60AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
61 fShareNumerator(0),
62 fShareDenominator(0),
63 fQualityNumerator(0),
64 fQualityDenominator(0)
65{
d92ed900 66 // copy constructor
0215f606 67 if (aCorrFctn.fShareNumerator)
68 fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
69 if (aCorrFctn.fShareDenominator)
70 fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
71 if (aCorrFctn.fQualityNumerator)
72 fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
73 if (aCorrFctn.fQualityDenominator)
74 fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
75}
76//____________________________
67427ff7 77AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
d92ed900 78 // destructor
67427ff7 79 delete fShareNumerator;
80 delete fShareDenominator;
81 delete fQualityNumerator;
82 delete fQualityDenominator;
83}
84//_________________________
0215f606 85AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
86{
d92ed900 87 // assignment operator
0215f606 88 if (this == &aCorrFctn)
89 return *this;
90
91 if (aCorrFctn.fShareNumerator)
92 fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
93 else
94 fShareNumerator = 0;
95 if (aCorrFctn.fShareDenominator)
96 fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
97 else
98 fShareDenominator = 0;
99 if (aCorrFctn.fQualityNumerator)
100 fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
101 else
102 fQualityNumerator = 0;
103 if (aCorrFctn.fQualityDenominator)
104 fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
105 else
106 fQualityDenominator = 0;
107
108 return *this;
109}
110//_________________________
67427ff7 111void AliFemtoShareQualityCorrFctn::Finish(){
112 // here is where we should normalize, fit, etc...
113 // we should NOT Draw() the histos (as I had done it below),
114 // since we want to insulate ourselves from root at this level
115 // of the code. Do it instead at root command line with browser.
116 // mShareNumerator->Draw();
117 //mShareDenominator->Draw();
118 //mRatio->Draw();
119
120}
121
122//____________________________
123AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
d92ed900 124 // create report
67427ff7 125 string stemp = "Qinv Correlation Function Report:\n";
126 char ctemp[100];
127 sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
128 stemp += ctemp;
129 sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
130 stemp += ctemp;
131 // stemp += mCoulombWeight->Report();
132 AliFemtoString returnThis = stemp;
133 return returnThis;
134}
135//____________________________
65423af9 136void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
d92ed900 137 // add real (effect) pair
cc5faabc 138 // double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
67427ff7 139 Int_t nh = 0;
140 Int_t an = 0;
141 Int_t ns = 0;
142
65423af9 143 for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
67427ff7 144 // If both have clusters in the same row
65423af9 145 if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
146 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 147 // Do they share it ?
65423af9 148 if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
149 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
67427ff7 150 {
cc5faabc 151// if (tQinv < 0.01) {
152// cout << "Shared cluster in row " << imap << endl;
153// }
67427ff7 154 an++;
155 nh+=2;
156 ns+=2;
157 }
158
159 // Different hits on the same padrow
160 else {
161 an--;
162 nh+=2;
163 }
164 }
65423af9 165 else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
166 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 167 // One track has a hit, the other does not
168 an++;
169 nh++;
170 }
171 }
cc5faabc 172// if (tQinv < 0.01) {
173// cout << "Qinv of the pair is " << tQinv << endl;
174// cout << "Clusters: " << endl;
175// for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
176// cout << imap ;
177// if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
178// else cout << " 0 " ;
179// if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
180// else cout << " 0 " ;
181// cout << " ";
182// if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
183// else cout << " X ";
184// if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
185// else cout << " X ";
186// cout << endl;
187// }
188// }
67427ff7 189
190 Float_t hsmval = 0.0;
191 Float_t hsfval = 0.0;
192
193 if (nh >0) {
194 hsmval = an*1.0/nh;
195 hsfval = ns*1.0/nh;
196 }
197
cc5faabc 198// if ((tQinv < 0.01) && (hsmval<-0.46)) {
199// cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
200// cout << "Qinv of the pair is " << tQinv << endl;
201// cout << "Clusters: " << endl;
202// for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
203// cout << imap ;
204// if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
205// else cout << " 0 " ;
206// if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
207// else cout << " 0 " ;
208// cout << " ";
209// if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
210// else cout << " X ";
211// if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
212// else cout << " X ";
213// cout << endl;
214// }
215// cout << "Momentum1 "
216// << pair->Track1()->Track()->P().x() << " "
217// << pair->Track1()->Track()->P().y() << " "
218// << pair->Track1()->Track()->P().z() << " "
219// << pair->Track1()->Track()->Label() << " "
220// << pair->Track1()->Track()->TrackId() << " "
221// << pair->Track1()->Track()->Flags() << " "
222// << pair->Track1()->Track()->KinkIndex(0) << " "
223// << pair->Track1()->Track()->KinkIndex(1) << " "
224// << pair->Track1()->Track()->KinkIndex(2) << " "
225// << endl;
226// cout << "Momentum2 "
227// << pair->Track2()->Track()->P().x() << " "
228// << pair->Track2()->Track()->P().y() << " "
229// << pair->Track2()->Track()->P().z() << " "
230// << pair->Track2()->Track()->Label() << " "
231// << pair->Track2()->Track()->TrackId() << " "
232// << pair->Track2()->Track()->Flags() << " "
233// << pair->Track2()->Track()->KinkIndex(0) << " "
234// << pair->Track2()->Track()->KinkIndex(1) << " "
235// << pair->Track2()->Track()->KinkIndex(2) << " "
236// << endl;
237// }
67427ff7 238
cc5faabc 239// fShareNumerator->Fill(tQinv, hsfval);
240// fQualityNumerator->Fill(tQinv, hsmval);
d92ed900 241 // cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
65423af9 242 //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
67427ff7 243}
244//____________________________
65423af9 245void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
d92ed900 246 // add mixed (background) pair
67427ff7 247 double weight = 1.0;
d92ed900 248 double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
67427ff7 249 Int_t nh = 0;
250 Int_t an = 0;
251 Int_t ns = 0;
252
65423af9 253 for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
67427ff7 254 // If both have clusters in the same row
65423af9 255 if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
256 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 257 // Do they share it ?
cc5faabc 258 if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
65423af9 259 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
67427ff7 260 {
261 // cout << "A shared cluster !!!" << endl;
262 // cout << "imap idx1 idx2 "
263 // << imap << " "
264 // << tP1idx[imap] << " " << tP2idx[imap] << endl;
265 an++;
266 nh+=2;
267 ns+=2;
268 }
269
270 // Different hits on the same padrow
271 else {
272 an--;
273 nh+=2;
274 }
275 }
65423af9 276 else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
277 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 278 // One track has a hit, the other does not
279 an++;
280 nh++;
281 }
282 }
283
284 Float_t hsmval = 0.0;
285 Float_t hsfval = 0.0;
286
287 if (nh >0) {
288 hsmval = an*1.0/nh;
289 hsfval = ns*1.0/nh;
290 }
291
d92ed900 292 fShareDenominator->Fill(tQinv,hsfval,weight);
293 fQualityDenominator->Fill(tQinv,hsmval,weight);
67427ff7 294}
295
296
297void AliFemtoShareQualityCorrFctn::WriteHistos()
298{
299 fShareNumerator->Write();
300 fShareDenominator->Write();
301 fQualityNumerator->Write();
302 fQualityDenominator->Write();
303
304}