Coding convention: RN11 violation -> suppression
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / CorrFctn / AliFemtoShareQualityCorrFctn.cxx
CommitLineData
67427ff7 1/***************************************************************************
2 *
3 * $Id$
4 *
5 * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6 ***************************************************************************
7 *
8 * Description: part of STAR HBT Framework: AliFemtoMaker package
9 * a simple Q-invariant correlation function
10 *
11 ***************************************************************************
12 *
13 * $Log$
b2f60a91 14 * Revision 1.1.1.1 2007/04/25 15:38:41 panos
15 * Importing the HBT code dir
16 *
67427ff7 17 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
18 * First version on CVS
19 *
20 * Revision 1.4 2000/01/25 17:34:45 laue
21 * I. In order to run the stand alone version of the AliFemtoMaker the following
22 * changes have been done:
23 * a) all ClassDefs and ClassImps have been put into #ifdef __ROOT__ statements
24 * b) unnecessary includes of StMaker.h have been removed
25 * c) the subdirectory AliFemtoMaker/doc/Make has been created including everything
26 * needed for the stand alone version
27 *
28 * II. To reduce the amount of compiler warning
29 * a) some variables have been type casted
30 * b) some destructors have been declared as virtual
31 *
32 * Revision 1.3 1999/07/29 02:47:09 lisa
33 * 1) add OpeningAngle correlation function 2) add AliFemtoMcEventReader 3) make histos in CorrFctns do errors correctly
34 *
35 * Revision 1.2 1999/07/06 22:33:20 lisa
36 * Adjusted all to work in pro and new - dev itself is broken
37 *
38 * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
39 * Installation of AliFemtoMaker
40 *
41 **************************************************************************/
42
b2f60a91 43#include "AliFemtoShareQualityCorrFctn.h"
67427ff7 44//#include "Infrastructure/AliFemtoHisto.hh"
45#include <cstdio>
46
47#ifdef __ROOT__
48ClassImp(AliFemtoShareQualityCorrFctn)
49#endif
50
51//____________________________
52AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi){
53 // set up numerator
54 // title = "Num Qinv (MeV/c)";
55 char TitNum[100] = "NumShare";
56 strcat(TitNum,title);
57 fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
58 // set up denominator
59 //title = "Den Qinv (MeV/c)";
60 char TitDen[100] = "DenShare";
61 strcat(TitDen,title);
62 fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
63
64 char Tit2Num[100] = "NumQuality";
65 strcat(Tit2Num,title);
66 fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
67 // set up denominator
68 //title = "Den Qinv (MeV/c)";
69 char Tit2Den[100] = "DenQuality";
70 strcat(Tit2Den,title);
71 fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
72 // set up ratio
73 //title = "Ratio Qinv (MeV/c)";
74 // this next bit is unfortunately needed so that we can have many histos of same "title"
75 // it is neccessary if we typedef TH2D to TH1d (which we do)
76 //mShareNumerator->SetDirectory(0);
77 //mShareDenominator->SetDirectory(0);
78 //mRatio->SetDirectory(0);
79
80 // to enable error bar calculation...
81 fShareNumerator->Sumw2();
82 fShareDenominator->Sumw2();
83
84 fQualityNumerator->Sumw2();
85 fQualityDenominator->Sumw2();
86}
87
88//____________________________
89AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
90 delete fShareNumerator;
91 delete fShareDenominator;
92 delete fQualityNumerator;
93 delete fQualityDenominator;
94}
95//_________________________
96void AliFemtoShareQualityCorrFctn::Finish(){
97 // here is where we should normalize, fit, etc...
98 // we should NOT Draw() the histos (as I had done it below),
99 // since we want to insulate ourselves from root at this level
100 // of the code. Do it instead at root command line with browser.
101 // mShareNumerator->Draw();
102 //mShareDenominator->Draw();
103 //mRatio->Draw();
104
105}
106
107//____________________________
108AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
109 string stemp = "Qinv Correlation Function Report:\n";
110 char ctemp[100];
111 sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
112 stemp += ctemp;
113 sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
114 stemp += ctemp;
115 // stemp += mCoulombWeight->Report();
116 AliFemtoString returnThis = stemp;
117 return returnThis;
118}
119//____________________________
120void AliFemtoShareQualityCorrFctn::AddRealPair(const AliFemtoPair* pair){
121 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
122 Int_t nh = 0;
123 Int_t an = 0;
124 Int_t ns = 0;
125
126 for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
127 // If both have clusters in the same row
128 if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
129 pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
130 // Do they share it ?
131 if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
132 pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
133 {
134 if (Qinv < 0.01) {
135 cout << "Shared cluster in row " << imap << endl;
136 }
137 an++;
138 nh+=2;
139 ns+=2;
140 }
141
142 // Different hits on the same padrow
143 else {
144 an--;
145 nh+=2;
146 }
147 }
148 else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
149 pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
150 // One track has a hit, the other does not
151 an++;
152 nh++;
153 }
154 }
155 if (Qinv < 0.01) {
156 cout << "Qinv of the pair is " << Qinv << endl;
157 cout << "Clusters: " << endl;
158 for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
159 cout << imap ;
160 if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
161 else cout << " 0 " ;
162 if (pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
163 else cout << " 0 " ;
164 cout << " ";
165 if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
166 else cout << " X ";
167 if (pair->track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
168 else cout << " X ";
169 cout << endl;
170 }
171 }
172
173 Float_t hsmval = 0.0;
174 Float_t hsfval = 0.0;
175
176 if (nh >0) {
177 hsmval = an*1.0/nh;
178 hsfval = ns*1.0/nh;
179 }
180
181 if (Qinv < 0.01) {
182 cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->track1()->Track() << " " << pair->track2()->Track() << endl;
183 }
184
185 fShareNumerator->Fill(Qinv, hsfval);
186 fQualityNumerator->Fill(Qinv, hsmval);
187 // cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
188 //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
189}
190//____________________________
191void AliFemtoShareQualityCorrFctn::AddMixedPair(const AliFemtoPair* pair){
192 double weight = 1.0;
193 double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs...
194 Int_t nh = 0;
195 Int_t an = 0;
196 Int_t ns = 0;
197
198 for (unsigned int imap=0; imap<pair->track1()->Track()->TPCclusters().GetNbits(); imap++) {
199 // If both have clusters in the same row
200 if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) &&
201 pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
202 // Do they share it ?
203 if (pair->track1()->Track()->TPCsharing().TestBitNumber(imap) ||
204 pair->track2()->Track()->TPCsharing().TestBitNumber(imap))
205 {
206 // cout << "A shared cluster !!!" << endl;
207 // cout << "imap idx1 idx2 "
208 // << imap << " "
209 // << tP1idx[imap] << " " << tP2idx[imap] << endl;
210 an++;
211 nh+=2;
212 ns+=2;
213 }
214
215 // Different hits on the same padrow
216 else {
217 an--;
218 nh+=2;
219 }
220 }
221 else if (pair->track1()->Track()->TPCclusters().TestBitNumber(imap) ||
222 pair->track2()->Track()->TPCclusters().TestBitNumber(imap)) {
223 // One track has a hit, the other does not
224 an++;
225 nh++;
226 }
227 }
228
229 Float_t hsmval = 0.0;
230 Float_t hsfval = 0.0;
231
232 if (nh >0) {
233 hsmval = an*1.0/nh;
234 hsfval = ns*1.0/nh;
235 }
236
237 fShareDenominator->Fill(Qinv,hsfval,weight);
238 fQualityDenominator->Fill(Qinv,hsmval,weight);
239}
240
241
242void AliFemtoShareQualityCorrFctn::WriteHistos()
243{
244 fShareNumerator->Write();
245 fShareDenominator->Write();
246 fQualityNumerator->Write();
247 fQualityDenominator->Write();
248
249}