This commit was generated by cvs2svn to compensate for changes in r18145,
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / CorrFctn / AliFemtoShareQualityCorrFctn.cxx
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__ 
45 ClassImp(AliFemtoShareQualityCorrFctn)
46 #endif
47
48 //____________________________
49 AliFemtoShareQualityCorrFctn::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 //____________________________
86 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
87   delete fShareNumerator;
88   delete fShareDenominator;
89   delete fQualityNumerator;
90   delete fQualityDenominator;
91 }
92 //_________________________
93 void 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 //____________________________
105 AliFemtoString 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 //____________________________
117 void 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 //____________________________
188 void 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
239 void AliFemtoShareQualityCorrFctn::WriteHistos()
240 {
241   fShareNumerator->Write();
242   fShareDenominator->Write();
243   fQualityNumerator->Write();
244   fQualityDenominator->Write();
245   
246 }