doxy: correctly consider macro blocks in methods
[u/mrichter/AliRoot.git] / VZERO / V0_efficiencies.C
CommitLineData
4da3d4d3 1void V0_efficiencies()
2
3{
4// Caution: time windows and cuts on signals are HARD-coded
5// they need to be adjusted to the current configuration
6
7 gROOT->Reset();
8 rl = AliRunLoader::Open("galice.root");
9 rl->LoadgAlice();
10 gAlice = rl->GetAliRun();
11 Int_t nevent = rl->GetNumberOfEvents();
12 cout<<" --- Number of events in file = "<< nevent <<" ---"<<endl;
13
14//___________________________________________________
15// Book HISTOGRAMS
16
17 TH2F *hMapV0L = new TH2F("hMapV0L","V0L",161,-0.8,0.8,161,-0.8,0.8);
18 TH2F *hMapV0R = new TH2F("hMapV0R","V0R",161,-0.8,0.8,161,-0.8,0.8);
19 TH1F *hpdgV0 = new TH1F("hpdgV0","pdgV0",1001,-500,500);
20 TH1F *hvertexZ = new TH1F("hvertexZ","vertex Z",972,-40.,40.);
21 TH1F *hTimC = new TH1F("hTimC","Time of Flight in V0C",500,0,49);
22 TH1F *hTimA = new TH1F("hTimA","Time of Flight in V0A",500,0,49);
23 TH1F *hCell = new TH1F("hCell","Cell Number",100,0,99);
24 TH2F *hCellADC = new TH2F("hCellADC","ADC vs Cell",100,0,99,1000,0,3999);
25
26 TH1F *hMul0 = new TH1F("hMul0 ","Multiplicity in V0",80,0,79);
27 TH1F *hMulC0 = new TH1F("hMulC0","Multiplicity in V0C",50,0,49);
28 TH1F *hMulA0 = new TH1F("hMulA0","Multiplicity in V0A",50,0,49);
29 TH1F *hMulAnd0 = new TH1F("hMulAnd0","Trigger and",50,0,49);
30
31//___________________________________________________
32
33 AliLoader *ld = rl->GetLoader("VZEROLoader");
34 ld->GetHitsDataLoader()->Load("READ");
35
36 rl->LoadHeader();
37
38 // to access the particle Stack
39 rl->LoadKinematics("READ");
40
41 AliVZERO *v0 = (AliVZERO*)gAlice->GetDetector("VZERO");
42
43 bool flagV0L = false;
44 bool flagV0R = false;
45 Float_t timeV0L = 1e12;
46 Float_t timeV0R = 1e12;
47
48 Int_t fNdigits = 0;
49 Int_t fMulA = 0;
50 Int_t fMulC = 0;
51 Int_t fDigits = 0;
52 Float_t fPhotoCathodeEfficiency = 0.18;
53
54
55 Int_t nVOL = 0;
56 Int_t nVOR = 0;
57 Int_t nVOLetR = 0;
58 Int_t nVOLouR = 0;
59
60 TDatabasePDG * pdgdb = TDatabasePDG::Instance();
61
62
63 Float_t fPMVoltage = 768.0;
64 Float_t fPMGain1 = TMath::Power((fPMVoltage / 112.5) ,7.04277);
65 Float_t fPMGain[64];
66 Float_t cPM[64];
67
68 for(Int_t ii=0; ii<64; ii++){
69 fPMGain[ii] = gRandom->Gaus(fPMGain1, fPMGain1/5); // 20% uncertainty on the PM gain
70 cPM[ii] = fPhotoCathodeEfficiency * fPMGain[ii];}
71
72 Float_t kC = 1e-11;
73 Float_t kthau = 2.1*1e-9;
74 Float_t ktheta = 50.0 * kC;
75 Float_t kQe = 1.6e-19;
76 Float_t coef = 200.0; // p-p
77// Float_t coef = 1.0; // Pb-Pb
78
79 Int_t map[80];
80 for(Int_t i=0; i<80; i++) map[i] = 0;
81 Int_t hit[80];
82 for(Int_t i=0; i<80; i++) hit[i] = 0;
83
84 for(Int_t event = 0; event < nevent; event++){
85 Float_t timeV0L = 1e12;
86 Float_t timeV0R = 1e12;
87
88 rl->GetEvent(event);
89 for(Int_t i=0; i<80; i++) map[i] = 0;
90 for(Int_t i=0; i<80; i++) hit[i] = 0;
91 ld->LoadHits();
92 v0->SetTreeAddress();
93 TTree *treeH = ld->TreeH();
94 Int_t ntracks = treeH->GetEntries();
95
96 for (Int_t itrack=0; itrack<ntracks;itrack++) {
97
98 v0->ResetHits();
99 treeH->GetEvent(itrack);
100 if(v0){
101 for (AliVZEROhit *vhit=(AliVZEROhit*)v0->FirstHit(itrack);
102 vhit; vhit=(AliVZEROhit*)v0->NextHit())
103 {
104// vhit->Dump();
105 hpdgV0->Fill(vhit->TrackPiD());
106 hvertexZ->Fill(vhit->Vz()/100.);
107
108 Float_t dt_scintillator = gRandom->Gaus(0,0.7);
109 Float_t time = dt_scintillator + 1e9*vhit->Tof();
110
111 if(vhit->Z() > 0){
112 if (time < 9 || time > 13) {continue;} //time window V0A : 11 +/- 2 ns
113
114 flagV0L = true;
115 hTimA->Fill(time);
116 if(time < timeV0L) timeV0L = time;
117 hMapV0L->Fill(vhit->X()/100.,vhit->Y()/100.);}
118
119 if(vhit->Z() < 0){
120 if (time < 1 ||time > 5) {continue;} //time window V0C : 3 +/- 2 ns
121
122 flagV0R = true;
123 hTimC->Fill(time);
124 if(time < timeV0R) timeV0R = time;
125 hMapV0R->Fill(vhit->X()/100.,vhit->Y()/100.);}
126
127 Int_t nPhot = vhit->Nphot();
128 Int_t cell = vhit->Cell();
129 map[cell] += nPhot;
130 hit[cell] ++;
131
132 }
133 }
134 }
135
136
137 Int_t map2[64]; // cell to digits
138 Int_t hit2[64]; // cell to digits
139 Int_t j;
140 Float_t time1[80];
141 Float_t time2[64];
142
143 for (j=0; j<16; j++){
144 map2[j] = map [j];
145 hit2[j] = hit [j];
146 //time2[j]=time1[j];
147 }
148
149 for (j=0; j<16; j++){
150 map2[j+16] = map [2*j+16]+map [2*j+17];
151 hit2[j+16] = hit [2*j+16]+hit [2*j+17];
152 //time2[j+16]= TMath::Min(time1 [16 + 2*j], time1 [16 + 2*j + 1]);
153 }
154
155 for (j=32; j<64; j++){
156 map2[j] =map[j+16];
157 hit2[j] =hit[j+16];
158 //time2[j] =time[j+16];
159 }
160
161 fNdigits = 0;
162 fMulC = 0;
163 fMulA = 0;
164 for (Int_t i=0; i<64; i++) {
165 Float_t q1 = Float_t ( map2[i] )* cPM[i] * kQe;
166 Float_t noise = gRandom->Gaus(10.5,3.22);
167 Float_t pmResponse = q1/kC*TMath::Power(ktheta/kthau,1/(1-ktheta/kthau))
168 + noise*1e-3;
169 map2[i] = Int_t( pmResponse * coef);
170 int test = 0;
171 if(map2[i] > 10) {map2[i] = Int_t( gRandom->Gaus(map2[i], 80));} // charge smearing of MIP/4 -> sigma = MIP/4 = 120 (MIP = 480)
172 if(map2[i] > 240) { // cut at MIP/2 = 240
173 hCell->Fill(float(i));
174 hCellADC->Fill(float(i),map2[i]);
175 fNdigits++;
176 if(i<32) fMulC++;
177 if(i>31) fMulA++;
178 }
179
180 }
181 hMul0->Fill(fNdigits);
182 hMulC0->Fill(fMulC);
183 hMulA0->Fill(fMulA);
184 hMulAnd0->Fill(TMath::Min(fMulA, fMulC));
185
186 if(fMulA > 0){
187 nVOL++;}
188
189 if(fMulC > 0){
190 nVOR++;}
191
192 if(fMulA > 0 && fMulC > 0){
193 nVOLetR++;}
194
195 if(fMulA > 0 || fMulC > 0){
196 nVOLouR++;}
197
198 if(event%100==0) cout <<" event = " << event <<endl;
199// cout <<" multi = " << fNdigits <<endl;
200 }
201
202 cout <<" nVOA = " << nVOL <<endl;
203 cout <<" nVOC = " << nVOR <<endl;
204 cout <<" nVOAandC = " << nVOLetR <<endl;
205 cout <<" nVOAorC = " << nVOLouR <<endl;
206
207//__________________________________________________
208// Fill root file
209
210 TFile *histoFile = new TFile("Efficiencies.root","RECREATE");
211
212 hMapV0L->Write();
213 hMapV0R->Write();
214 hpdgV0->Write();
215 hvertexZ->Write();
216 hTimC->Write();
217 hTimA->Write();
218 hCell->Write();
219 hCellADC->Write();
220
221 hMul0->Write();
222 hMulC0->Write();
223 hMulA0->Write();
224 hMulAnd0->Write();
225 histoFile->Close();
226
227}