Port of changes from v4-07-Release and additional rule conformance
[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,100,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,100,0.0,1.00001);
34
35   char tTit2Num[100] = "NumQuality";
36   strcat(tTit2Num,title);
37   fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,150,-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,150,-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) && (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 //   }
238
239 //   fShareNumerator->Fill(tQinv, hsfval);
240 //   fQualityNumerator->Fill(tQinv, hsmval);
241   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
242   //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
243 }
244 //____________________________
245 void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
246   // add mixed (background) pair
247   double weight = 1.0;
248   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
249   Int_t nh = 0;
250   Int_t an = 0;
251   Int_t ns = 0;
252   
253   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
254     // If both have clusters in the same row
255     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
256         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
257       // Do they share it ?
258       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
259           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
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     }
276     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
277              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
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
292   fShareDenominator->Fill(tQinv,hsfval,weight);
293   fQualityDenominator->Fill(tQinv,hsmval,weight);
294 }
295
296
297 void AliFemtoShareQualityCorrFctn::WriteHistos()
298 {
299   fShareNumerator->Write();
300   fShareDenominator->Write();
301   fQualityNumerator->Write();
302   fQualityDenominator->Write();
303   
304 }