Update
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / makeTFile4CutsD0toKpi.C
1 #include <Riostream.h>
2 #include <TFile.h>
3 #include <AliRDHFCutsD0toKpi.h>
4 #include <TClonesArray.h>
5 #include <TParameter.h>
6
7
8 //Use:
9 //Set hard coded commentet with //set this!!
10 // root[] .L makeInputD0tasks.C++
11 // root[] makeInputAliAnalysisTaskSED0Mass()
12 // root[] makeInputAliAnalysisTaskSESignificanceMaximization()
13 //similar macros for the other D mesons
14
15 //Author: Chiara Bianchin, cbianchi@pd.infn.it
16
17
18 //macro to make a .root file which contains an AliRDHFCutsD0toKpi for AliAnalysisTaskSED0Mass task
19
20 void makeInputAliAnalysisTaskSED0Mass(){
21
22   AliRDHFCutsD0toKpi* RDHFD0toKpi=new AliRDHFCutsD0toKpi();
23   RDHFD0toKpi->SetName("D0toKpiCuts");
24   RDHFD0toKpi->SetTitle("Cuts for D0 analysis");
25
26   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
27   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
28   //default
29   esdTrackCuts->SetRequireTPCRefit(kTRUE);
30   esdTrackCuts->SetRequireITSRefit(kTRUE);
31   esdTrackCuts->SetMinNClustersITS(4); // default is 5
32   //esdTrackCuts->SetMinNClustersTPC(70);
33   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
34                                          AliESDtrackCuts::kAny); 
35  // default is kBoth, otherwise kAny
36   esdTrackCuts->SetMinDCAToVertexXY(0.);
37   esdTrackCuts->SetPtRange(0.3,1.e10);
38
39
40   RDHFD0toKpi->AddTrackCuts(esdTrackCuts);
41
42   const Int_t nvars=9;
43
44   const Int_t nptbins=5;
45   Float_t* ptbins;
46   ptbins=new Float_t[nptbins+1];
47   ptbins[0]=0.;
48   ptbins[1]=1.;
49   ptbins[2]=2.;
50   ptbins[3]=3.;
51   ptbins[4]=5.;
52   ptbins[5]=10.;
53   
54   RDHFD0toKpi->SetPtBins(nptbins+1,ptbins);
55   
56
57   Float_t** rdcutsvalmine;
58   rdcutsvalmine=new Float_t*[nvars];
59   for(Int_t iv=0;iv<nvars;iv++){
60     rdcutsvalmine[iv]=new Float_t[nptbins];
61   }
62
63   //setting my cut values
64     //cuts order
65     //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
66     //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
67     //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
68     //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
69     //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
70     //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
71     //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
72     //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
73     //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
74
75   /*
76   //setting PPR cut values
77   rdcutsvalPPR[0][0]=0.7;
78   rdcutsvalPPR[1][0]=0.04;
79   rdcutsvalPPR[2][0]=0.8;
80   rdcutsvalPPR[3][0]=0.5;
81   rdcutsvalPPR[4][0]=0.5;
82   rdcutsvalPPR[5][0]=0.05;
83   rdcutsvalPPR[6][0]=0.05;
84   rdcutsvalPPR[7][0]=-0.0002;
85   rdcutsvalPPR[8][0]=0.5;
86
87   rdcutsvalPPR[0][1]=rdcutsvalPPR[0][2]=0.7;
88   rdcutsvalPPR[1][1]=rdcutsvalPPR[1][2]=0.02;
89   rdcutsvalPPR[2][1]=rdcutsvalPPR[2][2]=0.8;
90   rdcutsvalPPR[3][1]=rdcutsvalPPR[3][2]=0.7;
91   rdcutsvalPPR[4][1]=rdcutsvalPPR[4][2]=0.7;
92   rdcutsvalPPR[5][1]=rdcutsvalPPR[5][2]=0.05;
93   rdcutsvalPPR[6][1]=rdcutsvalPPR[6][2]=0.05;
94   rdcutsvalPPR[7][1]=rdcutsvalPPR[7][2]=-0.0002;
95   rdcutsvalPPR[8][1]=rdcutsvalPPR[8][2]=0.6;
96
97   rdcutsvalPPR[0][3]=0.7;
98   rdcutsvalPPR[1][3]=0.02;
99   rdcutsvalPPR[2][3]=0.8;
100   rdcutsvalPPR[3][3]=0.7;
101   rdcutsvalPPR[4][3]=0.7;
102   rdcutsvalPPR[5][3]=0.05;
103   rdcutsvalPPR[6][3]=0.05;
104   rdcutsvalPPR[7][3]=-0.0001;
105   rdcutsvalPPR[8][3]=0.8;
106
107   rdcutsvalPPR[0][4]=0.7;
108   rdcutsvalPPR[1][4]=0.02;
109   rdcutsvalPPR[2][4]=0.8;
110   rdcutsvalPPR[3][4]=0.7;
111   rdcutsvalPPR[4][4]=0.7;
112   rdcutsvalPPR[5][4]=0.05;
113   rdcutsvalPPR[6][4]=0.05;
114   rdcutsvalPPR[7][4]=-0.00005;
115   rdcutsvalPPR[8][4]=0.8;
116   */
117
118   //setting my cut values
119   rdcutsvalmine[0][0]=0.7;
120   rdcutsvalmine[1][0]=0.04;
121   rdcutsvalmine[2][0]=0.8;
122   rdcutsvalmine[3][0]=0.5;
123   rdcutsvalmine[4][0]=0.5;
124   rdcutsvalmine[5][0]=0.05;
125   rdcutsvalmine[6][0]=0.05;
126   rdcutsvalmine[7][0]=-0.00025;
127   rdcutsvalmine[8][0]=0.7;
128
129   rdcutsvalmine[0][1]=rdcutsvalmine[0][2]=0.7;
130   rdcutsvalmine[1][1]=rdcutsvalmine[1][2]=0.02;
131   rdcutsvalmine[2][1]=rdcutsvalmine[2][2]=0.8;
132   rdcutsvalmine[3][1]=rdcutsvalmine[3][2]=0.7;
133   rdcutsvalmine[4][1]=rdcutsvalmine[4][2]=0.7;
134   rdcutsvalmine[5][1]=rdcutsvalmine[5][2]=1.;
135   rdcutsvalmine[6][1]=rdcutsvalmine[6][2]=1.;
136   rdcutsvalmine[7][1]=rdcutsvalmine[7][2]=-0.00025;
137   rdcutsvalmine[8][1]=rdcutsvalmine[8][2]=0.8;
138
139   rdcutsvalmine[0][3]=0.7;
140   rdcutsvalmine[1][3]=0.02;
141   rdcutsvalmine[2][3]=0.8;
142   rdcutsvalmine[3][3]=0.7;
143   rdcutsvalmine[4][3]=0.7;
144   rdcutsvalmine[5][3]=0.05;
145   rdcutsvalmine[6][3]=0.05;
146   rdcutsvalmine[7][3]=-0.00015;
147   rdcutsvalmine[8][3]=0.8;
148
149   rdcutsvalmine[0][4]=0.7;
150   rdcutsvalmine[1][4]=0.02;
151   rdcutsvalmine[2][4]=0.8;
152   rdcutsvalmine[3][4]=0.7;
153   rdcutsvalmine[4][4]=0.7;
154   rdcutsvalmine[5][4]=0.05;
155   rdcutsvalmine[6][4]=0.05;
156   rdcutsvalmine[7][4]=-0.00015;
157   rdcutsvalmine[8][4]=0.9;
158
159   RDHFD0toKpi->SetCuts(nvars,nptbins,rdcutsvalmine);
160   cout<<"This is the odject I'm going to save:"<<endl;
161   RDHFD0toKpi->PrintAll();
162   TFile* fout=new TFile("D0toKpiCuts.root","recreate");   //set this!! 
163   fout->cd();
164   RDHFD0toKpi->Write();
165   fout->Close();
166
167 }
168
169 //macro to make a .root file (for significance maximization) which contains an AliRDHFCutsD0toKpi with loose set of cuts  and TParameter with the tighest value of these cuts
170
171 void makeInputAliAnalysisTaskSESignificanceMaximization(){
172
173   AliRDHFCutsD0toKpi* RDHFD0toKpi=new AliRDHFCutsD0toKpi();
174   RDHFD0toKpi->SetName("loosercuts");
175   RDHFD0toKpi->SetTitle("Cuts for significance maximization");
176
177   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
178   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
179   //default
180   esdTrackCuts->SetRequireTPCRefit(kTRUE);
181   esdTrackCuts->SetRequireITSRefit(kTRUE);
182   esdTrackCuts->SetMinNClustersITS(4);
183   
184   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
185   esdTrackCuts->SetMinDCAToVertexXY(0.);
186   esdTrackCuts->SetEtaRange(-0.9,0.9);
187   esdTrackCuts->SetPtRange(0.1,1.e10);
188   
189   RDHFD0toKpi->AddTrackCuts(esdTrackCuts);
190
191   const Int_t nvars=9;
192
193   const Int_t nptbins=5;
194   Float_t ptbins[nptbins+1];
195   ptbins[0]=0.;
196   ptbins[1]=1.;
197   ptbins[2]=2.;
198   ptbins[3]=3.;
199   ptbins[4]=5.;
200   ptbins[5]=10.;
201   
202   RDHFD0toKpi->SetPtBins(nptbins+1,ptbins);
203   
204
205   Float_t** rdcutsvalmine;
206   rdcutsvalmine=new Float_t*[nvars];
207   for(Int_t iv=0;iv<nvars;iv++){
208     rdcutsvalmine[iv]=new Float_t[nptbins];
209   }
210
211   //setting my cut values
212     //cuts order
213     //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
214     //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
215     //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
216     //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
217     //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
218     //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
219     //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
220     //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
221     //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
222
223   /*
224   //setting PPR cut values
225   rdcutsvalPPR[0][0]=0.7;
226   rdcutsvalPPR[1][0]=0.04;
227   rdcutsvalPPR[2][0]=0.8;
228   rdcutsvalPPR[3][0]=0.5;
229   rdcutsvalPPR[4][0]=0.5;
230   rdcutsvalPPR[5][0]=0.05;
231   rdcutsvalPPR[6][0]=0.05;
232   rdcutsvalPPR[7][0]=-0.0002;
233   rdcutsvalPPR[8][0]=0.5;
234
235   rdcutsvalPPR[0][1]=rdcutsvalPPR[0][2]=0.7;
236   rdcutsvalPPR[1][1]=rdcutsvalPPR[1][2]=0.02;
237   rdcutsvalPPR[2][1]=rdcutsvalPPR[2][2]=0.8;
238   rdcutsvalPPR[3][1]=rdcutsvalPPR[3][2]=0.7;
239   rdcutsvalPPR[4][1]=rdcutsvalPPR[4][2]=0.7;
240   rdcutsvalPPR[5][1]=rdcutsvalPPR[5][2]=0.05;
241   rdcutsvalPPR[6][1]=rdcutsvalPPR[6][2]=0.05;
242   rdcutsvalPPR[7][1]=rdcutsvalPPR[7][2]=-0.0002;
243   rdcutsvalPPR[8][1]=rdcutsvalPPR[8][2]=0.6;
244
245   rdcutsvalPPR[0][3]=0.7;
246   rdcutsvalPPR[1][3]=0.02;
247   rdcutsvalPPR[2][3]=0.8;
248   rdcutsvalPPR[3][3]=0.7;
249   rdcutsvalPPR[4][3]=0.7;
250   rdcutsvalPPR[5][3]=0.05;
251   rdcutsvalPPR[6][3]=0.05;
252   rdcutsvalPPR[7][3]=-0.0001;
253   rdcutsvalPPR[8][3]=0.8;
254
255   rdcutsvalPPR[0][4]=0.7;
256   rdcutsvalPPR[1][4]=0.02;
257   rdcutsvalPPR[2][4]=0.8;
258   rdcutsvalPPR[3][4]=0.7;
259   rdcutsvalPPR[4][4]=0.7;
260   rdcutsvalPPR[5][4]=0.05;
261   rdcutsvalPPR[6][4]=0.05;
262   rdcutsvalPPR[7][4]=-0.00005;
263   rdcutsvalPPR[8][4]=0.8;
264   */
265   //setting my cut values
266
267   rdcutsvalmine[0][0]=0.065;
268   rdcutsvalmine[1][0]=0.04;
269   rdcutsvalmine[2][0]=0.8;
270   rdcutsvalmine[3][0]=0.7;
271   rdcutsvalmine[4][0]=0.7;
272   rdcutsvalmine[5][0]=1000.*1E-4;
273   rdcutsvalmine[6][0]=1000.*1E-4;
274   rdcutsvalmine[7][0]=-0.0004;
275   rdcutsvalmine[8][0]=0.71;
276
277   rdcutsvalmine[0][1]=0.065;
278   rdcutsvalmine[1][1]=0.04;
279   rdcutsvalmine[2][1]=0.8;
280   rdcutsvalmine[3][1]=0.7;
281   rdcutsvalmine[4][1]=0.7;
282   rdcutsvalmine[5][1]=1.;
283   rdcutsvalmine[6][1]=1.;
284   rdcutsvalmine[7][1]=-0.0003;
285   rdcutsvalmine[8][1]=0.79;
286
287   rdcutsvalmine[0][2]=0.65;
288   rdcutsvalmine[1][2]=200.*1E-4;
289   rdcutsvalmine[2][2]=0.8;
290   rdcutsvalmine[3][2]=0.7;
291   rdcutsvalmine[4][2]=0.7;
292   rdcutsvalmine[5][2]=1000.*1E-4;
293   rdcutsvalmine[6][2]=1000.*1E-4;
294   rdcutsvalmine[7][2]=-0.0001;
295   rdcutsvalmine[8][2]=0.83;
296
297   rdcutsvalmine[0][3]=0.65;
298   rdcutsvalmine[1][3]=200.*1E-4;
299   rdcutsvalmine[2][3]=0.8;
300   rdcutsvalmine[3][3]=0.7;
301   rdcutsvalmine[4][3]=0.7;
302   rdcutsvalmine[5][3]=1000.*1E-4;
303   rdcutsvalmine[6][3]=1000.*1E-4;
304   rdcutsvalmine[7][3]=-0.00005;
305   rdcutsvalmine[8][3]=0.78;
306
307   rdcutsvalmine[0][4]=0.65;
308   rdcutsvalmine[1][4]=200.*1E-4;
309   rdcutsvalmine[2][4]=0.8;
310   rdcutsvalmine[3][4]=0.7;
311   rdcutsvalmine[4][4]=0.7;
312   rdcutsvalmine[5][4]=1000.*1E-4;
313   rdcutsvalmine[6][4]=1000.*1E-4;
314   rdcutsvalmine[7][4]=-0.00001;
315   rdcutsvalmine[8][4]=0.79;
316
317
318   RDHFD0toKpi->SetCuts(nvars,nptbins,rdcutsvalmine);
319
320   Int_t nvarsforopt=RDHFD0toKpi->GetNVarsForOpt();
321   Int_t dim=2; //set this!!
322   Bool_t *boolforopt;
323   boolforopt=new Bool_t[nvars];
324   if(dim>nvarsforopt){
325     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
326     return;
327   } else {
328     if(dim==nvarsforopt){
329       boolforopt=RDHFD0toKpi->GetVarsForOpt();
330     }else{
331       TString *names;
332       names=new TString[nvars];
333       TString answer="";
334       Int_t checktrue=0;
335       names=RDHFD0toKpi->GetVarNames();
336       for(Int_t i=0;i<nvars;i++){
337         cout<<names[i]<<" for opt? (y/n)"<<endl;
338         cin>>answer;
339         if(answer=="y") {
340           boolforopt[i]=kTRUE;
341           checktrue++;
342         }
343         else boolforopt[i]=kFALSE;
344       }
345       if (checktrue!=dim) {
346         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
347         return;
348       }
349       RDHFD0toKpi->SetVarsForOpt(dim,boolforopt);
350     }
351   }
352
353
354   Float_t tighterval[dim][nptbins];
355   //dca
356   //costhetastar
357   //d0d0 <-this 
358   //costhetapoint <-this 
359
360   
361   //number of steps for each variable is 4 now
362   //set this!!
363   // tighterval[0][0]=0.01;
364   // tighterval[1][0]=0.8;
365   tighterval[0][0]=-0.0007;
366   tighterval[1][0]=0.99;
367
368   // tighterval[0][1]=0.01;
369   // tighterval[1][1]=0.8;
370   tighterval[0][1]=-0.0006;
371   tighterval[1][1]=0.99;
372  
373  // tighterval[0][2]=0.01;
374   // tighterval[1][2]=0.8;
375   tighterval[0][2]=-0.0004;
376   tighterval[1][2]=0.99;
377  
378   // tighterval[0][3]=0.01;
379   // tighterval[1][3]=0.8;
380   tighterval[0][3]=-0.00035;
381   tighterval[1][3]=0.98;
382
383   // tighterval[0][4]=0.01;
384   // tighterval[1][4]=0.8;
385   tighterval[0][4]=-0.0003;
386   tighterval[1][4]=0.99;
387
388
389   TString name=""; 
390   Int_t arrdim=dim*nptbins;
391   cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
392   TClonesArray max("TParameter<float>",arrdim);
393   for(Int_t ival=0;ival<dim;ival++){
394     for(Int_t jpt=0;jpt<nptbins;jpt++){
395       name=Form("par%dptbin%d",ival,jpt);
396       cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
397       new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
398     }
399   }
400  
401   TFile* fout=new TFile("cuts4SignifMaxim.root","recreate");   //set this!! 
402   fout->cd();
403   RDHFD0toKpi->Write();
404   max.Write();
405   fout->Close();
406  
407 }
408