Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoTPCInnerCorrFctn.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 ///                                                                          ///
3 /// AliFemtoTPCInnerCorrFctn - A correlation function that saves the         ///
4 /// distance at the entrance to the TPC between two tracks as a function     ///
5 /// of qinv                                                                  ///
6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
7 ///                                                                          ///
8 ////////////////////////////////////////////////////////////////////////////////
9
10 #include "AliFemtoTPCInnerCorrFctn.h"
11 //#include "AliFemtoHisto.hh"
12 #include <cstdio>
13
14 #ifdef __ROOT__ 
15 ClassImp(AliFemtoTPCInnerCorrFctn)
16 #endif
17
18 //____________________________
19 AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
20   fDTPCNumerator(0),
21   fDTPCDenominator(0),
22   fRadDNumerator(0),
23   fRadDDenominator(0),
24   fRadius(100)
25 {
26   // set up numerator
27   //  title = "Num Qinv (MeV/c)";
28   char tTitNum[101] = "NumDTPC";
29   strncat(tTitNum,title, 100);
30   fDTPCNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
31   // set up denominator
32   //title = "Den Qinv (MeV/c)";
33   char tTitDen[101] = "DenDTPC";
34   strncat(tTitDen,title, 100);
35   fDTPCDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,20.0);
36
37   char tTitNumR[101] = "NumRadD";
38   strncat(tTitNumR,title, 100);
39   fRadDNumerator = new TH2D(tTitNumR,title,50,-0.1,0.1,50,-0.1,0.1);
40   // set up denominator
41   //title = "Den Qinv (MeV/c)";
42   char tTitDenR[101] = "DenRadD";
43   strncat(tTitDenR,title, 100);
44   fRadDDenominator = new TH2D(tTitDenR,title,50,-0.1,0.1,50,-0.1,0.1);
45
46   // to enable error bar calculation...
47   fDTPCNumerator->Sumw2();
48   fDTPCDenominator->Sumw2();
49   fRadDNumerator->Sumw2();
50   fRadDDenominator->Sumw2();
51 }
52
53 //____________________________
54 AliFemtoTPCInnerCorrFctn::AliFemtoTPCInnerCorrFctn(const AliFemtoTPCInnerCorrFctn& aCorrFctn) :
55   AliFemtoCorrFctn(),
56   fDTPCNumerator(0),
57   fDTPCDenominator(0),
58   fRadDNumerator(0),
59   fRadDDenominator(0),
60   fRadius(100)
61 {
62   // copy constructor
63   if (aCorrFctn.fDTPCNumerator)
64     fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
65   if (aCorrFctn.fDTPCDenominator)
66     fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
67   if (aCorrFctn.fRadDNumerator)
68     fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
69   if (aCorrFctn.fRadDDenominator)
70     fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
71 }
72 //____________________________
73 AliFemtoTPCInnerCorrFctn::~AliFemtoTPCInnerCorrFctn(){
74   // destructor
75   delete fDTPCNumerator;
76   delete fDTPCDenominator;
77   delete fRadDNumerator;
78   delete fRadDDenominator;
79 }
80 //_________________________
81 AliFemtoTPCInnerCorrFctn& AliFemtoTPCInnerCorrFctn::operator=(const AliFemtoTPCInnerCorrFctn& aCorrFctn)
82 {
83   // assignment operator
84   if (this == &aCorrFctn)
85     return *this;
86
87   if (aCorrFctn.fDTPCNumerator)
88     fDTPCNumerator = new TH2D(*aCorrFctn.fDTPCNumerator);
89   else
90     fDTPCNumerator = 0;
91   if (aCorrFctn.fDTPCDenominator)
92     fDTPCDenominator = new TH2D(*aCorrFctn.fDTPCDenominator);
93   else
94     fDTPCDenominator = 0;
95
96   if (aCorrFctn.fRadDNumerator)
97     fRadDNumerator = new TH2D(*aCorrFctn.fRadDNumerator);
98   else
99     fRadDNumerator = 0;
100   if (aCorrFctn.fRadDDenominator)
101     fRadDDenominator = new TH2D(*aCorrFctn.fRadDDenominator);
102   else
103     fRadDDenominator = 0;
104
105   fRadius = aCorrFctn.fRadius;
106
107   return *this;
108 }
109 //_________________________
110 void AliFemtoTPCInnerCorrFctn::Finish(){
111   // here is where we should normalize, fit, etc...
112   // we should NOT Draw() the histos (as I had done it below),
113   // since we want to insulate ourselves from root at this level
114   // of the code.  Do it instead at root command line with browser.
115   //  mShareNumerator->Draw();
116   //mShareDenominator->Draw();
117   //mRatio->Draw();
118
119 }
120
121 //____________________________
122 AliFemtoString AliFemtoTPCInnerCorrFctn::Report(){
123   // create report
124   string stemp = "Entrace TPC distance Correlation Function Report:\n";
125   char ctemp[100];
126   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDTPCNumerator->GetEntries());
127   stemp += ctemp;
128   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDTPCDenominator->GetEntries());
129   stemp += ctemp;
130   //  stemp += mCoulombWeight->Report();
131   AliFemtoString returnThis = stemp;
132   return returnThis;
133 }
134 //____________________________
135 void AliFemtoTPCInnerCorrFctn::AddRealPair( AliFemtoPair* pair){
136   // add real (effect) pair
137   if (fPairCut)
138     if (!fPairCut->Pass(pair)) return;
139
140   double pih = TMath::Pi();
141   double pit = TMath::Pi()*2;
142
143   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
144   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
145   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
146   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
147   double dist = sqrt(distx*distx + disty*disty + distz*distz);
148
149   fDTPCNumerator->Fill(tQinv, dist);
150
151   if (tQinv < 0.1) {
152     double phi1 = pair->Track1()->Track()->P().Phi();
153     double phi2 = pair->Track2()->Track()->P().Phi();
154     double chg1 = pair->Track1()->Track()->Charge();
155     double chg2 = pair->Track2()->Track()->Charge();
156     double ptv1 = pair->Track1()->Track()->Pt();
157     double ptv2 = pair->Track2()->Track()->Pt();
158     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
159     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
160     double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
161     double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
162     double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
163
164     while (phid>pih) phid -= pit;
165     while (phid<-pih) phid += pit;
166     //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
167     double etad = eta2 - eta1;
168     fRadDNumerator->Fill(phid, etad);
169   }
170 }
171 //____________________________
172 void AliFemtoTPCInnerCorrFctn::AddMixedPair( AliFemtoPair* pair){
173   // add mixed (background) pair
174   if (fPairCut)
175     if (!fPairCut->Pass(pair)) return;
176
177   double pih = TMath::Pi();
178   double pit = TMath::Pi()*2;
179
180   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
181   double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
182   double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
183   double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
184   double dist = sqrt(distx*distx + disty*disty + distz*distz);
185
186   fDTPCDenominator->Fill(tQinv,dist);
187
188   if (tQinv < 0.1) {
189     double phi1 = pair->Track1()->Track()->P().Phi();
190     double phi2 = pair->Track2()->Track()->P().Phi();
191     double chg1 = pair->Track1()->Track()->Charge();
192     double chg2 = pair->Track2()->Track()->Charge();
193     double ptv1 = pair->Track1()->Track()->Pt();
194     double ptv2 = pair->Track2()->Track()->Pt();
195     double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
196     double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
197     double arg1 = -0.3 * 0.5 * chg1 * fRadius/(2*ptv1);
198     double arg2 = -0.3 * 0.5 * chg2 * fRadius/(2*ptv2);
199     double phid = phi2 - phi1 + TMath::ASin(arg2) - TMath::ASin(arg1);
200
201     while (phid>pih) phid -= pit;
202     while (phid<-pih) phid += pit;
203     //    dist = phi2 - phi1 + TMath::ASin(-0.3 * 0.5 * chg2 * fRadius/(2*ptv2)) - TMath::ASin(-0.3 * 0.5 * chg1 * fRadius/(2*ptv1));
204     double etad = eta2 - eta1;
205     fRadDDenominator->Fill(phid, etad);
206   }
207 }
208
209
210 void AliFemtoTPCInnerCorrFctn::WriteHistos()
211 {
212   // Write out result histograms
213   fDTPCNumerator->Write();
214   fDTPCDenominator->Write();
215   fRadDNumerator->Write();
216   fRadDDenominator->Write();
217 }
218 //______________________________
219 TList* AliFemtoTPCInnerCorrFctn::GetOutputList()
220 {
221   // Prepare the list of objects to be written to the output
222   TList *tOutputList = new TList();
223
224   tOutputList->Add(fDTPCNumerator); 
225   tOutputList->Add(fDTPCDenominator);  
226   tOutputList->Add(fRadDNumerator); 
227   tOutputList->Add(fRadDDenominator);  
228
229   return tOutputList;
230 }
231
232 void AliFemtoTPCInnerCorrFctn::SetRadius(double rad)
233 {
234   fRadius = rad;
235 }