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