Various changes and fixes
[u/mrichter/AliRoot.git] / PWG2 / DIFFRACTIVE / AliAnalysisTaskDDMeson.cxx
CommitLineData
878f37d6 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
47void 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
77AliAnalysisTaskDDMeson::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
136AliAnalysisTaskDDMeson::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
189AliAnalysisTaskDDMeson & 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
243AliAnalysisTaskDDMeson::~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
277void 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
297void 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
399void 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
497void AliAnalysisTaskDDMeson::Terminate(Option_t *)
498{
499 //
500 //fillbit before leaving
501 //
502 FillBit();
503}
504
505void 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
521void 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 strcpy(fullname[0],"V0A-ONLY");
568 strcpy(fullname[1],"V0C-ONLY");
569 strcpy(fullname[2],"V0A & V0C");
570 strcpy(fullname[3],"!V0A & !V0C");
571 strcpy(fullname[4],"TOTAL");
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
587Bool_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
621void 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 AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
635 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
636 printf("xlulog DDMeson %s Geom Loaded for run %d !\n", fOpt.Data(), fRun);
637}
638
639Bool_t AliAnalysisTaskDDMeson::SPDLoc2Glo(const Int_t id, const Double_t *loc, Double_t *glo) const
640{
641 //
642 //SPDLoc2Glo
643 //
644
645 static TGeoHMatrix mat;
646 Int_t vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
647 if (vid<0) {
648 printf("xlulog Did not find module with such ID %d\n",id);
649 return kFALSE;
650 }
651 AliITSAlignMille2Module::SensVolMatrix(vid,&mat);
652 mat.LocalToMaster(loc,glo);
653 return kTRUE;
654}
655
656Bool_t AliAnalysisTaskDDMeson::CheckChipEta(const Int_t chipKey) const
657{
658 //
659 //CheckChipEta
660 //
661
662 //no eta cut
663 if(fOpt.Contains("SPD0"))
664 return kTRUE;
665
666 const Double_t etacut = 1.0;
667
668 UInt_t module=999, offchip=999;
669 AliSPDUtils::GetOfflineFromOfflineChipKey(chipKey,module,offchip);
670 UInt_t hs = AliSPDUtils::GetOnlineHSFromOffline(module);
671 if(hs<2) offchip = 4 - offchip; // inversion in the inner layer...
672
673 const Int_t col[]={
674 hs<2? 0 : 31,
675 hs<2? 31 : 0,
676 hs<2? 31 : 0,
677 hs<2? 0 : 31};
678 const Int_t aa[]={0, 0, 255, 255};
679 const AliITSsegmentationSPD seg;
680
681 for(Int_t ic=0; ic<4; ic++){
682 Float_t localchip[3]={0.,0.,0.};
683 seg.DetToLocal(aa[ic],col[ic]+32*offchip,localchip[0],localchip[2]); // local coordinate of the chip center
684 //printf("local coordinates %d %d: %f %f \n",chipKey, ic, localchip[0],localchip[2]);
685 const Double_t local[3] = {localchip[0],localchip[1],localchip[2]};
686 Double_t glochip[3]={0.,0.,0.};
687 if(!SPDLoc2Glo(module,local,glochip)){
688 return kFALSE;
689 }
690
691 const TVector3 pos(glochip[0],glochip[1],glochip[2]);
692 //pos.Print();
693
694 if(fabs(pos.Eta()) > etacut)
695 return kFALSE;
696 }
697
698 return kTRUE;
699}
700
701void AliAnalysisTaskDDMeson::GetNFO(Int_t &ni, Int_t &no) const
702{
703 //
704 //GetNFO
705 //
706
707 const AliMultiplicity *mult = fESD->GetMultiplicity();
708 Int_t nFOinner=0;
709 Int_t nFOouter=0;
710 for(Int_t iChipKey=0; iChipKey < 1200; iChipKey++){
711 if(mult->TestFastOrFiredChips(iChipKey) && CheckChipEta(iChipKey)){ // here you check if the FastOr bit is 1 or 0
712 if(iChipKey<400)
713 nFOinner++; // here you count the FastOr bits in the inner layer
714 else
715 nFOouter++; // here you count the FastOr bits in the outer layer
716 }
717 }
718
719 ni = nFOinner;
720 no = nFOouter;
721}
722
723Bool_t AliAnalysisTaskDDMeson::CheckBit()
724{
725 //
726 //get offline trigger
727 //
728 fhps->Fill(1);
729
730 AliTriggerAnalysis triggerAnalysis;
731
732 //hardware 1: hw 0: offline
733 const Bool_t khw = kFALSE;
734 const Bool_t v0A = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BB);
735 const Bool_t v0C = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BB);
736
737 /*
738 const Bool_t v0ABG = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kASide, khw) == AliTriggerAnalysis::kV0BG);
739 const Bool_t v0CBG = (triggerAnalysis.V0Trigger(fESD, AliTriggerAnalysis::kCSide, khw) == AliTriggerAnalysis::kV0BG);
740 const Bool_t v0BG = v0ABG || v0CBG;
741 const Int_t fastOROffline = triggerAnalysis.SPDFiredChips(fESD, 0); // SPD number of chips from clusters (!)
742
743 //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliPhysicsSelection.cxx?view=markup&root=AliRoot
744 //Bool_t mb1 = (fastOROffline > 0 || v0A || v0C) && (!v0BG);
745 */
746
747 const Bool_t fmdA = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kASide);
748 const Bool_t fmdC = triggerAnalysis.FMDTrigger(fESD, AliTriggerAnalysis::kCSide);
749
750 Bool_t bitA = kFALSE, bitC = kFALSE;
751 if(fOpt.Contains("V0")){
752 bitA = bitA || v0A;
753 bitC = bitC || v0C;
754 }
755 if(fOpt.Contains("FMD")){
756 bitA = bitA || fmdA;
757 bitC = bitC || fmdC;
758 }
759
760 //---------------------------------------------------------------------------------------------
761 //---------------------------------------------------------------------------------------------
762 //considering the efficiency per chip, nfo and fastORhw is required only to have non-0 count!!
763
764 //simulate SPD _hardware_ trigger in eta range
765 if(fOpt.Contains("SPD")){
766 Int_t ni=-999, no = -999;
767 GetNFO(ni, no);
768 if(!bitA && !bitC){
769 fhfo->Fill(ni, no);
770 }
771 if(ni<2 || no<2)
772 return 0;
773 }
774
775 //simulate MB condition, since for double gap SPD is definitely fired, fastORhw>0 to reduce GA, GC, NG
776 //http://alisoft.cern.ch/viewvc/trunk/ANALYSIS/AliTriggerAnalysis.cxx?view=markup&root=AliRoot
777 //Int_t AliTriggerAnalysis::SPDFiredChips(const AliESDEvent* aEsd, Int_t origin, Bool_t fillHists, Int_t layer)
778 // returns the number of fired chips in the SPD
779 // origin = 0 --> aEsd->GetMultiplicity()->GetNumberOfFiredChips() (filled from clusters)
780 // origin = 1 --> aEsd->GetMultiplicity()->TestFastOrFiredChips() (from hardware bits)
781 const Int_t fastORhw = triggerAnalysis.SPDFiredChips(fESD, 1);
782 fhspd->Fill(fastORhw);
783 if(fastORhw<1)
784 return 0;
785
786 //---------------------------------------------------------------------------------------------
787 //---------------------------------------------------------------------------------------------
788
789 const Int_t iv0 = v0A + 2*v0C;
790 const Int_t ifmd = fmdA + 2*fmdC;
791 fhv0fmd->Fill(iv0, ifmd);
792
793 //---------------------------------------------------------------------------------------------
794 //---------------------------------------------------------------------------------------------
795
796 //need a base number to distringuish from null entry; fHBit starts from 1; important!!
797 fBitcg=1;
798
799 if(bitA)
800 fBitcg += fCHECKVBA;
801 if(bitC)
802 fBitcg += fCHECKVBC;
803
804 fHBIT->Fill( fBitcg );
805
806 return 1;
807}
808
809Int_t AliAnalysisTaskDDMeson::CutESD(const AliESDtrack *outtrk[])
810{
811 //
812 //track cuts
813 //
814
815 //collision vertex cut
816 //A cut in XY is implicitly done during the reconstruction by constraining the vertex to the beam diamond.
817
818 // Primary vertex
819 Bool_t kpr0 = kTRUE;
820 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
821 if(vertex->GetNContributors()<1) {
822 // SPD vertex
823 vertex = fESD->GetPrimaryVertexSPD();
824 if(vertex->GetNContributors()<1) {
825 // NO GOOD VERTEX, SKIP EVENT
826 kpr0 = kFALSE;
827 }
828 }
829 const Bool_t kpriv = kpr0 && ( fabs(fESD->GetPrimaryVertex()->GetZ())<10 );
830 fhpriv->Fill(kpriv);
831 if(!kpriv)
832 return 0;
833
834 AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts;//::GetStandardITSTPCTrackCuts2010(kTRUE);
835
836 //--->
837 // TPC
838 esdTrackCuts->SetMinNClustersTPC(70);
839 esdTrackCuts->SetMaxChi2PerClusterTPC(4);
840 esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
841 esdTrackCuts->SetRequireTPCRefit(kTRUE);
842 // ITS
843 esdTrackCuts->SetRequireITSRefit(kTRUE);
844 esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
845 AliESDtrackCuts::kAny);
846 // 7*(0.0026+0.0050/pt^1.01)
847 esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
848
849 esdTrackCuts->SetMaxDCAToVertexZ(2);
850 esdTrackCuts->SetDCAToVertex2D(kFALSE);
851 esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
852 //---<
853
854 esdTrackCuts->SetEtaRange(-0.9, 0.9);
855
856 const TObjArray* seltrk = esdTrackCuts->GetAcceptedTracks(fESD);
857 const Int_t n0=seltrk->GetEntries();
858 const AliESDtrack * trks[n0];
859 Int_t nsel=0;
860 for(Int_t ii=0; ii<n0; ii++){
861 trks[ii]=0x0;
862 const AliESDtrack *esdtrack = (AliESDtrack *)(seltrk->At(ii));
863
864 /*
865 if(!CutTrack(esdtrack))
866 continue;
867 */
868
869 trks[nsel++]=esdtrack;
870 }
871
872 delete seltrk;
873 delete esdTrackCuts;
874
875 fv0ntrk->Fill(nsel,GetV0());
876
877 fhntrk->Fill(nsel);
878
879 frsntrk->Fill(fESD->GetNumberOfTracks(), nsel);
880
881 if(nsel>=fnsel1)
882 return 0;
883
884 for(Int_t ii=0; ii<nsel; ii++){
885 outtrk[ii]=trks[ii];
886 }
887
888 return nsel;
889}
890
891/*
892Bool_t AliAnalysisTaskDDMeson::CutTrack(const AliESDtrack * esdtrack) const
893{
894 //
895 //track cut
896 //
897
898 return kTRUE;
899}
900*/
901
902//=================================================================================
903
904void AliAnalysisTaskDDMeson::SwapTrack(const AliESDtrack * trks[]) const
905{
906 //
907 //swap two esdtracks
908 //
909 const AliESDtrack * tmpt = trks[0];
910 trks[0]=trks[1];
911 trks[1]=tmpt;
912}
913
914Int_t AliAnalysisTaskDDMeson::GetV0() const
915{
916 //
917 //V0 bit combination
918 //
919
920 const Bool_t ka = (fBitcg & fCHECKVBA);
921 const Bool_t kc = (fBitcg & fCHECKVBC);
922
923
924 //a1c0 a0c1 must be adjacent so that a cut for signal gas can be eaily applied
925 const Int_t a0c0 = 1;
926 const Int_t a1c0 = 2;
927 const Int_t a0c1 = 3;
928 const Int_t a1c1 = 4;
929
930 if(!ka && !kc)
931 return a0c0;
932 else{
933 if(ka && kc)
934 return a1c1;
935 else if(ka)
936 return a1c0;
937 else
938 return a0c1;
939 }
940}
941
942Int_t AliAnalysisTaskDDMeson::GetCombCh(const Double_t s1, const Double_t s2) const
943{
944 //
945 //get combination of charges
946 //
947 const Int_t combPP = 1;
948 const Int_t combMM = 2;
949 const Int_t combPM = 3;
950
951 if(s1*s2<0){
952 return combPM;
953 }
954 else if(s1>0)
955 return combPP;
956 else
957 return combMM;
958}
959
960//==========================================================================
961
962TLorentzVector AliAnalysisTaskDDMeson::GetKinematics(const Double_t *pa, const Double_t *pb, const Double_t ma, const Double_t mb, Double_t & cts) const
963{
964 //
965 //get kinematics, cts = cos(theta_{#star})
966 //
967 TLorentzVector va, vb, sumv;
968 va.SetXYZM(pa[0], pa[1], pa[2], ma);
969 vb.SetXYZM(pb[0], pb[1], pb[2], mb);
970 sumv = va+vb;
971
972 const TVector3 bv = -sumv.BoostVector();
973
974 va.Boost(bv);
975 vb.Boost(bv);
976 //printf("test a %f %f %f -- %f %f\n", va.X(), va.Y(), va.Z(), va.Theta(), va.Phi());
977 //printf("test b %f %f %f -- %f %f\n", vb.X(), vb.Y(), vb.Z(), vb.Theta(), vb.Phi());
978
979 cts = fabs(cos(va.Theta()));
980
981 return sumv;
982}
983
984Double_t AliAnalysisTaskDDMeson::GetCtlab(const Double_t *pa, const Double_t *pb) const
985{
986 //
987 //cos(theat_lab)
988 //
989
990 TVector3 va, vb;
991 va.SetXYZ(pa[0], pa[1], pa[2]);
992 vb.SetXYZ(pb[0], pb[1], pb[2]);
993
994 const TVector3 ua = va.Unit();
995 const TVector3 ub = vb.Unit();
996
997 return ua.Dot(ub);
998}