fixing Coverity defects 16118 and 15959
[u/mrichter/AliRoot.git] / PWG2 / DIFFRACTIVE / AliAnalysisTaskDDMeson.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 // Reconstruct mesons 
17 // in central diffractive event 
18 // in pp collisions
19 //
20 // Authors:
21 //  Xianguo Lu <lu@physi.uni-heidelberg.de>
22 //
23
24 #include <TH1I.h>
25 #include <TH2I.h>
26 #include <TH2D.h>
27 #include <TList.h>
28 #include <TLorentzVector.h>
29 #include <TRandom3.h>
30 #include <THnSparse.h>
31
32 #include "AliCDBManager.h"
33 #include "AliCDBEntry.h"
34 #include "AliESDInputHandler.h"
35 #include "AliESDtrackCuts.h"
36 #include "AliGeomManager.h"
37 #include "AliITSAlignMille2Module.h"
38 #include "AliITSsegmentationSPD.h"
39 #include "AliKFParticle.h"
40 #include "AliMultiplicity.h"
41 #include "AliPID.h"
42 #include "AliSPDUtils.h"
43 #include "AliTriggerAnalysis.h"
44
45 #include "AliAnalysisTaskDDMeson.h"
46
47 void AliAnalysisTaskDDMeson::IniTask()
48 {
49   //
50   //initialize values used in many places
51   //
52   fnmass = 512;
53   fmass1 = 0.01 * fnmass;
54
55   fnptv = 10;
56   fptv1 = 0.2 * fnptv;
57
58   fnpta = 51*2;
59   fpta1 = 0.05 * fnpta;
60
61   fneta = 64;
62   feta = 0.04 * fneta/2.;
63
64   fnnsel = 2;
65   fnsel1 = 2+fnnsel;
66
67   fncts = 10;
68   fnctlab = 20;
69
70   fCHECKVBA = 2;
71   fCHECKVBC = 4;
72
73   //------------------
74   fOpt.ToUpper();
75 }
76
77 AliAnalysisTaskDDMeson::AliAnalysisTaskDDMeson(const TString opt):
78   AliAnalysisTaskSE("DDMeson")
79   , fOpt(opt)
80   , fESD(0x0)
81
82   , fnmass(-999)
83   , fmass1(-999)
84   , fnptv(-999)
85   , fptv1(-999)
86   , fnpta(-999)
87   , fpta1(-999)
88   , fneta(-999)
89   , feta(-999)
90   , fnnsel(-999)
91   , fnsel1(-999)
92   , fncts(-999)
93   , fnctlab(-999)
94   , fCHECKVBA(-999)
95   , fCHECKVBC(-999)
96
97   , fHBIT(0x0)
98   , fBitcg(0)
99
100   , fRun(-999)
101   , fat(0x0)
102   , fct(0x0)
103   , fbt(0x0)
104   , fnt(0x0)
105   , ftt(0x0)
106
107   , fv0ntrk(0x0)
108   , frsntrk(0x0)
109
110   , fhps(0x0)
111   , fhfo(0x0)
112   , fhspd(0x0)
113   , fhv0fmd(0x0)
114   , fhpriv(0x0)
115   , fhntrk(0x0)
116
117   , fHist(0x0)
118   , fThnMass(0x0)
119   , fThnDPt(0x0)
120   , fThnDEta(0x0)
121   , fThnKF(0x0)
122 {
123   //
124   // Dummy constructor
125   //
126   //slot in TaskSE must start from 1
127   DefineOutput(1, TList::Class());
128   DefineOutput(2, THnSparse::Class());
129   DefineOutput(3, THnSparse::Class());
130   DefineOutput(4, THnSparse::Class());
131   DefineOutput(5, THnSparse::Class());
132  
133   IniTask();
134 }
135
136 AliAnalysisTaskDDMeson::AliAnalysisTaskDDMeson(const AliAnalysisTaskDDMeson  &p):
137   AliAnalysisTaskSE("DDMeson")
138   , fOpt(p.fOpt)
139   , fESD(p.fESD)
140
141   , fnmass(p.fnmass)
142   , fmass1(p.fmass1)
143   , fnptv(p.fnptv)
144   , fptv1(p.fptv1)
145   , fnpta(p.fnpta)
146   , fpta1(p.fpta1)
147   , fneta(p.fneta)
148   , feta(p.feta)
149   , fnnsel(p.fnnsel)
150   , fnsel1(p.fnsel1)
151   , fncts(p.fncts)
152   , fnctlab(p.fnctlab)
153   , fCHECKVBA(p.fCHECKVBA)
154   , fCHECKVBC(p.fCHECKVBC)
155
156   , fHBIT(p.fHBIT)
157   , fBitcg(p.fBitcg)
158
159   , fRun(p.fRun)
160   , fat(p.fat)
161   , fct(p.fct)
162   , fbt(p.fbt)
163   , fnt(p.fnt)
164   , ftt(p.ftt)
165
166   , fv0ntrk(p.fv0ntrk)
167   , frsntrk(p.frsntrk)
168
169   , fhps(p.fhps)
170   , fhfo(p.fhfo)
171   , fhspd(p.fhspd)
172   , fhv0fmd(p.fhv0fmd)
173   , fhpriv(p.fhpriv)
174   , fhntrk(p.fhntrk)
175
176   , fHist(p.fHist)
177   , fThnMass(p.fThnMass)
178   , fThnDPt(p.fThnDPt)
179   , fThnDEta(p.fThnDEta)
180   , fThnKF(p.fThnKF)
181 {
182   //
183   // Copy constructor
184   //
185
186   IniTask();
187 }
188
189 AliAnalysisTaskDDMeson & AliAnalysisTaskDDMeson::operator=(const AliAnalysisTaskDDMeson  &p)
190 {
191   //
192   // overload =
193   //
194   fESD = p.fESD;
195
196   fnmass = p.fnmass;
197   fmass1 = p.fmass1;
198   fnptv = p.fnptv;
199   fptv1 = p.fptv1;
200   fnpta = p.fnpta;
201   fpta1 = p.fpta1;
202   fneta = p.fneta;
203   feta = p.feta;
204   fnnsel = p.fnnsel;
205   fnsel1 = p.fnsel1;
206   fncts = p.fncts;
207   fnctlab = p.fnctlab;
208   fCHECKVBA = p.fCHECKVBA;
209   fCHECKVBC = p.fCHECKVBC;
210   
211   fHBIT = p.fHBIT;
212   fBitcg = p.fBitcg;
213
214   fRun = p.fRun;
215
216   fat = p.fat;
217   fct = p.fct;
218   fbt = p.fbt;
219   fnt = p.fnt;
220   ftt = p.ftt;
221
222   fv0ntrk = p.fv0ntrk;
223   frsntrk = p.frsntrk;
224
225   fhps = p.fhps;
226   fhfo = p.fhfo;
227   fhspd = p.fhspd;
228   fhv0fmd = p.fhv0fmd;
229   fhpriv = p.fhpriv;
230   fhntrk = p.fhntrk;
231
232   fHist = p.fHist;
233   fThnMass = p.fThnMass;
234   fThnDPt = p.fThnDPt;
235   fThnDEta = p.fThnDEta;
236   fThnKF = p.fThnKF;
237
238   IniTask();
239
240   return *this;
241 }
242
243 AliAnalysisTaskDDMeson::~AliAnalysisTaskDDMeson()
244 {
245   //
246   //Destructor
247   //
248   if(fESD) delete fESD;
249
250   delete fHBIT;
251   delete fat;
252   delete fct;
253   delete fbt;
254   delete fnt;
255   delete ftt;
256
257   delete fv0ntrk;
258   delete frsntrk;
259
260   delete fhps;
261   delete fhfo;
262   delete fhspd;
263   delete fhv0fmd;
264   delete fhpriv;
265   delete fhntrk;
266
267   delete fThnMass;
268   delete fThnDPt;
269   delete fThnDEta;
270   delete fThnKF;
271
272   fHist->Clear();
273   delete fHist;  
274  
275 }
276
277 void AliAnalysisTaskDDMeson::CheckRange(Double_t &ptv, Double_t &pta, Double_t &etaa
278                                         , Double_t &mpi
279                                         ) const 
280 {
281   //
282   //save over/under flow bins
283   //
284
285   const Double_t eps = 1e-6;
286   if(ptv>=fptv1) ptv = fptv1 - eps;
287   if(pta>=fpta1) pta = fpta1 - eps;
288
289   if(etaa >= feta) etaa = feta - eps;
290
291   if(etaa <= -feta) etaa = -feta + eps;
292
293   if(mpi >= fmass1) mpi = fmass1 - eps;
294
295 }
296
297 void AliAnalysisTaskDDMeson::UserCreateOutputObjects()
298 {
299   //
300   //CreateOutputObjects
301   //
302   //=======================================
303   const Double_t kvz0=1, kvz1=5;
304   const Int_t nkvz=(Int_t)(kvz1-kvz0);
305
306   const Double_t charge0=1, charge1=4;
307   const Int_t ncharge=(Int_t)(charge1-charge0);
308   
309   const Double_t mass0 = 0;
310   //=======================================
311   //total momentum and pt in lab frame:
312   const Double_t p0=0;
313   //=======================================
314   const Double_t nsel0 = 2;
315
316   //=======================================
317   const Double_t cts0=0, cts1=1;
318
319   //=======================================
320   const Double_t ctlab0 = -1, ctlab1 = 1;
321
322   //=======================================
323   //------- Mass
324   const Int_t    nparPtMass = 7;
325   const Int_t    nbinPtMass[]     ={fnnsel, nkvz, ncharge,  fnmass, fnptv, fncts, fnctlab};
326   const Double_t xminPtMass[]     ={nsel0,  kvz0, charge0,   mass0,    p0,  cts0,  ctlab0};
327   const Double_t xmaxPtMass[]     ={fnsel1, kvz1, charge1,  fmass1, fptv1,  cts1,  ctlab1};
328   fThnMass = new THnSparseD("DDMeson_Mass","nsel, kvz, charge, Mpipi, ptv, cts, ctlab", nparPtMass, nbinPtMass, xminPtMass, xmaxPtMass);
329  
330   //------- DPt
331   const Int_t    nparDPt = 4;
332   const Int_t    nbinDPt[] ={fnnsel, nkvz, ncharge,  fnpta};
333   const Double_t xminDPt[] ={nsel0, kvz0, charge0,  p0};
334   const Double_t xmaxDPt[] ={fnsel1, kvz1, charge1,  fpta1};
335   fThnDPt = new THnSparseD("DDMeson_DPt","nsel, kvz, charge, pt1", nparDPt, nbinDPt, xminDPt, xmaxDPt);
336   //------- DEta
337   const Int_t    nparDEta = 4;
338   const Int_t    nbinDEta[] ={fnnsel, nkvz, ncharge,  fneta};
339   const Double_t xminDEta[] ={nsel0, kvz0, charge0,  -feta};
340   const Double_t xmaxDEta[] ={fnsel1, kvz1, charge1,  feta};
341   fThnDEta = new THnSparseD("DDMeson_DEta"," nsel, kvz, charge, eta1", nparDEta, nbinDEta, xminDEta, xmaxDEta);
342
343   //------- KF
344   const Double_t mks = 0.494;
345   const Double_t dks = 0.024;
346   const Int_t nparKF = 4;
347   const Int_t nbinKF[] = {nkvz, 200, 200, 200};
348   const Double_t xminKF[] ={kvz0, mks-dks, mks-dks, 0};
349   const Double_t xmaxKF[] ={kvz1, mks+dks, mks+dks, 10};
350   fThnKF = new THnSparseD("DDMeson_KF","kvz, rawks, kfks, chi2", nparKF, nbinKF, xminKF, xmaxKF);
351   //=======================================
352   //=======================================
353   fHBIT = new TH1I("HBIT","",32,1,33);
354
355   const Int_t runb0 = 114650; //runa1=114649, runb1=117630, runc0=117631, runc1=121526, rund0=121527, rund1=126460, rune0=126461, rune1 = 130930; //runf0=130931, all checked on elog
356   const Int_t runf1 = 133900;//02OCT2010
357   const Int_t run0 = runb0-1, run1 = runf1+1;
358   const Int_t nrun = (Int_t)(run1-run0);
359
360   fat = new TH1I("at","",nrun,run0,run1);
361   fct = new TH1I("ct","",nrun,run0,run1);
362   fbt = new TH1I("bt","",nrun,run0,run1);
363   fnt = new TH1I("nt","",nrun,run0,run1);
364   ftt = new TH1I("tt","",nrun,run0,run1);
365   
366   fv0ntrk = new TH2D("c_v0ntrk","",80,0,80, 4, 1, 5);//x: ntrk; y: GetV0
367   frsntrk = new TH2D("c_rsntrk","",1000,0,1000, 200, 0,200);
368  
369   fhps =    new TH1I("a000_ps",    "", 20, 0,20);
370   fhfo =    new TH2I("a010_fo",    "", 50, 0, 50, 50, 0, 50);
371   fhspd =   new TH1I("a011_spd",   "", 50, 0, 50);
372   fhv0fmd = new TH2I("a020_v0fmd", "", 4, 0, 4, 4, 0, 4);
373   fhpriv =  new TH1I("a100_priv",  "", 2, 0,2);
374   fhntrk =  new TH1I("a110_ntrk",  "", (Int_t)fnsel1, 0, (Int_t)fnsel1);
375
376   //------
377
378   fHist = new TList;
379
380   //------>>> Very important!!!!!!!
381   fHist->SetOwner();
382
383   fHist->Add(fat);
384   fHist->Add(fct);
385   fHist->Add(fbt);
386   fHist->Add(fnt);
387   fHist->Add(ftt);
388
389   fHist->Add(fv0ntrk);
390   fHist->Add(frsntrk);
391   fHist->Add(fhps);
392   fHist->Add(fhfo);
393   fHist->Add(fhspd);
394   fHist->Add(fhv0fmd);
395   fHist->Add(fhpriv);
396   fHist->Add(fhntrk);
397 }
398
399 void AliAnalysisTaskDDMeson::UserExec(Option_t *)
400 {
401   //
402   //Execute
403   // 
404   //----------------------------------- general protection
405   if(!CheckESD())
406     return;
407
408   //----------------------------------- trigger bits
409   if(!CheckBit())
410     return;  
411
412   //----------------------------------- track cuts
413   const Int_t ntrk0 = fESD->GetNumberOfTracks();
414   const AliESDtrack * trks[ntrk0];
415   for(Int_t ii=0; ii<ntrk0; ii++){
416     trks[ii] = 0x0;
417   }
418   const Int_t nsel = CutESD(trks);
419   if(nsel<2)
420     return;
421
422   for(Int_t ii=0; ii<nsel; ii++){
423     for(Int_t jj=ii+1; jj<nsel; jj++){
424       const AliESDtrack * trkpair[2]={trks[ii], trks[jj]};
425
426       //----------------------------------- tracks initailization
427       TRandom3 tmprd(0);
428       if(tmprd.Rndm()>0.5){
429         SwapTrack(trkpair);
430       }
431   
432       //------------------------------------------------------------------------
433       const Int_t pipdg = 211;
434       const AliKFParticle daughter0(*(trkpair[0]), pipdg);
435       const AliKFParticle daughter1(*(trkpair[1]), pipdg);
436
437       const AliKFParticle parent(daughter0, daughter1);
438       const Double_t kfchi2 = parent.GetChi2();
439       const Double_t kfmass = parent.GetMass();
440
441       //------------------------------------------------------------------------
442       const Int_t kvz = GetV0();
443
444       const AliESDtrack * trk1=trkpair[0];
445       const AliESDtrack * trk2=trkpair[1];
446   
447       const Double_t ch1 = trk1->GetSign();
448       const Double_t ch2 = trk2->GetSign();
449       const Int_t combch = GetCombCh(ch1, ch2);
450       
451       Double_t pta  = trk1->Pt();
452       Double_t etaa = trk1->Eta();
453       
454       Double_t tmpp1[3], tmpp2[3];
455       trk1->GetPxPyPz(tmpp1);
456       trk2->GetPxPyPz(tmpp2);
457   
458       //------------------------------------------------------------------------
459       //------------------------------------------------------------------------
460       const Double_t ctlab = GetCtlab(tmpp1, tmpp2);
461
462       Double_t cts = -999;
463       const Double_t masspi = AliPID::ParticleMass(AliPID::kPion);
464       const TLorentzVector sumv = GetKinematics(tmpp1, tmpp2, masspi, masspi, cts);
465       Double_t pipi = sumv.Mag();
466       Double_t ptv = sumv.Pt();
467
468       CheckRange(ptv, pta, etaa, pipi);
469       
470       const Double_t varMass[]  ={nsel, kvz, combch
471                                        , pipi, ptv, cts
472                                        , ctlab
473       };
474       const Double_t varDPt[]        ={nsel, kvz, combch, pta};
475       const Double_t varDEta[]       ={nsel, kvz, combch, etaa};
476       fThnMass->Fill(varMass);
477       fThnDPt->Fill(varDPt);
478       fThnDEta->Fill(varDEta);
479
480       const Double_t varKF[] = {kvz, pipi, kfmass, kfchi2};
481       if(nsel==2 && combch==3)
482         fThnKF->Fill(varKF);
483     }
484   }
485   //------------------------------------------------------------------------
486   //------------------------------------------------------------------------
487
488   PostData(1, fHist);
489   PostData(2, fThnMass);
490   PostData(3, fThnDPt);
491   PostData(4, fThnDEta);
492   PostData(5, fThnKF);
493 }
494 //=======================================================================
495 //=======================================================================
496     
497 void AliAnalysisTaskDDMeson::Terminate(Option_t *)
498 {
499   //
500   //fillbit before leaving
501   //
502   FillBit();
503 }
504
505 void AliAnalysisTaskDDMeson::FillBit()
506 {
507   //
508   //save the bit count in f*t 
509   //
510   Double_t tot[5];
511   CalcBit(fHBIT, tot);
512
513   const Int_t bin = fat->FindBin(fRun);
514   fat->AddBinContent(bin, tot[0]);
515   fct->AddBinContent(bin, tot[1]);
516   fbt->AddBinContent(bin, tot[2]);
517   fnt->AddBinContent(bin, tot[3]);
518   ftt->AddBinContent(bin, tot[4]);
519 }
520
521 void AliAnalysisTaskDDMeson::CalcBit(TH1I *hc, Double_t tot[]) const
522 {
523   //
524   //calculate the bit count in f*t 
525   //
526   if(hc->GetBinLowEdge(1)!=1){
527     printf("xlulog hc defined incorrectly!! %f\n",hc->GetBinLowEdge(1));
528     return;
529   }
530   if(hc->GetBinContent(0)){
531     printf("\nxlulog hc !!!!!!!!!!!!!!!!!!! underflow!!!!!!!!!!!!!!!!! %f\n\n", hc->GetBinContent(0));
532     return;
533   }
534   if(hc->GetBinContent(hc->GetNbinsX())){
535     printf("\nxlulog hc !!!!!!!!!!!!!!!!! OVERflow!!!!!!!!!!!!!! %f\n\n", hc->GetBinContent(hc->GetNbinsX()));
536     return;
537   }
538
539
540   Double_t atot=0, ctot=0, btot=0, ntot=0;
541   for(Int_t ib=1; ib<=hc->GetNbinsX(); ib++){
542     const Double_t cont = hc->GetBinContent(ib);
543     if(       (ib & fCHECKVBA) && !(ib & fCHECKVBC) ){
544       atot += cont;
545     }
546     else if( !(ib & fCHECKVBA) &&  (ib & fCHECKVBC) ){
547       ctot += cont;
548     }
549     else if(  (ib & fCHECKVBA) &&  (ib & fCHECKVBC) ){
550       btot += cont;
551     }
552     else{
553       ntot += cont;
554     }
555   }
556
557   const Double_t ttot   = atot + ctot + btot + ntot;
558   
559   tot[0]=atot;
560   tot[1]=ctot;
561   tot[2]=btot;
562   tot[3]=ntot;
563   tot[4]=ttot;
564
565   //---
566   char fullname[5][20];
567   strncpy(fullname[0],"V0A-ONLY",8);
568   strncpy(fullname[1],"V0C-ONLY",8);
569   strncpy(fullname[2],"V0A & V0C",9);
570   strncpy(fullname[3],"!V0A & !V0C",11);
571   strncpy(fullname[4],"TOTAL",5);
572
573   for(Int_t ii=0;ii<5;ii++){
574     printf("xlulog %s CTP-L0: %s\tTOT: %10.0f\n",fOpt.Data(), fullname[ii],tot[ii]);
575   }
576
577   //---
578
579   const Int_t nb=hc->GetNbinsX();
580   for(Int_t ib=0; ib<=nb+1; ib++){
581     hc->SetBinContent(ib, 0);
582   }
583 }
584
585 //====================================================================================
586
587 Bool_t AliAnalysisTaskDDMeson::CheckESD()
588 {
589   //
590   //general protection
591   //
592   const AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
593   if(esdH)  
594     fESD = esdH->GetEvent();
595   if(!fESD){
596     AliError("xlulog No ESD Event");
597     return 0;
598   }
599   
600   if(fabs(fESD->GetMagneticField())<1){
601     printf("xlulog strange Bfield! %f\n", fESD->GetMagneticField());
602     return 0;
603   }
604
605   const Int_t tmprun = fESD->GetRunNumber();
606
607   if(fRun!=tmprun){
608     if(fRun>0){
609       printf("xlulog run changed during exec!! %d %d\n", tmprun, fRun);
610       FillBit();
611     }
612     fRun = tmprun;
613     if(fOpt.Contains("SPD")){
614       SPDLoadGeom();
615     }
616   }
617
618   return 1;  
619 }
620
621 void AliAnalysisTaskDDMeson::SPDLoadGeom() const 
622 {
623   //method to get the gGeomanager
624   // it is called at the CreatedOutputObject stage
625   // to comply with the CAF environment
626   AliCDBManager *man = AliCDBManager::Instance();
627   TString cdbpath("local:///lustre/alice/alien/alice/data/2010/OCDB");
628   man->SetDefaultStorage(cdbpath);
629   man->SetRun(fRun);  
630   man->SetSpecificStorage("ITS/Align/Data",cdbpath);
631   man->SetSpecificStorage("GRP/Geometry/Data",cdbpath);
632   
633   AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
634   if (obj != NULL) {
635     AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
636     AliGeomManager::ApplyAlignObjsFromCDB("ITS"); 
637     printf("xlulog DDMeson %s Geom Loaded for run %d !\n", fOpt.Data(), fRun);
638   }
639   else {
640     printf("xlulog DDMeson %s was unable to load Geom for run %d !\n", fOpt.Data(), fRun);
641   }
642 }
643
644 Bool_t AliAnalysisTaskDDMeson::SPDLoc2Glo(const Int_t id, const Double_t *loc, Double_t *glo) const 
645 {
646   //
647   //SPDLoc2Glo
648   //
649
650   static TGeoHMatrix mat;
651   Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
652   if (vid<0) {
653     printf("xlulog Did not find module with such ID %d\n",id);
654     return kFALSE;
655   }
656   AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
657   mat.LocalToMaster(loc,glo);
658   return kTRUE;
659 }
660
661 Bool_t AliAnalysisTaskDDMeson::CheckChipEta(const Int_t chipKey) const 
662 {
663   //
664   //CheckChipEta
665   //
666
667   //no eta cut 
668   if(fOpt.Contains("SPD0"))
669     return kTRUE;
670
671   const Double_t etacut = 1.0;
672
673   UInt_t module=999, offchip=999;
674   AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
675   UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
676   if(hs<2) offchip = 4 - offchip; // inversion  in the inner layer... 
677   
678   const Int_t col[]={
679     hs<2? 0 : 31, 
680     hs<2? 31 : 0, 
681     hs<2? 31 : 0, 
682     hs<2? 0 : 31};
683   const Int_t aa[]={0, 0, 255, 255};
684   const AliITSsegmentationSPD seg;
685
686   for(Int_t ic=0; ic<4; ic++){
687     Float_t localchip[3]={0.,0.,0.};
688     seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]); // local coordinate of the chip center
689     //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
690     const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
691     Double_t glochip[3]={0.,0.,0.};
692     if(!SPDLoc2Glo(module,local,glochip)){
693       return kFALSE;
694     }
695
696     const TVector3 pos(glochip[0],glochip[1],glochip[2]);
697     //pos.Print();
698
699     if(fabs(pos.Eta()) > etacut)
700       return kFALSE;
701   }
702
703   return kTRUE;
704 }
705
706 void AliAnalysisTaskDDMeson::GetNFO(Int_t &ni, Int_t &no) const 
707 {
708   //
709   //GetNFO
710   //
711
712   const AliMultiplicity *mult = fESD->GetMultiplicity();
713   Int_t nFOinner=0;
714   Int_t nFOouter=0;
715   for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
716     if(mult->TestFastOrFiredChips(iChipKey) && CheckChipEta(iChipKey)){ // here you check if the FastOr bit is 1 or 0
717       if(iChipKey<400) 
718         nFOinner++;  // here you count the FastOr bits in the inner layer
719       else 
720         nFOouter++;  // here you count the FastOr bits in the outer layer
721     }
722   }
723
724   ni = nFOinner;
725   no = nFOouter;
726 }
727
728 Bool_t AliAnalysisTaskDDMeson::CheckBit()
729 {
730   //
731   //get offline trigger
732   //
733   fhps->Fill(1);
734
735   AliTriggerAnalysis triggerAnalysis;
736
737   //hardware 1: hw 0: offline
738   const Bool_t khw = kFALSE;
739   const Bool_t v0A       = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BB);
740   const Bool_t v0C       = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BB);
741
742   /*
743     const Bool_t v0ABG     = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BG);
744     const Bool_t v0CBG     = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BG);
745     const Bool_t v0BG      = v0ABG || v0CBG;
746     const Int_t fastOROffline = triggerAnalysis.SPDFiredChips(fESD, 0); // SPD number of chips from clusters (!)        
747
748     //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliPhysicsSelection.cxx?view=markup&root=AliRoot
749     //Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
750    */
751
752   const Bool_t fmdA = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kASide);
753   const Bool_t fmdC = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kCSide);
754
755   Bool_t bitA = kFALSE, bitC = kFALSE;
756   if(fOpt.Contains("V0")){
757     bitA = bitA || v0A;
758     bitC = bitC || v0C;
759   }
760   if(fOpt.Contains("FMD")){
761     bitA = bitA || fmdA;
762     bitC = bitC || fmdC;
763   }
764   
765   //---------------------------------------------------------------------------------------------
766   //---------------------------------------------------------------------------------------------
767   //considering the efficiency per chip, nfo and fastORhw is required only to have non-0 count!!  
768
769   //simulate SPD _hardware_ trigger in eta range
770   if(fOpt.Contains("SPD")){
771     Int_t ni=-999, no = -999;
772     GetNFO(ni, no);
773     if(!bitA && !bitC){
774       fhfo->Fill(ni, no);
775     }
776     if(ni<2 || no<2)
777       return 0;
778   }
779
780   //simulate MB condition, since for double gap SPD is definitely fired, fastORhw>0 to reduce GA, GC, NG
781   //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
782   //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
783   // returns the number of fired chips in the SPD
784   // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
785   // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
786   const Int_t fastORhw = triggerAnalysis.SPDFiredChips(fESD, 1);
787   fhspd->Fill(fastORhw);
788   if(fastORhw<1)
789     return 0;
790
791   //---------------------------------------------------------------------------------------------
792   //---------------------------------------------------------------------------------------------
793     
794   const Int_t iv0 = v0A + 2*v0C;
795   const Int_t ifmd = fmdA + 2*fmdC;
796   fhv0fmd->Fill(iv0, ifmd);
797
798   //---------------------------------------------------------------------------------------------
799   //---------------------------------------------------------------------------------------------
800   
801   //need a base number to distringuish from null entry; fHBit starts from 1; important!!
802   fBitcg=1; 
803
804   if(bitA)
805     fBitcg += fCHECKVBA;
806   if(bitC)
807     fBitcg += fCHECKVBC;
808
809   fHBIT->Fill( fBitcg );
810   
811   return 1;
812 }
813
814 Int_t AliAnalysisTaskDDMeson::CutESD(const AliESDtrack *outtrk[])
815 {
816   //
817   //track cuts
818   //
819   
820   //collision vertex cut
821   //A cut in XY is implicitly done during the reconstruction by constraining the vertex to the beam diamond.
822
823   // Primary vertex
824   Bool_t kpr0 = kTRUE;
825   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
826   if(vertex->GetNContributors()<1) {
827     // SPD vertex
828     vertex = fESD->GetPrimaryVertexSPD();
829     if(vertex->GetNContributors()<1) {
830       // NO GOOD VERTEX, SKIP EVENT 
831       kpr0 = kFALSE;
832     }
833   }
834   const Bool_t kpriv = kpr0 && ( fabs(fESD->GetPrimaryVertex()->GetZ())<10 ); 
835   fhpriv->Fill(kpriv);
836   if(!kpriv)  
837     return 0;
838   
839   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts;//::GetStandardITSTPCTrackCuts2010(kTRUE);
840
841   //--->
842   // TPC  
843   esdTrackCuts->SetMinNClustersTPC(70);
844   esdTrackCuts->SetMaxChi2PerClusterTPC(4);
845   esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
846   esdTrackCuts->SetRequireTPCRefit(kTRUE);
847   // ITS
848   esdTrackCuts->SetRequireITSRefit(kTRUE);
849   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
850                                          AliESDtrackCuts::kAny);
851   // 7*(0.0026+0.0050/pt^1.01)
852   esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
853
854   esdTrackCuts->SetMaxDCAToVertexZ(2);
855   esdTrackCuts->SetDCAToVertex2D(kFALSE);
856   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
857   //---<
858
859   esdTrackCuts->SetEtaRange(-0.9, 0.9);
860   
861   const TObjArray* seltrk = esdTrackCuts->GetAcceptedTracks(fESD);
862   const Int_t n0=seltrk->GetEntries();
863   const AliESDtrack * trks[n0];
864   Int_t nsel=0;
865   for(Int_t ii=0; ii<n0; ii++){
866     trks[ii]=0x0;
867     const AliESDtrack *esdtrack = (AliESDtrack *)(seltrk->At(ii));
868
869     /*
870     if(!CutTrack(esdtrack))
871       continue;
872     */
873
874     trks[nsel++]=esdtrack;
875   }
876
877   delete seltrk;
878   delete esdTrackCuts;
879
880   fv0ntrk->Fill(nsel,GetV0());
881
882   fhntrk->Fill(nsel);
883
884   frsntrk->Fill(fESD->GetNumberOfTracks(), nsel);
885
886   if(nsel>=fnsel1)
887     return 0;
888
889   for(Int_t ii=0; ii<nsel; ii++){
890     outtrk[ii]=trks[ii];
891   }
892
893   return nsel;
894 }
895
896 /*
897 Bool_t AliAnalysisTaskDDMeson::CutTrack(const AliESDtrack *  esdtrack) const 
898 {
899   //
900   //track cut
901   //  
902
903   return kTRUE;
904 }
905 */
906
907 //=================================================================================
908
909 void AliAnalysisTaskDDMeson::SwapTrack(const AliESDtrack * trks[]) const 
910 {
911   //
912   //swap two esdtracks
913   //
914   const AliESDtrack * tmpt = trks[0];
915   trks[0]=trks[1];
916   trks[1]=tmpt;
917 }
918
919 Int_t AliAnalysisTaskDDMeson::GetV0() const
920 {
921   //
922   //V0 bit combination
923   //
924
925   const Bool_t ka = (fBitcg & fCHECKVBA);
926   const Bool_t kc = (fBitcg & fCHECKVBC);
927   
928
929   //a1c0 a0c1 must be adjacent so that a cut for signal gas can be eaily applied
930   const Int_t a0c0 = 1;
931   const Int_t a1c0 = 2;
932   const Int_t a0c1 = 3;
933   const Int_t a1c1 = 4;
934
935   if(!ka && !kc)
936     return a0c0;
937   else{
938     if(ka && kc)
939       return a1c1;
940     else if(ka)
941       return a1c0;
942     else 
943       return a0c1;
944   }
945 }
946
947 Int_t AliAnalysisTaskDDMeson::GetCombCh(const Double_t s1, const Double_t s2) const
948 {
949   //
950   //get combination of charges
951   //
952   const Int_t combPP = 1;
953   const Int_t combMM = 2;
954   const Int_t combPM = 3;
955
956   if(s1*s2<0){
957     return combPM;
958   }
959   else if(s1>0)
960     return combPP;
961   else
962     return combMM;
963 }
964
965 //==========================================================================
966
967 TLorentzVector AliAnalysisTaskDDMeson::GetKinematics(const Double_t *pa, const Double_t *pb, const Double_t ma, const Double_t mb, Double_t & cts) const
968 {
969   //
970   //get kinematics, cts = cos(theta_{#star})
971   //
972   TLorentzVector va, vb, sumv;
973   va.SetXYZM(pa[0], pa[1], pa[2], ma);
974   vb.SetXYZM(pb[0], pb[1], pb[2], mb);
975   sumv = va+vb;
976
977   const TVector3 bv = -sumv.BoostVector();
978
979   va.Boost(bv);
980   vb.Boost(bv);
981   //printf("test a %f %f %f -- %f %f\n", va.X(), va.Y(), va.Z(), va.Theta(), va.Phi());
982   //printf("test b %f %f %f -- %f %f\n", vb.X(), vb.Y(), vb.Z(), vb.Theta(), vb.Phi());
983
984   cts = fabs(cos(va.Theta()));
985
986   return sumv;
987 }
988
989 Double_t AliAnalysisTaskDDMeson::GetCtlab(const Double_t *pa, const Double_t *pb) const 
990 {
991   //
992   //cos(theat_lab)
993   //
994
995   TVector3 va, vb;
996   va.SetXYZ(pa[0], pa[1], pa[2]);
997   vb.SetXYZ(pb[0], pb[1], pb[2]);
998
999   const TVector3 ua = va.Unit();
1000   const TVector3 ub = vb.Unit();
1001
1002   return ua.Dot(ub);
1003 }