Update (Zaida, Rossella)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / makeTFile4CutsLctopKpi.C
1 #include <Riostream.h>
2 #include <TFile.h>
3 #include <AliRDHFCutsLctopKpi.h>
4 #include <AliAODPidHF.h>
5 #include <TClonesArray.h>
6 #include <TParameter.h>
7
8
9 //Use:
10 //Set hard coded commentet with //set this!!
11 // root[] .L makeInput...C++
12 // root[] makeInputAliAnalysisTaskSE...()
13 //similar macros for the other D mesons
14
15 //Author: Rosa Romita, r.romita@gsi.de
16
17
18 //macro to make a .root file which contains an AliRDHFCutsLctopKpi for AliAnalysisTaskSELambdac task
19
20 void makeInputAliAnalysisTaskSELctopKpi(){
21
22   AliRDHFCutsLctopKpi* RDHFLctopKpiProd=new AliRDHFCutsLctopKpi();
23   RDHFLctopKpiProd->SetName("LctopKpiProdCuts");
24   RDHFLctopKpiProd->SetTitle("Production cuts for Lc analysis");
25
26   AliRDHFCutsLctopKpi* RDHFLctopKpiAn=new AliRDHFCutsLctopKpi();
27   RDHFLctopKpiAn->SetName("LctopKpiAnalysisCuts");
28   RDHFLctopKpiAn->SetTitle("Analysis cuts for Lc analysis");
29
30   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
31   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
32   //default
33   esdTrackCuts->SetRequireTPCRefit(kTRUE);
34   esdTrackCuts->SetRequireITSRefit(kTRUE);
35   esdTrackCuts->SetMinNClustersITS(4); // default is 5
36   esdTrackCuts->SetMinNClustersTPC(70);
37   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
38                                          AliESDtrackCuts::kAny); 
39  // default is kBoth, otherwise kAny
40   esdTrackCuts->SetMinDCAToVertexXY(0.);
41   esdTrackCuts->SetPtRange(0.3,1.e10);
42
43
44   RDHFLctopKpiProd->AddTrackCuts(esdTrackCuts);
45   RDHFLctopKpiAn->AddTrackCuts(esdTrackCuts);
46
47   const Int_t nvars=12;
48
49   const Int_t nptbins=4;
50   Float_t* ptbins;
51   ptbins=new Float_t[nptbins+1];
52   ptbins[0]=0.;
53   ptbins[1]=2.;
54   ptbins[2]=3.;
55   ptbins[3]=4.;
56   ptbins[4]=9999.;
57   
58
59   Float_t** prodcutsval;
60   prodcutsval=new Float_t*[nvars];
61   for(Int_t iv=0;iv<nvars;iv++){
62     prodcutsval[iv]=new Float_t[nptbins];
63   }
64
65   //  inv. mass [GeV]
66   // pTK [GeV/c]
67   // pTPi [GeV/c]
68   // d0K [cm]   lower limit!
69   // d0Pi [cm]  lower limit!
70   // dist12 (cm)
71   // sigmavert (cm)
72   // dist prim-sec (cm)
73   // pM=Max{pT1,pT2,pT3} (GeV/c)
74   // cosThetaPoint
75   // Sum d0^2 (cm^2)
76   // dca cut (cm)
77   for(Int_t ipt=0;ipt<nptbins;ipt++){
78     prodcutsval[0][ipt]=0.18;
79     prodcutsval[1][ipt]=0.4;
80     prodcutsval[2][ipt]=0.5;
81     prodcutsval[3][ipt]=0.;
82     prodcutsval[4][ipt]=0.;
83     prodcutsval[5][ipt]=0.01;
84     prodcutsval[6][ipt]=0.06;
85     prodcutsval[7][ipt]=0.005;
86     prodcutsval[8][ipt]=0.7;
87     prodcutsval[9][ipt]=0.;
88     prodcutsval[10][ipt]=0.;
89     prodcutsval[11][ipt]=0.05;
90   }
91
92   RDHFLctopKpiProd->SetPtBins(nptbins+1,ptbins);
93   RDHFLctopKpiProd->SetCuts(nvars,nptbins,prodcutsval);
94
95   Float_t** anacutsval;
96   anacutsval=new Float_t*[nvars];
97   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
98   for(Int_t ipt2=0;ipt2<nptbins;ipt2++){
99    anacutsval[0][ipt2]=0.18;
100    anacutsval[1][ipt2]=0.6;
101    anacutsval[3][ipt2]=0.;
102    anacutsval[4][ipt2]=0.;
103    anacutsval[5][ipt2]=0.01;
104    anacutsval[6][ipt2]=0.03;
105    anacutsval[9][ipt2]=0.;
106    anacutsval[10][ipt2]=0.;
107    anacutsval[11][ipt2]=0.05;
108   }
109
110   anacutsval[2][0]=0.6;
111   anacutsval[2][1]=0.8;
112   anacutsval[2][2]=1.;
113   anacutsval[2][3]=1.2;
114
115   anacutsval[7][0]=0.005;
116   anacutsval[7][1]=0.015;
117   anacutsval[7][2]=0.018;
118   anacutsval[7][3]=0.018;
119
120   anacutsval[8][0]=0.6;
121   anacutsval[8][1]=0.8;
122   anacutsval[8][2]=1.;
123   anacutsval[8][3]=1.2;
124
125   anacutsval[11][0]=0.04;
126   anacutsval[11][1]=0.04;
127   anacutsval[11][2]=0.03;
128   anacutsval[11][3]=0.03;
129
130
131   RDHFLctopKpiAn->SetPtBins(nptbins+1,ptbins);
132   RDHFLctopKpiAn->SetCuts(nvars,nptbins,anacutsval);
133
134  //  RDHFLc->SetRecoKF(); //set this if you want to recompute the secondary vertex with the KF package
135     //pid settings
136     //1. kaon: default one
137   AliAODPidHF* pidObjK=new AliAODPidHF();
138   Double_t sigmasK[5]={3.,1.,1.,3.,2.};
139   pidObjK->SetSigma(sigmasK);
140   pidObjK->SetAsym(kTRUE);
141   pidObjK->SetMatch(1);
142   pidObjK->SetTPC(kTRUE);
143   pidObjK->SetTOF(kTRUE);
144   pidObjK->SetITS(kTRUE);
145   Double_t plimK[2]={0.5,0.8};
146   pidObjK->SetPLimit(plimK,2);
147   
148   RDHFLctopKpiProd->SetPidHF(pidObjK);
149   RDHFLctopKpiAn->SetPidHF(pidObjK);
150
151     //2. pion 
152   AliAODPidHF* pidObjpi=new AliAODPidHF();
153   pidObjpi->SetTPC(kTRUE);
154   Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
155   pidObjpi->SetSigma(sigmaspi);
156
157   RDHFLctopKpiProd->SetPidpion(pidObjpi);
158   RDHFLctopKpiAn->SetPidpion(pidObjpi);
159
160   // 3. proton
161   AliAODPidHF* pidObjp=new AliAODPidHF();
162   Double_t sigmasp[5]={3.,1.,1.,3.,2.};
163   pidObjp->SetSigma(sigmasp);
164   pidObjp->SetAsym(kTRUE);
165   pidObjp->SetMatch(1);
166   pidObjp->SetTPC(kTRUE);
167   pidObjp->SetTOF(kTRUE);
168   pidObjp->SetITS(kTRUE);
169   Double_t plimp[2]={1.,2.};
170   pidObjp->SetPLimit(plimp,2);
171
172   RDHFLctopKpiProd->SetPidprot(pidObjp);
173   RDHFLctopKpiAn->SetPidprot(pidObjp);
174
175   Bool_t pidflag=kTRUE;
176   RDHFLctopKpiAn->SetUsePID(pidflag);
177   RDHFLctopKpiProd->SetUsePID(pidflag);
178   if(pidflag) cout<<"PID is used"<<endl;
179   else cout<<"PID is not used"<<endl;
180
181   cout<<"This is the object I'm going to save:"<<endl;
182   RDHFLctopKpiProd->PrintAll();
183   RDHFLctopKpiAn->PrintAll();
184   TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE"); 
185   fout->cd();
186   RDHFLctopKpiProd->Write();
187   RDHFLctopKpiAn->Write();
188   fout->Close();
189   delete fout;
190   delete pidObjp;
191   delete pidObjpi;
192   delete pidObjK;
193   delete RDHFLctopKpiProd;
194   delete RDHFLctopKpiAn;
195
196 }
197
198
199 //macro to make a .root file (for significance maximization) which contains an AliRDHFCutsLctopKpi with loose set of cuts  and TParameter with the tighest value of these cuts
200
201 void makeInputAliAnalysisTaskSESignificanceMaximization(){
202
203   AliRDHFCutsLctopKpi* RDHFLctopKpi=new AliRDHFCutsLctopKpi();
204   RDHFLctopKpi->SetName("loosercuts");
205   RDHFLctopKpi->SetTitle("Cuts for significance maximization");
206
207   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
208   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
209   //default
210   esdTrackCuts->SetRequireTPCRefit(kTRUE);
211   esdTrackCuts->SetMinNClustersTPC(70);
212   esdTrackCuts->SetRequireITSRefit(kTRUE);
213   esdTrackCuts->SetMinNClustersITS(4);
214
215   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
216   esdTrackCuts->SetMinDCAToVertexXY(0.);
217   esdTrackCuts->SetEtaRange(-0.8,0.8);
218   esdTrackCuts->SetPtRange(0.3,1.e10);
219   
220   RDHFLctopKpi->AddTrackCuts(esdTrackCuts);
221
222   const Int_t nvars=12;
223
224   const Int_t nptbins=4; //change this when adding pt bins!
225   Float_t ptbins[nptbins+1];
226   ptbins[0]=0.;
227   ptbins[1]=2.;
228   ptbins[2]=3.;
229   ptbins[3]=4.;
230   ptbins[4]=9999.;
231
232   RDHFLctopKpi->SetPtBins(nptbins+1,ptbins);
233
234   Float_t** rdcutsvalmine;
235   rdcutsvalmine=new Float_t*[nvars];
236   for(Int_t iv=0;iv<nvars;iv++){
237     rdcutsvalmine[iv]=new Float_t[nptbins];
238   }
239
240   //setting my cut values
241   //  inv. mass [GeV]
242   // pTK [GeV/c]
243   // pTPi [GeV/c]
244   // d0K [cm]   lower limit!
245   // d0Pi [cm]  lower limit!
246   // dist12 (cm)
247   // sigmavert (cm)
248   // dist prim-sec (cm)
249   // pM=Max{pT1,pT2,pT3} (GeV/c)
250   // cosThetaPoint
251   // Sum d0^2 (cm^2)
252   // dca cut (cm)
253   Float_t cutsMatrixLctopKpiStand[nptbins][nvars]=
254   {{0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
255    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
256    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05},
257    {0.18,0.4,0.5,0.,0.,0.01,0.06,0.005,0.7,0.,0.,0.05}};
258
259   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
260   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
261   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
262   for (Int_t ibin=0;ibin<nptbins;ibin++){
263     for (Int_t ivar = 0; ivar<nvars; ivar++){
264       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixLctopKpiStand[ibin][ivar];
265     }
266   }
267   RDHFLctopKpi->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
268
269
270   Int_t nvarsforopt=RDHFLctopKpi->GetNVarsForOpt();
271   Int_t dim=4; //set this!!
272   Bool_t *boolforopt;
273   boolforopt=new Bool_t[nvars];
274   if(dim>nvarsforopt){
275     cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
276     return;
277   } else {
278     if(dim==nvarsforopt){
279       boolforopt=RDHFLctopKpi->GetVarsForOpt();
280     }else{
281       TString *names;
282       names=new TString[nvars];
283       TString answer="";
284       Int_t checktrue=0;
285       names=RDHFLctopKpi->GetVarNames();
286       for(Int_t i=0;i<nvars;i++){
287         cout<<names[i]<<" for opt? (y/n)"<<endl;
288         cin>>answer;
289         if(answer=="y") {
290           boolforopt[i]=kTRUE;
291           checktrue++;
292         }
293         else boolforopt[i]=kFALSE;
294       }
295       if (checktrue!=dim) {
296         cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
297         return;
298       }
299       RDHFLctopKpi->SetVarsForOpt(dim,boolforopt);
300     }
301   }
302
303
304   Float_t tighterval[dim][nptbins];
305   // 5 dist12 (cm)        <---
306   // 7 dist prim-sec (cm) <---
307   // 9 cosThetaPoint      <---
308   // 10 Sum d0^2 (cm^2)   <---
309   // 11 dca cut (cm)      
310   // {0.18,0.4,0.5,0.,0.,0.01 <--- ,0.06,0.005 <--- ,0.7,0. <--- ,0. <--- ,0.05 }
311
312   //number of steps for each variable is set in the AddTask arguments (default=8)
313   //set this!!
314   for(Int_t ipt=0;ipt<nptbins;ipt++){
315     tighterval[0][ipt]=0.03;
316     tighterval[1][ipt]=0.02;
317     tighterval[2][ipt]=1.;
318     tighterval[3][ipt]=0.01;
319   }
320
321
322
323   TString name=""; 
324   Int_t arrdim=dim*nptbins;
325   cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
326   TClonesArray max("TParameter<float>",arrdim);
327   for(Int_t ival=0;ival<dim;ival++){
328     for(Int_t jpt=0;jpt<nptbins;jpt++){
329       name=Form("par%dptbin%d",ival,jpt);
330       cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
331       new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
332     }
333   }
334
335   Bool_t flagPID=kFALSE;
336   RDHFLctopKpi->SetUsePID(flagPID);
337
338   RDHFLctopKpi->PrintAll();
339   printf("Use PID? %s\n",flagPID ? "yes" : "no");
340
341     //pid settings
342     //1. kaon: default one
343   AliAODPidHF* pidObjK=new AliAODPidHF();
344   Double_t sigmasK[5]={3.,1.,1.,3.,2.};
345   pidObjK->SetSigma(sigmasK);
346   pidObjK->SetAsym(kTRUE);
347   pidObjK->SetMatch(1);
348   pidObjK->SetTPC(kTRUE);
349   pidObjK->SetTOF(kTRUE);
350   pidObjK->SetITS(kTRUE);
351   Double_t plimK[2]={0.5,0.8};
352   pidObjK->SetPLimit(plimK,2);
353   
354   RDHFLctopKpi->SetPidHF(pidObjK);
355
356     //2. pion 
357   AliAODPidHF* pidObjpi=new AliAODPidHF();
358   pidObjpi->SetTPC(kTRUE);
359   Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
360   pidObjpi->SetSigma(sigmaspi);
361
362   RDHFLctopKpi->SetPidpion(pidObjpi);
363
364   // 3. proton
365   AliAODPidHF* pidObjp=new AliAODPidHF();
366   Double_t sigmasp[5]={3.,1.,1.,3.,2.};
367   pidObjp->SetSigma(sigmasp);
368   pidObjp->SetAsym(kTRUE);
369   pidObjp->SetMatch(1);
370   pidObjp->SetTPC(kTRUE);
371   pidObjp->SetTOF(kTRUE);
372   pidObjp->SetITS(kTRUE);
373   Double_t plimp[2]={1.,2.};
374   pidObjp->SetPLimit(plimp,2);
375
376   RDHFLctopKpi->SetPidprot(pidObjp);
377
378   //activate pileup rejection (for pp)
379   //RDHFLctopKpi->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
380
381   //Do not recalculate the vertex
382   RDHFLctopKpi->SetRemoveDaughtersFromPrim(kFALSE); //activate for pp
383
384   TString cent="";
385   //centrality selection (Pb-Pb)
386   Float_t minc=20,maxc=80;
387   RDHFLctopKpi->SetMinCentrality(minc);
388   RDHFLctopKpi->SetMaxCentrality(maxc);
389   cent=Form("%.0f%.0f",minc,maxc);
390   //  RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentV0M); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
391   RDHFLctopKpi->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
392
393   //temporary
394   //  RDHFLctopKpi->SetFixRefs();
395
396   TFile* fout=new TFile(Form("cuts4SignifMaxim%s%s%sRecVtx%sPileupRej.root", RDHFLctopKpi->GetUseCentrality()==0 ? "pp" : "PbPb",cent.Data(),RDHFLctopKpi->GetIsPrimaryWithoutDaughters() ? "" : "No",RDHFLctopKpi->GetOptPileUp() ? "" : "No"),"recreate");   //set this!! 
397
398   fout->cd();
399   RDHFLctopKpi->Write();
400   max.Write();
401   fout->Close();
402  
403 }
404