c9dc4f730778300a8c18e74eb8ad9149fe8438a7
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoShareQualityCorrFctn.cxx
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 ////////////////////////////////////////////////////////////////////////////////
8
9 #include "AliFemtoShareQualityCorrFctn.h"
10 //#include "AliFemtoHisto.hh"
11 #include <cstdio>
12
13 #ifdef __ROOT__ 
14 ClassImp(AliFemtoShareQualityCorrFctn)
15 #endif
16
17 //____________________________
18 AliFemtoShareQualityCorrFctn::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 {
24   // set up numerator
25   //  title = "Num Qinv (MeV/c)";
26   char tTitNum[100] = "NumShare";
27   strcat(tTitNum,title);
28   fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
29   // set up denominator
30   //title = "Den Qinv (MeV/c)";
31   char tTitDen[100] = "DenShare";
32   strcat(tTitDen,title);
33   fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001);
34
35   char tTit2Num[100] = "NumQuality";
36   strcat(tTit2Num,title);
37   fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001);
38   // set up denominator
39   //title = "Den Qinv (MeV/c)";
40   char tTit2Den[100] = "DenQuality";
41   strcat(tTit2Den,title);
42   fQualityDenominator = new TH2D(tTit2Den,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 //____________________________
60 AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
61   fShareNumerator(0),
62   fShareDenominator(0),
63   fQualityNumerator(0),
64   fQualityDenominator(0)
65 {
66   // copy constructor
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 //____________________________
77 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
78   // destructor
79   delete fShareNumerator;
80   delete fShareDenominator;
81   delete fQualityNumerator;
82   delete fQualityDenominator;
83 }
84 //_________________________
85 AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
86 {
87   // assignment operator
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 //_________________________
111 void 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 //____________________________
123 AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
124   // create report
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 //____________________________
136 void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
137   // add real (effect) pair
138   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
139   Int_t nh = 0;
140   Int_t an = 0;
141   Int_t ns = 0;
142   
143   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
144     // If both have clusters in the same row
145     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
146         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
147       // Do they share it ?
148       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
149           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
150         {
151           if (tQinv < 0.01) {
152             cout << "Shared cluster in row " << imap << endl; 
153           }
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     }
165     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
166              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
167       // One track has a hit, the other does not
168       an++;
169       nh++;
170     }
171   }
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   }
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
198   if (tQinv < 0.01) {
199     cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
200   }
201
202   fShareNumerator->Fill(tQinv, hsfval);
203   fQualityNumerator->Fill(tQinv, hsmval);
204   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
205   //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
206 }
207 //____________________________
208 void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
209   // add mixed (background) pair
210   double weight = 1.0;
211   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
212   Int_t nh = 0;
213   Int_t an = 0;
214   Int_t ns = 0;
215   
216   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
217     // If both have clusters in the same row
218     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
219         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
220       // Do they share it ?
221       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) ||
222           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
223         {
224           //      cout << "A shared cluster !!!" << endl;
225           //    cout << "imap idx1 idx2 " 
226           //         << imap << " "
227           //         << tP1idx[imap] << " " << tP2idx[imap] << endl;
228           an++;
229           nh+=2;
230           ns+=2;
231         }
232       
233       // Different hits on the same padrow
234       else {
235         an--;
236         nh+=2;
237       }
238     }
239     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
240              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
241       // One track has a hit, the other does not
242       an++;
243       nh++;
244     }
245   }
246   
247   Float_t hsmval = 0.0;
248   Float_t hsfval = 0.0;
249
250   if (nh >0) {
251     hsmval = an*1.0/nh;
252     hsfval = ns*1.0/nh;
253   }
254
255   fShareDenominator->Fill(tQinv,hsfval,weight);
256   fQualityDenominator->Fill(tQinv,hsmval,weight);
257 }
258
259
260 void AliFemtoShareQualityCorrFctn::WriteHistos()
261 {
262   fShareNumerator->Write();
263   fShareDenominator->Write();
264   fQualityNumerator->Write();
265   fQualityDenominator->Write();
266   
267 }