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