Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[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   AliFemtoCorrFctn(),
20   fShareNumerator(0),
21   fShareDenominator(0),
22   fQualityNumerator(0),
23   fQualityDenominator(0),
24   fTPCSepNumerator(0),
25   fTPCSepDenominator(0)
26 {
27   // set up numerator
28   //  title = "Num Qinv (MeV/c)";
29   char tTitNum[101] = "NumShare";
30   strncat(tTitNum,title, 100);
31   fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
32   // set up denominator
33   //title = "Den Qinv (MeV/c)";
34   char tTitDen[101] = "DenShare";
35   strncat(tTitDen,title, 100);
36   fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,1.00001);
37
38   char tTit2Num[101] = "NumQuality";
39   strncat(tTit2Num,title, 100);
40   fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
41   // set up denominator
42   //title = "Den Qinv (MeV/c)";
43   char tTit2Den[101] = "DenQuality";
44   strncat(tTit2Den,title, 100);
45   fQualityDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001);
46   // set up ratio
47   //title = "Ratio Qinv (MeV/c)";
48   // this next bit is unfortunately needed so that we can have many histos of same "title"
49   // it is neccessary if we typedef TH2D to TH1d (which we do)
50   //mShareNumerator->SetDirectory(0);
51   //mShareDenominator->SetDirectory(0);
52   //mRatio->SetDirectory(0);
53
54   char tTit3Num[101] = "NumTPCSep";
55   strncat(tTit3Num,title, 100);
56   fTPCSepNumerator = new TH2D(tTit3Num,title,nbins,QinvLo,QinvHi,150,0.0,100.0);
57   // set up denominator
58   //title = "Den Qinv (MeV/c)";
59   char tTit3Den[101] = "DenTPCSep";
60   strncat(tTit3Den,title, 100);
61   fTPCSepDenominator = new TH2D(tTit3Den,title,nbins,QinvLo,QinvHi,150,0.0,100.0);
62
63   // to enable error bar calculation...
64   fShareNumerator->Sumw2();
65   fShareDenominator->Sumw2();
66
67   fQualityNumerator->Sumw2();
68   fQualityDenominator->Sumw2();
69   
70   fTPCSepNumerator->Sumw2();
71   fTPCSepDenominator->Sumw2();
72 }
73
74 //____________________________
75 AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(const AliFemtoShareQualityCorrFctn& aCorrFctn) :
76   AliFemtoCorrFctn(),
77   fShareNumerator(0),
78   fShareDenominator(0),
79   fQualityNumerator(0),
80   fQualityDenominator(0),
81   fTPCSepNumerator(0),
82   fTPCSepDenominator(0)
83 {
84   // copy constructor
85   if (aCorrFctn.fShareNumerator)
86     fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
87   if (aCorrFctn.fShareDenominator)
88     fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
89   if (aCorrFctn.fQualityNumerator)
90     fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
91   if (aCorrFctn.fQualityDenominator)
92     fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
93   if (aCorrFctn.fTPCSepNumerator)
94     fTPCSepNumerator = new TH2D(*aCorrFctn.fTPCSepNumerator);
95   if (aCorrFctn.fTPCSepDenominator)
96     fTPCSepDenominator = new TH2D(*aCorrFctn.fTPCSepDenominator);
97 }
98 //____________________________
99 AliFemtoShareQualityCorrFctn::~AliFemtoShareQualityCorrFctn(){
100   // destructor
101   delete fShareNumerator;
102   delete fShareDenominator;
103   delete fQualityNumerator;
104   delete fQualityDenominator;
105   delete fTPCSepNumerator;
106   delete fTPCSepDenominator;
107 }
108 //_________________________
109 AliFemtoShareQualityCorrFctn& AliFemtoShareQualityCorrFctn::operator=(const AliFemtoShareQualityCorrFctn& aCorrFctn)
110 {
111   // assignment operator
112   if (this == &aCorrFctn)
113     return *this;
114
115   if (aCorrFctn.fShareNumerator)
116     fShareNumerator = new TH2D(*aCorrFctn.fShareNumerator);
117   else
118     fShareNumerator = 0;
119   if (aCorrFctn.fShareDenominator)
120     fShareDenominator = new TH2D(*aCorrFctn.fShareDenominator);
121   else
122     fShareDenominator = 0;
123   if (aCorrFctn.fQualityNumerator)
124     fQualityNumerator = new TH2D(*aCorrFctn.fQualityNumerator);
125   else
126     fQualityNumerator = 0;
127   if (aCorrFctn.fQualityDenominator)
128     fQualityDenominator = new TH2D(*aCorrFctn.fQualityDenominator);
129   else
130     fQualityDenominator = 0;
131   if (aCorrFctn.fTPCSepNumerator)
132     fTPCSepNumerator = new TH2D(*aCorrFctn.fTPCSepNumerator);
133   else
134     fTPCSepNumerator = 0;
135   if (aCorrFctn.fTPCSepDenominator)
136     fTPCSepDenominator = new TH2D(*aCorrFctn.fTPCSepDenominator);
137   else
138     fTPCSepDenominator = 0;
139
140   return *this;
141 }
142 //_________________________
143 void AliFemtoShareQualityCorrFctn::Finish(){
144   // here is where we should normalize, fit, etc...
145   // we should NOT Draw() the histos (as I had done it below),
146   // since we want to insulate ourselves from root at this level
147   // of the code.  Do it instead at root command line with browser.
148   //  mShareNumerator->Draw();
149   //mShareDenominator->Draw();
150   //mRatio->Draw();
151
152 }
153
154 //____________________________
155 AliFemtoString AliFemtoShareQualityCorrFctn::Report(){
156   // create report
157   string stemp = "Qinv Correlation Function Report:\n";
158   char ctemp[100];
159   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fShareNumerator->GetEntries());
160   stemp += ctemp;
161   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fShareDenominator->GetEntries());
162   stemp += ctemp;
163   //  stemp += mCoulombWeight->Report();
164   AliFemtoString returnThis = stemp;
165   return returnThis;
166 }
167 //____________________________
168 void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){
169   // add real (effect) pair
170   if (fPairCut)
171     if (!fPairCut->Pass(pair)) return;
172
173   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
174   Int_t nh = 0;
175   Int_t an = 0;
176   Int_t ns = 0;
177   
178   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
179     // If both have clusters in the same row
180     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
181         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
182       // Do they share it ?
183       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
184           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
185         {
186 //        if (tQinv < 0.01) {
187 //          cout << "Shared cluster in row " << imap << endl; 
188 //        }
189           an++;
190           nh+=2;
191           ns+=2;
192         }
193       
194       // Different hits on the same padrow
195       else {
196         an--;
197         nh+=2;
198       }
199     }
200     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
201              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
202       // One track has a hit, the other does not
203       an++;
204       nh++;
205     }
206   }
207 //    if (tQinv < 0.01) {
208 //     cout << "Qinv of the pair is " << tQinv << endl;
209 //     cout << "Clusters: " << endl;
210 //     for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
211 //       cout << imap ;
212 //       if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
213 //       else cout << " 0 " ;
214 //       if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
215 //       else cout << " 0 " ;
216 //       cout << "     ";
217 //       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
218 //       else cout << " X ";
219 //       if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
220 //       else cout << " X ";
221 //       cout << endl;
222 //     }
223 //   }
224
225   Float_t hsmval = 0.0;
226   Float_t hsfval = 0.0;
227
228   if (nh >0) {
229     hsmval = an*1.0/nh;
230     hsfval = ns*1.0/nh;
231   }
232
233 //   if ((tQinv < 0.005) && (hsmval<-0.0)) {
234 //     cout << "Quality  Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl;
235 //     cout << "Qinv of the pair is " << tQinv << endl;
236 //     cout << "Clusters: " << endl;
237 //     for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
238 //       cout << imap ;
239 //       if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
240 //       else cout << " 0 " ;
241 //       if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 ";
242 //       else cout << " 0 " ;
243 //       cout << "     ";
244 //       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
245 //       else cout << " X ";
246 //       if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S ";
247 //       else cout << " X ";
248 //       cout << endl;
249 //     }
250 //     cout << "Momentum1 " 
251 //       << pair->Track1()->Track()->P().x() << " " 
252 //       << pair->Track1()->Track()->P().y() << " "  
253 //       << pair->Track1()->Track()->P().z() << " "  
254 //       << pair->Track1()->Track()->Label() << " "  
255 //       << pair->Track1()->Track()->TrackId() << " "  
256 //       << pair->Track1()->Track()->Flags() << " "
257 //       << pair->Track1()->Track()->KinkIndex(0) << " "
258 //       << pair->Track1()->Track()->KinkIndex(1) << " "
259 //       << pair->Track1()->Track()->KinkIndex(2) << " "
260 //       << pair->Track1()->Track()->ITSchi2() << " "
261 //       << pair->Track1()->Track()->ITSncls() << " "
262 //       << pair->Track1()->Track()->TPCchi2() << " "
263 //       << pair->Track1()->Track()->TPCncls() << " "
264 //       << endl;
265 //     cout << "Momentum2 " 
266 //       << pair->Track2()->Track()->P().x() << " "  
267 //       << pair->Track2()->Track()->P().y() << " "  
268 //       << pair->Track2()->Track()->P().z() << " "  
269 //       << pair->Track2()->Track()->Label() << " "  
270 //       << pair->Track2()->Track()->TrackId() << " "  
271 //       << pair->Track2()->Track()->Flags() << " " 
272 //       << pair->Track2()->Track()->KinkIndex(0) << " "
273 //       << pair->Track2()->Track()->KinkIndex(1) << " "
274 //       << pair->Track2()->Track()->KinkIndex(2) << " "
275 //       << pair->Track2()->Track()->ITSchi2() << " "
276 //       << pair->Track2()->Track()->ITSncls() << " "
277 //       << pair->Track2()->Track()->TPCchi2() << " "
278 //       << pair->Track2()->Track()->TPCncls() << " "
279 //       << endl;
280 //   }
281
282   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
283   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
284   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
285   double dist = sqrt(distx*distx + disty*disty + distz*distz);
286
287   fShareNumerator->Fill(tQinv, hsfval);
288   fQualityNumerator->Fill(tQinv, hsmval);
289   fTPCSepNumerator->Fill(tQinv, dist);
290
291   //  cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
292   //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl;
293 }
294 //____________________________
295 void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){
296   // add mixed (background) pair
297   if (fPairCut)
298     if (!fPairCut->Pass(pair)) return;
299
300   double weight = 1.0;
301   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
302   Int_t nh = 0;
303   Int_t an = 0;
304   Int_t ns = 0;
305   
306   for (unsigned int imap=0; imap<pair->Track1()->Track()->TPCclusters().GetNbits(); imap++) {
307     // If both have clusters in the same row
308     if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && 
309         pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
310       // Do they share it ?
311       if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) &&
312           pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
313         {
314           //      cout << "A shared cluster !!!" << endl;
315           //    cout << "imap idx1 idx2 " 
316           //         << imap << " "
317           //         << tP1idx[imap] << " " << tP2idx[imap] << endl;
318           an++;
319           nh+=2;
320           ns+=2;
321         }
322       
323       // Different hits on the same padrow
324       else {
325         an--;
326         nh+=2;
327       }
328     }
329     else if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) ||
330              pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
331       // One track has a hit, the other does not
332       an++;
333       nh++;
334     }
335   }
336   
337   Float_t hsmval = 0.0;
338   Float_t hsfval = 0.0;
339
340   if (nh >0) {
341     hsmval = an*1.0/nh;
342     hsfval = ns*1.0/nh;
343   }
344
345   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
346   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
347   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
348   double dist = sqrt(distx*distx + disty*disty + distz*distz);
349
350   fShareDenominator->Fill(tQinv,hsfval,weight);
351   fQualityDenominator->Fill(tQinv,hsmval,weight);
352   fTPCSepDenominator->Fill(tQinv, dist, weight);
353
354 }
355
356
357 void AliFemtoShareQualityCorrFctn::WriteHistos()
358 {
359   // Write out result histograms
360   fShareNumerator->Write();
361   fShareDenominator->Write();
362   fQualityNumerator->Write();
363   fQualityDenominator->Write();
364   fTPCSepNumerator->Write();
365   fTPCSepDenominator->Write();
366   
367 }
368 //______________________________
369 TList* AliFemtoShareQualityCorrFctn::GetOutputList()
370 {
371   // Prepare the list of objects to be written to the output
372   TList *tOutputList = new TList();
373
374   tOutputList->Add(fShareNumerator); 
375   tOutputList->Add(fShareDenominator);  
376   tOutputList->Add(fQualityNumerator);  
377   tOutputList->Add(fQualityDenominator);  
378   tOutputList->Add(fTPCSepNumerator);  
379   tOutputList->Add(fTPCSepDenominator);  
380
381   return tOutputList;
382 }