Making the directory structure of AliFemtoUser flat. All files go into one common...
[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)";
26 char TitNum[100] = "NumShare";
27 strcat(TitNum,title);
28 fShareNumerator = new TH2D(TitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
29 // set up denominator
30 //title = "Den Qinv (MeV/c)";
31 char TitDen[100] = "DenShare";
32 strcat(TitDen,title);
33 fShareDenominator = new TH2D(TitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
34
35 char Tit2Num[100] = "NumQuality";
36 strcat(Tit2Num,title);
37 fQualityNumerator = new TH2D(Tit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
38 // set up denominator
39 //title = "Den Qinv (MeV/c)";
40 char Tit2Den[100] = "DenQuality";
41 strcat(Tit2Den,title);
42 fQualityDenominator = new TH2D(Tit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
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{
66 if (aCorrFctn.fShareNumerator)
67 fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
68 if (aCorrFctn.fShareDenominator)
69 fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
70 if (aCorrFctn.fQualityNumerator)
71 fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
72 if (aCorrFctn.fQualityDenominator)
73 fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
74}
75//____________________________
67427ff7 76AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
77 delete fShareNumerator;
78 delete fShareDenominator;
79 delete fQualityNumerator;
80 delete fQualityDenominator;
81}
82//_________________________
0215f606 83AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
84{
85 if (this == &aCorrFctn)
86 return *this;
87
88 if (aCorrFctn.fShareNumerator)
89 fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
90 else
91 fShareNumerator = 0;
92 if (aCorrFctn.fShareDenominator)
93 fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
94 else
95 fShareDenominator = 0;
96 if (aCorrFctn.fQualityNumerator)
97 fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
98 else
99 fQualityNumerator = 0;
100 if (aCorrFctn.fQualityDenominator)
101 fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
102 else
103 fQualityDenominator = 0;
104
105 return *this;
106}
107//_________________________
67427ff7 108void AliFemtoShareQualityCorrFctn::Finish(){
109 // here is where we should normalize, fit, etc...
110 // we should NOT Draw() the histos (as I had done it below),
111 // since we want to insulate ourselves from root at this level
112 // of the code. Do it instead at root command line with browser.
113 // mShareNumerator->Draw();
114 //mShareDenominator->Draw();
115 //mRatio->Draw();
116
117}
118
119//____________________________
120AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
121 string stemp = "Qinv Correlation Function Report:\n";
122 char ctemp[100];
123 sprintf(ctemp,"Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
124 stemp += ctemp;
125 sprintf(ctemp,"Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
126 stemp += ctemp;
127 // stemp += mCoulombWeight->Report();
128 AliFemtoString returnThis = stemp;
129 return returnThis;
130}
131//____________________________
65423af9 132void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
133 double Qinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
67427ff7 134 Int_t nh = 0;
135 Int_t an = 0;
136 Int_t ns = 0;
137
65423af9 138 for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
67427ff7 139 // If both have clusters in the same row
65423af9 140 if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
141 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 142 // Do they share it ?
65423af9 143 if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
144 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
67427ff7 145 {
146 if (Qinv < 0.01) {
147 cout << "Shared cluster in row " << imap << endl;
148 }
149 an++;
150 nh+=2;
151 ns+=2;
152 }
153
154 // Different hits on the same padrow
155 else {
156 an--;
157 nh+=2;
158 }
159 }
65423af9 160 else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
161 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 162 // One track has a hit, the other does not
163 an++;
164 nh++;
165 }
166 }
167 if (Qinv < 0.01) {
168 cout << "Qinv of the pair is " << Qinv << endl;
169 cout << "Clusters: " << endl;
65423af9 170 for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
67427ff7 171 cout << imap ;
65423af9 172 if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
67427ff7 173 else cout << " 0 " ;
65423af9 174 if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
67427ff7 175 else cout << " 0 " ;
176 cout << " ";
65423af9 177 if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
67427ff7 178 else cout << " X ";
65423af9 179 if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
67427ff7 180 else cout << " X ";
181 cout << endl;
182 }
183 }
184
185 Float_t hsmval = 0.0;
186 Float_t hsfval = 0.0;
187
188 if (nh >0) {
189 hsmval = an*1.0/nh;
190 hsfval = ns*1.0/nh;
191 }
192
193 if (Qinv < 0.01) {
65423af9 194 cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
67427ff7 195 }
196
197 fShareNumerator->Fill(Qinv, hsfval);
198 fQualityNumerator->Fill(Qinv, hsmval);
199 // cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
65423af9 200 //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
67427ff7 201}
202//____________________________
65423af9 203void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
67427ff7 204 double weight = 1.0;
65423af9 205 double Qinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
67427ff7 206 Int_t nh = 0;
207 Int_t an = 0;
208 Int_t ns = 0;
209
65423af9 210 for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
67427ff7 211 // If both have clusters in the same row
65423af9 212 if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) &&
213 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 214 // Do they share it ?
65423af9 215 if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) ||
216 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
67427ff7 217 {
218 // cout << "A shared cluster !!!" << endl;
219 // cout << "imap idx1 idx2 "
220 // << imap << " "
221 // << tP1idx[imap] << " " << tP2idx[imap] << endl;
222 an++;
223 nh+=2;
224 ns+=2;
225 }
226
227 // Different hits on the same padrow
228 else {
229 an--;
230 nh+=2;
231 }
232 }
65423af9 233 else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
234 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
67427ff7 235 // One track has a hit, the other does not
236 an++;
237 nh++;
238 }
239 }
240
241 Float_t hsmval = 0.0;
242 Float_t hsfval = 0.0;
243
244 if (nh >0) {
245 hsmval = an*1.0/nh;
246 hsfval = ns*1.0/nh;
247 }
248
249 fShareDenominator->Fill(Qinv,hsfval,weight);
250 fQualityDenominator->Fill(Qinv,hsmval,weight);
251}
252
253
254void AliFemtoShareQualityCorrFctn::WriteHistos()
255{
256 fShareNumerator->Write();
257 fShareDenominator->Write();
258 fQualityNumerator->Write();
259 fQualityDenominator->Write();
260
261}