Fix in the last caall to CleanOwnPrimaryVertex
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsD0toKpipipi.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 /////////////////////////////////////////////////////////////\r
17 //\r
18 // Class for cuts on AOD reconstructed D0->Kpipipi\r
19 //\r
20 // Author: r.romita@gsi.de, andrea.dainese@pd.infn.it,
21 //         fabio.colamaria@ba.infn.it\r
22 /////////////////////////////////////////////////////////////\r
23 \r
24 #include <TDatabasePDG.h>\r
25 #include <Riostream.h>\r
26 \r
27 #include "AliRDHFCutsD0toKpipipi.h"\r
28 #include "AliAODRecoDecayHF4Prong.h"\r
29 #include "AliAODTrack.h"\r
30 #include "AliESDtrack.h"\r
31 #include "AliAODPidHF.h"
32 \r
33 ClassImp(AliRDHFCutsD0toKpipipi)\r
34 \r
35 //--------------------------------------------------------------------------\r
36 AliRDHFCutsD0toKpipipi::AliRDHFCutsD0toKpipipi(const char* name) : \r
37 AliRDHFCuts(name)\r
38 {\r
39   //\r
40   // Default Constructor\r
41   //\r
42   Int_t nvars=9;\r
43   SetNVars(nvars);\r
44   TString varNames[9]={"inv. mass [GeV]",   \r
45                        "dca [cm]",\r
46                        "Dist 2-trk Vtx to PrimVtx [cm]",\r
47                        "Dist 3-trk Vtx to PrimVtx [cm]",\r
48                        "Dist 4-trk Vtx to PrimVtx [cm]",\r
49                        "cosThetaPoint",\r
50                        "pt [GeV/c]",\r
51                        "rho mass [GeV]",\r
52                        "PID cut"};\r
53   Bool_t isUpperCut[9]={kTRUE,\r
54                         kTRUE,\r
55                         kFALSE,\r
56                         kFALSE,\r
57                         kFALSE,\r
58                         kFALSE,\r
59                         kFALSE,\r
60                         kTRUE,\r
61                         kFALSE};\r
62   SetVarNames(nvars,varNames,isUpperCut);\r
63   Bool_t forOpt[9]={kFALSE,\r
64                     kTRUE,\r
65                     kTRUE,\r
66                     kTRUE,\r
67                     kTRUE,\r
68                     kTRUE,\r
69                     kFALSE,\r
70                     kFALSE,\r
71                     kFALSE};\r
72   SetVarsForOpt(5,forOpt);\r
73   Float_t limits[2]={0,999999999.};\r
74   SetPtBins(2,limits);\r
75 }\r
76 //--------------------------------------------------------------------------\r
77 AliRDHFCutsD0toKpipipi::AliRDHFCutsD0toKpipipi(const AliRDHFCutsD0toKpipipi &source) :\r
78   AliRDHFCuts(source)\r
79 {\r
80   //\r
81   // Copy constructor\r
82   //\r
83 \r
84 }\r
85 //--------------------------------------------------------------------------\r
86 AliRDHFCutsD0toKpipipi &AliRDHFCutsD0toKpipipi::operator=(const AliRDHFCutsD0toKpipipi &source)\r
87 {\r
88   //\r
89   // assignment operator\r
90   //\r
91   if(&source == this) return *this;\r
92 \r
93   AliRDHFCuts::operator=(source);\r
94 \r
95   return *this;\r
96 }\r
97 \r
98 \r
99 //---------------------------------------------------------------------------\r
100 void AliRDHFCutsD0toKpipipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent* aod) {\r
101   // \r
102   // Fills in vars the values of the variables \r
103   //\r
104 \r
105   if(nvars!=fnVarsForOpt) {\r
106     printf("AliRDHFCutsD0toKpipipi::GetCutsVarsForOpt: wrong number of variables\n");\r
107     return;\r
108   }\r
109 \r
110   AliAODRecoDecayHF4Prong *dd = (AliAODRecoDecayHF4Prong*)d;\r
111 \r
112   //recalculate vertex w/o daughters
113   Bool_t cleanvtx=kFALSE;
114   AliAODVertex *origownvtx=0x0;
115   if(fRemoveDaughtersFromPrimary) {
116     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
117     cleanvtx=kTRUE;
118     if(!RecalcOwnPrimaryVtx(dd,aod)) {
119       CleanOwnPrimaryVtx(dd,aod,origownvtx);
120       cleanvtx=kFALSE;
121     }
122   }
123
124   Int_t iter=-1;\r
125 \r
126   if(fVarsForOpt[0]) {\r
127     iter++;\r
128     Double_t mD0[2],mD0bar[2];\r
129     if(TMath::Abs(pdgdaughters[1])==321 || TMath::Abs(pdgdaughters[3])==321) {\r
130       dd->InvMassD0(mD0);\r
131       if(TMath::Abs(pdgdaughters[1])==321) {\r
132        vars[iter]=mD0[0];\r
133       }else{\r
134        vars[iter]=mD0[1];\r
135       }\r
136     } else {\r
137       dd->InvMassD0bar(mD0bar);\r
138       if(TMath::Abs(pdgdaughters[0])==321) {\r
139        vars[iter]=mD0bar[0];\r
140       }else{\r
141        vars[iter]=mD0bar[1];\r
142       }\r
143    }\r
144   }\r
145 \r
146   if(fVarsForOpt[1]){\r
147     iter++;\r
148     vars[iter]=dd->GetDCA();\r
149   }\r
150 \r
151   if(fVarsForOpt[2]){\r
152     iter++;\r
153     vars[iter]=dd->GetDist12toPrim();\r
154   }\r
155   if(fVarsForOpt[3]){\r
156     iter++;\r
157     vars[iter]=dd->GetDist3toPrim();\r
158   }\r
159   if(fVarsForOpt[4]){\r
160     iter++;\r
161     vars[iter]=dd->GetDist4toPrim();\r
162   }\r
163   if(fVarsForOpt[5]){\r
164     iter++;\r
165     vars[iter]=dd->CosPointingAngle();\r
166   }\r
167   if(fVarsForOpt[6]){\r
168     iter++;\r
169     vars[iter]=dd->Pt();\r
170   }\r
171   if(fVarsForOpt[7]){\r
172     iter++;\r
173     vars[iter]=999999999.;\r
174     printf("ERROR: optmization for rho mass cut not implemented\n");\r
175   }\r
176   if(fVarsForOpt[8]){\r
177     iter++;\r
178     vars[iter]=999999999.;\r
179     printf("ERROR: optmization for PID cut not implemented\n");\r
180   }\r
181   \r
182     if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
183   return;\r
184 }\r
185 //---------------------------------------------------------------------------\r
186 Int_t AliRDHFCutsD0toKpipipi::IsSelected(TObject* obj,Int_t selectionLevel) {\r
187   //\r
188   // Apply selection\r
189   //\r
190 \r
191   if(!fCutsRD){\r
192     cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
193     return 0;\r
194   }\r
195   //PrintAll();\r
196   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
197 \r
198   if(!d){\r
199     cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
200     return 0;\r
201   }\r
202 \r
203   Double_t ptD=d->Pt();\r
204   if(ptD<fMinPtCand) return 0;\r
205   if(ptD>fMaxPtCand) return 0;
206 \r
207   // selection on daughter tracks \r
208   if(selectionLevel==AliRDHFCuts::kAll || \r
209      selectionLevel==AliRDHFCuts::kTracks) {\r
210     if(!AreDaughtersSelected(d)) return 0;\r
211   }\r
212 \r
213 \r
214   Int_t returnvalue=1;\r
215 \r
216   // selection on candidate\r
217   if(selectionLevel==AliRDHFCuts::kAll || \r
218      selectionLevel==AliRDHFCuts::kCandidate) {\r
219 \r
220     Int_t ptbin=PtBin(d->Pt());\r
221   
222     Int_t okD0=1,okD0bar=1;       \r
223     Double_t mD0[2],mD0bar[2];\r
224     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
225 \r
226     d->InvMassD0(mD0);\r
227     if(TMath::Abs(mD0[0]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)] &&\r
228        TMath::Abs(mD0[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;\r
229     d->InvMassD0bar(mD0bar);\r
230     if(TMath::Abs(mD0bar[0]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)] &&\r
231        TMath::Abs(mD0bar[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;\r
232     if(!okD0 && !okD0bar) return 0;\r
233     \r
234     if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]
235        || d->GetDCA(3) > fCutsRD[GetGlobalIndex(1,ptbin)]
236        || d->GetDCA(2) > fCutsRD[GetGlobalIndex(1,ptbin)]
237        || d->GetDCA(5) > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;\r
238     if(d->GetDist12toPrim() < fCutsRD[GetGlobalIndex(2,ptbin)]) return 0;\r
239     if(d->GetDist3toPrim() < fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;\r
240     if(d->GetDist4toPrim() < fCutsRD[GetGlobalIndex(4,ptbin)]) return 0;\r
241     if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;\r
242     if(d->Pt() < fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;\r
243     if(!d->CutRhoMass(mD0,mD0bar,fCutsRD[GetGlobalIndex(0,ptbin)],fCutsRD[GetGlobalIndex(7,ptbin)])) return 0;\r
244
245     if (okD0) returnvalue=1; //cuts passed as D0\r
246     if (okD0bar) returnvalue=2; //cuts passed as D0bar\r
247     if (okD0 && okD0bar) returnvalue=3; //cuts passed as D0 and D0bar\r
248   }\r
249
250   // selection on PID (from AliAODPidHF)\r
251   if(selectionLevel==AliRDHFCuts::kAll || \r
252      selectionLevel==AliRDHFCuts::kPID) {\r
253
254     Int_t selD01 = D01Selected(d,AliRDHFCuts::kCandidate);
255     Int_t selD02 = D02Selected(d,AliRDHFCuts::kCandidate);  
256     Int_t selD0bar1 = D0bar1Selected(d,AliRDHFCuts::kCandidate);
257     Int_t selD0bar2 = D0bar2Selected(d,AliRDHFCuts::kCandidate);
258
259     Int_t d01PID = 0, d02PID = 0, d0bar1PID = 0, d0bar2PID = 0;
260     returnvalue = IsSelectedFromPID(d, &d01PID, &d02PID, &d0bar1PID, &d0bar2PID);  //This returnvalue is dummy! Now it's modified as it must be!
261
262 returnvalue = 0;
263
264     if((selD01 == 1 && d01PID == 1)||(selD02 == 1 && d02PID == 1)||(selD0bar1 == 1 && d0bar1PID == 1)||(selD0bar2 == 1 && d0bar2PID == 1)) returnvalue = 1;\r
265   }
266
267   return returnvalue;\r
268 }
269 \r
270 //---------------------------------------------------------------------------\r
271 Int_t AliRDHFCutsD0toKpipipi::IsSelectedFromPID(AliAODRecoDecayHF4Prong *d, Int_t *hyp1, Int_t *hyp2, Int_t *hyp3, Int_t *hyp4) {\r
272   //\r
273   // Apply selection (using AliAODPidHF methods)
274   // Mass hypothesis true if each particle is at least compatible with specie of hypothesis\r
275   // \r
276
277   Int_t output=0;
278
279   Int_t matchK[4], matchPi[4];
280   Double_t ptlimit[2] = {0.6,0.8};\r
281   AliAODTrack* trk[4];
282   trk[0] = (AliAODTrack*)d->GetDaughter(0);
283   trk[1] = (AliAODTrack*)d->GetDaughter(1);
284   trk[2] = (AliAODTrack*)d->GetDaughter(2);
285   trk[3] = (AliAODTrack*)d->GetDaughter(3);
286 //  if(!trk[0] || !trk[1] || !trk[2] || !trk[3]) {          //REMOVED (needs also AliAODEvent to be passed, here and in IsSelected method)
287 //    trk[0]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
288 //    trk[1]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
289 //    trk[2]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(2)]);
290 //    trk[3]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(3)]);}
291
292   AliAODPidHF* PidObj = new AliAODPidHF();
293   PidObj->SetAsym(kTRUE);
294   PidObj->SetPLimit(ptlimit);
295   PidObj->SetSigma(0,2.);  //TPC sigma, in three pT ranges
296   PidObj->SetSigma(1,1.);
297   PidObj->SetSigma(2,0.);  
298   PidObj->SetSigma(3,2.);  //TOF sigma, whole pT range
299   PidObj->SetTPC(kTRUE);
300   PidObj->SetTOF(kTRUE);
301    
302   for(Int_t ii=0; ii<4; ii++) {
303     PidObj->SetSigma(0,2.);
304     matchK[ii] = PidObj->MatchTPCTOF(trk[ii],1,3,kTRUE); //arguments: track, mode, particle#, compatibility allowed
305     PidObj->SetSigma(0,2.);
306     matchPi[ii] = PidObj->MatchTPCTOF(trk[ii],1,2,kTRUE); 
307   }
308
309   //Rho invariant mass under various hypotheses (to be matched with PID infos in order to accet the candidate)
310   Int_t d01rho03 = 0, d01rho23 = 0, d02rho01 = 0, d02rho12 = 0, d0bar1rho12 = 0, d0bar1rho23 = 0, d0bar2rho01 = 0, d0bar2rho03 = 0;
311   if(TMath::Abs(0.775 - d->InvMassRho(0,3))<0.1)  {d01rho03 = 1; d0bar2rho03 = 1;}
312   if(TMath::Abs(0.775 - d->InvMassRho(2,3))<0.1)  {d01rho23 = 1; d0bar1rho23 = 1;}
313   if(TMath::Abs(0.775 - d->InvMassRho(0,1))<0.1)  {d02rho01 = 1; d0bar2rho01 = 1;}
314   if(TMath::Abs(0.775 - d->InvMassRho(1,2))<0.1)  {d02rho12 = 1; d0bar1rho12 = 1;}
315   Int_t d01rho = 0, d02rho = 0, d0bar1rho = 0, d0bar2rho = 0;
316   if(d01rho03==1||d01rho23==1) d01rho = 1;
317   if(d02rho01==1||d02rho12==1) d02rho = 1;
318   if(d0bar1rho12==1||d0bar1rho23==1) d0bar1rho = 1;
319   if(d0bar2rho01==1||d0bar2rho03==1) d0bar2rho = 1;
320
321   //This way because there could be multiple hypotheses accepted\r
322   if(d01rho==1 && (matchK[1]>=0 && matchPi[0]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp1 = 1; output = 1;} //d01 hyp
323   if(d02rho==1 && (matchK[3]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0)) {*hyp2 = 1; output = 1;} //d02 hyp
324   if(d0bar1rho==1 && (matchK[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp3 = 1; output = 1;} //d0bar1 hyp
325   if(d0bar2rho==1 && (matchK[2]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[3]>=0)) {*hyp4 = 1; output = 1;} //d0bar2 hyp
326
327   return output;\r
328 }\r
329 //---------------------------------------------------------------------------\r
330 Int_t AliRDHFCutsD0toKpipipi::D01Selected(TObject* obj,Int_t selectionLevel) {\r
331   //\r
332   // Apply selection\r
333   //\r
334 \r
335   if(!fCutsRD){\r
336     cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
337     return 0;\r
338   }\r
339   //PrintAll();\r
340   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
341 \r
342   if(!d){\r
343     cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
344     return 0;\r
345   }\r
346 \r
347   Int_t returnvalue=0;\r
348 \r
349   // selection on candidate\r
350   if(selectionLevel==AliRDHFCuts::kAll || \r
351      selectionLevel==AliRDHFCuts::kCandidate) {\r
352 \r
353     Int_t ptbin=PtBin(d->Pt());\r
354     \r
355     Double_t mD0[2];\r
356     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
357 \r
358     d->InvMassD0(mD0);\r
359     if(TMath::Abs(mD0[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
360   }\r
361 \r
362   return returnvalue;\r
363 }
364 //---------------------------------------------------------------------------\r
365 Int_t AliRDHFCutsD0toKpipipi::D02Selected(TObject* obj,Int_t selectionLevel) {\r
366   //\r
367   // Apply selection\r
368   //\r
369 \r
370   if(!fCutsRD){\r
371     cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
372     return 0;\r
373   }\r
374   //PrintAll();\r
375   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
376 \r
377   if(!d){\r
378     cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
379     return 0;\r
380   }\r
381 \r
382   Int_t returnvalue=0;\r
383 \r
384   // selection on candidate\r
385   if(selectionLevel==AliRDHFCuts::kAll || \r
386      selectionLevel==AliRDHFCuts::kCandidate) {\r
387 \r
388     Int_t ptbin=PtBin(d->Pt());\r
389       \r
390     Double_t mD0[2];\r
391     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
392 \r
393     d->InvMassD0(mD0);\r
394     if(TMath::Abs(mD0[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
395   }\r
396 \r
397   return returnvalue;\r
398 }\r//---------------------------------------------------------------------------\r
399 Int_t AliRDHFCutsD0toKpipipi::D0bar1Selected(TObject* obj,Int_t selectionLevel) {\r
400   //\r
401   // Apply selection\r
402   //\r
403 \r
404   if(!fCutsRD){\r
405     cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
406     return 0;\r
407   }\r
408   //PrintAll();\r
409   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
410 \r
411   if(!d){\r
412     cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
413     return 0;\r
414   }\r
415 \r
416   Int_t returnvalue=0;\r
417 \r
418   // selection on candidate\r
419   if(selectionLevel==AliRDHFCuts::kAll || \r
420      selectionLevel==AliRDHFCuts::kCandidate) {\r
421 \r
422     Int_t ptbin=PtBin(d->Pt());\r
423  \r
424     Double_t mD0bar[2];\r
425     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
426
427     d->InvMassD0bar(mD0bar);\r
428     if(TMath::Abs(mD0bar[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
429   }\r
430 \r
431   return returnvalue;\r
432 }
433 //---------------------------------------------------------------------------\r
434 Int_t AliRDHFCutsD0toKpipipi::D0bar2Selected(TObject* obj,Int_t selectionLevel) {\r
435   //\r
436   // Apply selection\r
437   //\r
438 \r
439   if(!fCutsRD){\r
440     cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
441     return 0;\r
442   }\r
443   //PrintAll();\r
444   AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
445 \r
446   if(!d){\r
447     cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
448     return 0;\r
449   }\r
450 \r
451   Int_t returnvalue=0;\r
452 \r
453   // selection on candidate\r
454   if(selectionLevel==AliRDHFCuts::kAll || \r
455      selectionLevel==AliRDHFCuts::kCandidate) {\r
456 \r
457     Int_t ptbin=PtBin(d->Pt());\r
458     \r
459     Double_t mD0bar[2];\r
460     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
461
462     d->InvMassD0bar(mD0bar);\r
463     if(TMath::Abs(mD0bar[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
464   }\r
465 \r
466   return returnvalue;\r
467 }
468 //---------------------------------------------------------------------------
469 Bool_t AliRDHFCutsD0toKpipipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
470 {
471   //
472   // Checking if D0 is in fiducial acceptance region 
473   //
474
475   if(pt > 5.) {
476     // applying cut for pt > 5 GeV
477     AliDebug(4,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt)); 
478     if (TMath::Abs(y) > 0.8){
479       return kFALSE;
480     }
481   } else {
482     // appliying smooth cut for pt < 5 GeV
483     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
484     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;         
485     AliDebug(4,Form("pt of D0 = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
486     if (y < minFiducialY || y > maxFiducialY){
487       return kFALSE;
488     }
489   }
490
491   return kTRUE;
492 }
493