Making the directory structure of AliFemto flat. All files go into one common directory
[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 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 //____________________________
60 AliFemtoShareQualityCorrFctn::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 //____________________________
76 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
77   delete fShareNumerator;
78   delete fShareDenominator;
79   delete fQualityNumerator;
80   delete fQualityDenominator;
81 }
82 //_________________________
83 AliFemtoShareQualityCorrFctn& 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 //_________________________
108 void 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 //____________________________
120 AliFemtoString 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 //____________________________
132 void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
133   double Qinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
134   Int_t nh = 0;
135   Int_t an = 0;
136   Int_t ns = 0;
137   
138   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
139     // If both have clusters in the same row
140     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
141         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
142       // Do they share it ?
143       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
144           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
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     }
160     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
161              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
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;
170     for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
171       cout << imap ;
172       if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
173       else cout << " 0 " ;
174       if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
175       else cout << " 0 " ;
176       cout << "     ";
177       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
178       else cout << " X ";
179       if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
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) {
194     cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
195   }
196
197   fShareNumerator->Fill(Qinv, hsfval);
198   fQualityNumerator->Fill(Qinv, hsmval);
199   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << Qinv <<
200   //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
201 }
202 //____________________________
203 void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
204   double weight = 1.0;
205   double Qinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
206   Int_t nh = 0;
207   Int_t an = 0;
208   Int_t ns = 0;
209   
210   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
211     // If both have clusters in the same row
212     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
213         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
214       // Do they share it ?
215       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) ||
216           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
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     }
233     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
234              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
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
254 void AliFemtoShareQualityCorrFctn::WriteHistos()
255 {
256   fShareNumerator->Write();
257   fShareDenominator->Write();
258   fQualityNumerator->Write();
259   fQualityDenominator->Write();
260   
261 }