Faster reweighting method implemented.
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTriggerEff.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 #include "AliFastMuonTriggerEff.h"
19 #include <TFile.h>
20 #include <TTree.h>
21 #include <TROOT.h>
22 #include <stdlib.h>
23
24 // Debugging flag
25 //#define MYTRIGDEBUG
26
27 ClassImp(AliFastMuonTriggerEff);
28
29 AliFastMuonTriggerEff::AliFastMuonTriggerEff():
30     AliFastResponse("Efficiency", "Muon Trigger Efficiency")
31 {
32 //
33 // Default constructor
34 //
35     SetCut(kLow);
36     fZones=0;
37 }
38
39
40 void AliFastMuonTriggerEff::Init()
41 {
42 //
43 //  Initialization
44 //
45     printf("Initializing %s / %s", fName.Data(), fTitle.Data());
46 //
47
48     fPhiMin   = -90. ;
49     fPhiMax   =  90. ;
50     fDphi     =   9. ;
51     fThetaMin =   2. ;
52     fThetaMax =   9. ;
53     fDtheta   =   0.7;
54     fDpt      =   0.1;
55
56     InitTree();
57 }
58
59 void AliFastMuonTriggerEff::InitTree()
60 {
61 //
62 //  Initialize tables from tree
63 //
64     TTree          *chain;   // pointer to the analyzed TTree or TChain
65     Int_t           nSim, nZona;
66     Double_t        pmin, pmax, tmin, tmax, phimin, phimax, bkg;
67     Double_t        len50, hen50, leff, heff;
68     Double_t        vLPt[50];
69     Double_t        vHPt[50];
70     Char_t file[100]="$(ALICE_ROOT)/FASTSIM/data/vettorpara.root";
71 //
72 //  Avoid memory leak in case of reinitialization
73     if(fZones!=0) {
74         printf("\nWarning: reinitialization of an object of class: AliFastMuonTriggerEff\n");
75         for (Int_t i=0; i<fZones; i++) {
76             if(fEffLow [i])delete[] fEffLow[i];
77             if(fEffHigh[i])delete[] fEffHigh[i];
78         }
79         if(fEffLow) {
80             delete[] fEffLow;
81             fEffLow=0;
82         }
83         if(fEffHigh) {
84             delete[] fEffHigh;
85             fEffHigh=0;
86         }
87     }
88     printf("AliFastMuonTriggerEff: Initialization\n");
89     TFile *f = new TFile(file);
90     if(f->IsZombie()) {
91         printf("Cannot open file: %s\n",file);
92         return;
93     }
94     f->ls();
95     
96     TTree* tree = (TTree*)gDirectory->Get("fitPar");
97
98 //
99 //
100     chain = tree;
101     chain->SetMakeClass(1);
102     
103     chain->SetBranchAddress("nSim",&nSim);
104     chain->SetBranchAddress("nZona",&nZona);
105     chain->SetBranchAddress("ptmin",&pmin);
106     chain->SetBranchAddress("ptmax",&pmax);
107     chain->SetBranchAddress("Thetamin",&tmin);
108     chain->SetBranchAddress("Thetamax",&tmax);
109     chain->SetBranchAddress("Phimin",&phimin);
110     chain->SetBranchAddress("Phimax",&phimax);
111     chain->SetBranchAddress("Bkg",&bkg);
112     chain->SetBranchAddress("EffLPt",vLPt);
113     chain->SetBranchAddress("EffHPt",vHPt);
114     chain->SetBranchAddress("Pt0.5Low",&len50);
115     chain->SetBranchAddress("Pt0.5High",&hen50);
116     chain->SetBranchAddress("EffLowMax",&leff);
117     chain->SetBranchAddress("EffHighMax",&heff);
118 // 
119 //
120     Int_t nentries = Int_t(chain->GetEntries());
121     Int_t nbytes = 0, nb = 0;
122
123 // Count the number of zones of the parametrization
124     Int_t nzone0=0, nzone1=0;
125     for (Int_t jentry=0; jentry<nentries; jentry++) {
126         nb = chain->GetEntry(jentry);
127         if(bkg==0.) {
128              if(nSim==0)nzone0++;
129              if(nSim==1)nzone1++;
130         }
131     }
132     
133     printf("Trigger parametrisation for %d zones for Pt: 0. 5. GeV/c\n",nzone0);
134     printf("and %d zones extended to 10. GeV/c\n",nzone1);
135     fZones=nzone0+nzone1;
136 //    printf("Ciao\n");
137     if(fZones<=0){
138         printf("Number of zones must be greater than 0\n");
139         exit(6);
140     }
141     
142     fEffLow =new Float_t*[fZones];
143     fEffHigh=new Float_t*[fZones];
144     for (Int_t i=0; i<fZones; i++) {
145         fEffLow [i]=new Float_t[50];
146         fEffHigh[i]=new Float_t[50];
147     }
148     
149 //  Initialize look-up table to standard values
150     Int_t isim, itheta, iphi;
151     for (isim=0; isim<2; isim++) {
152         for (itheta=0; itheta<10; itheta++) {
153             for (iphi=0; iphi<20; iphi++) {
154                 fLook[isim][itheta][iphi]=0;
155             }
156         }
157     }
158
159 //  Loading Trigger efficiencies
160     Int_t myzone=0;
161 #ifdef MYTRIGDEBUG
162             printf("Filling nSim nZona pmin pmax tmin tmax phimin phimax: ....\n");
163 #endif
164     for (Int_t jentry=0; jentry<nentries; jentry++) {
165 //      Int_t ientry = LoadTree(jentry); 
166         nb = chain->GetEntry(jentry);   
167         nbytes += nb;
168 #ifdef MYTRIGDEBUG
169        printf("Getting entry %d... ",jentry);
170 #endif
171 // For the time being it works with background 0
172         if ((nSim == 0 || nSim == 1)&&bkg==0.) {
173 #ifdef MYTRIGDEBUG
174             printf("Filling %d %d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f: ",
175             nSim,nZona,pmin,pmax,tmin,tmax,phimin,phimax);
176 #endif
177             for (Int_t k = 0; k < 50; k++) {
178                 fEffLow [myzone][k] = vLPt[k];
179                 fEffHigh[myzone][k] = vHPt[k];
180 #ifdef MYTRIGDEBUG
181                 if(k<15)printf(" %5.3f",vLPt[k]);
182 #endif
183             }
184 #ifdef MYTRIGDEBUG
185             printf("\n");
186 #endif
187             myzone++;
188             iphi=Int_t(nZona/10)-10;
189             itheta=nZona%10;
190             if(iphi<0||iphi>19||itheta<0||itheta>9) {
191                 printf("The format of data file is not consistent\nwith this version of the code\n");
192                 printf("This should never happen: iphi %d, itheta: %d\n",iphi,itheta);
193                 exit(7);
194             }
195             fLook[nSim][itheta][iphi]=myzone;
196         } else {
197             printf("Skipping entry with nSim=%d, bkg=%f\n",nSim,bkg);
198         }
199     }
200 #ifdef MYTRIGDEBUG
201     printf("This is the content of the LUT after first step of initialization\n");
202     for(isim=0; isim<2; isim++) {
203         printf("isim=%d\n",isim);
204         for(iphi=0; iphi<20; iphi++) {
205             printf("iphi: %2d:",iphi);
206             for(itheta=0; itheta<10; itheta++) {
207                 printf(" %4d",fLook[isim][itheta][iphi]);
208             }
209             printf("\n");
210         }
211     }
212 #endif
213 // Filling look=up table for the zones where the extended simulation does
214 // not exists
215     for(iphi=0; iphi<20; iphi++) {
216         for(itheta=0; itheta<10; itheta++) {
217             if(fLook[0][itheta][iphi]==0) {
218                 printf("Missing entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
219                 exit(8);
220             }
221             if(fLook[0][itheta][iphi]<0||fLook[0][itheta][iphi]>fZones) {
222                 printf("Problem with entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
223                 exit(9);
224             }
225             if(fLook[1][itheta][iphi]==0) {
226                  fLook[1][itheta][iphi]=-fLook[0][itheta][iphi];
227             }
228         }
229     }
230 #ifdef MYTRIGDEBUG
231     for(isim=0; isim<2; isim++) {
232         printf("isim=%d\n",isim);
233         for(iphi=0; iphi<20; iphi++) {
234             printf("iphi: %2d:",iphi);
235             for(itheta=0; itheta<10; itheta++) {
236                 printf(" %4d",fLook[isim][itheta][iphi]);
237             }
238             printf("\n");
239         }
240     }
241         for(iphi=0; iphi<20; iphi++) {
242             for(itheta=0; itheta<10; itheta++) {
243     for(isim=0; isim<2; isim++) {
244                 printf("%d %2d %2d %4d:",isim,iphi,itheta,fLook[isim][itheta][iphi]);
245                 if(fLook[isim][itheta][iphi]>0) {
246                     myzone=fLook[isim][itheta][iphi]-1;
247                     for(Int_t ipt=0; ipt<20; ipt++) {
248                         //printf(" %5.3f",fEffLow[myzone][ipt]);
249                         printf(" %5.3f",fEffHigh[myzone][ipt]);
250                     }
251                     printf(" ...");
252                     for(Int_t ipt=40; ipt<50; ipt++) {
253                         //printf(" %5.3f",fEffLow[myzone][ipt]);
254                         printf(" %5.3f",fEffHigh[myzone][ipt]);
255                     }               
256                 }
257                 printf("\n");
258             }
259         }
260     }    
261 #endif
262     f->Close();
263 }
264
265 void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
266                 Float_t theta, Float_t phi,Float_t& effLow, Float_t& effHigh)
267 {
268 //
269 //  Trigger efficiency for pt, theta, phi (low and high cut)
270 //
271     effLow=0.;
272     effHigh=0.;
273     if(fZones==0) {
274         printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
275         return;
276     }
277     if(pt<0) {
278         printf("Warning: pt: %f < 0. GeV/c\n",pt);
279         return; 
280     }
281     Int_t iPt   = Int_t(pt/fDpt)%50;
282     Int_t iSim  = Int_t(pt/fDpt)/50;
283     Int_t iPhi  = Int_t((phi-fPhiMin)/fDphi);
284     if(phi<fPhiMin)iPhi=iPhi-1;
285     Int_t iTheta = Int_t((theta-fThetaMin)/fDtheta);
286 #ifdef MYTRIGDEBUG
287     printf("iSim iPt iTheta iPhi: %d %d %d %d\n",iSim,iPt,iTheta,iPhi);
288 #endif
289     iPhi=iPhi-40*(iPhi/40);
290 #ifdef MYTRIGDEBUG
291     printf("1: iPhi converted to: %d for angle equivalence\n",iPhi);
292 #endif
293     if(iPhi<0)iPhi=-iPhi-1;
294     if(iPhi>19)iPhi=39-iPhi;
295 #ifdef MYTRIGDEBUG
296     printf("2: iPhi converted to: %d for the symmetry of the spectrometer\n",iPhi);
297 #endif
298     if(charge==1.){
299     } else if(charge==-1.) {
300     iPhi=19-iPhi;
301 #ifdef MYTRIGDEBUG
302     printf("3: iPhi converted to: %d for charge symmetry\n",iPhi);
303 #endif
304     } else {
305         printf("Warning: not understand charge: %f\n",charge);
306         return;
307     }
308     if(iTheta<0||iTheta>9) {
309         printf("Warning: theta: %f outside acceptance\n",theta);
310         return;
311     }
312     if(iPt<0) {
313         printf("Warning: what do you mean with pt: %f <0?\n",pt);
314         return;
315     }
316     if(iSim>=fSim) {
317         iSim=fSim-1;
318         iPt=49;
319 #ifdef MYTRIGDEBUG
320     printf("4: iSim iPt converted to: %d %d (last zone)\n",iSim,iPt);
321 #endif
322     }
323     Int_t iLook=fLook[iSim][iTheta][iPhi];
324     if(iLook<0) {
325 #ifdef MYTRIGDEBUG
326     printf("5: iLook iPt: %d %d converted to: ",iLook,iPt);
327 #endif
328         iLook=-iLook-1;
329         iPt=49;
330 #ifdef MYTRIGDEBUG
331     printf("%d %d from look up table contents\n",iLook,iPt);
332 #endif
333     } else {
334         iLook=iLook-1;
335     }
336     effLow=fEffLow[iLook][iPt];
337     effHigh=fEffHigh[iLook][iPt];
338 #ifdef MYTRIGDEBUG
339     printf("6: charge, iSim, iTheta, iPhi, iPt: %f %d %d %d %d effLow: %f, effHigh: %f\n",
340                charge,iSim,iTheta,iPhi,iPt,effLow,effHigh);
341 #endif
342     
343     //fEffLow [iPhi][iTheta][iPt];
344     //fEffHigh[iPhi][iTheta][iPt];    
345 }
346
347 Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
348                    Float_t theta, Float_t phi)
349 {
350     //
351     // Trigger efficiency for pt, theta, phi depending of fCut
352     // 
353     if(fZones==0) {
354         printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
355         return 0.;
356     }
357     Float_t eff;
358     Float_t effLow, effHigh;
359     
360     Evaluate(charge,pt,theta,phi,effLow,effHigh);
361     if (fCut == kLow) 
362         eff  = effLow;
363     else
364         eff  = effHigh;
365
366     return eff;
367 }