PWGGA/PHOSTasks/omega3pi/AliAnalysisTaskOmegaPi0PiPi.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Boris Polishchuk (Boris.Polishchuk@cern.ch)                    *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 ////////////////////////////////////////////////////////////////////////////
  16 //------------------------------------------------------------------------//
  17 // Class used to do omega(782) -> pi0 pi+ pi- analysis.                   //
  18 //------------------------------------------------------------------------//
  19 ////////////////////////////////////////////////////////////////////////////
  20
  21 #include "AliAnalysisTaskOmegaPi0PiPi.h"
  22 #include "TList.h"
  23 #include "TH1.h"
  24 #include "TH2.h"
  25 #include "AliESDEvent.h"
  26 #include "TRefArray.h"
  27 #include "TLorentzVector.h"
  28 #include "AliESDCaloCluster.h"
  29 #include "AliESDv0.h"
  30 #include "AliESDtrack.h"
  31 #include "AliLog.h"
  32 #include "TString.h"
  33 #include "AliPHOSGeoUtils.h"
  34
  35 ClassImp(AliAnalysisTaskOmegaPi0PiPi)
  36
  37 AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi() :
  38   AliAnalysisTaskSE(),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
  39   fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
  40   fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0),fModules("123")
  41 {
  42   //Default constructor.
  43 }
  44
  45
  46 AliAnalysisTaskOmegaPi0PiPi::AliAnalysisTaskOmegaPi0PiPi(const char* name) :
  47   AliAnalysisTaskSE(name),fOutputContainer(0x0),fhM2piSel(0x0),fhDxy(0x0),fhMggSel(0x0),
  48   fhImpXY(0x0),fhM3pi0to2(0x0),fhM3pi2to4(0x0),fhM3pi4to6(0x0),fhM3pi6to8(0x0),
  49   fhM3pi8to10(0x0),fhM3pi10to12(0x0),fhM3pi12(0x0),fhM3pi0to8(0x0),fModules("123")
  50 {
  51   //Constructor used to create a named task.
  52
  53   DefineOutput(1,TList::Class());
  54 }
  55
  56
  57 AliAnalysisTaskOmegaPi0PiPi::~AliAnalysisTaskOmegaPi0PiPi()
  58 {
  59   //Destructor
  60
  61   if(fOutputContainer){
  62     fOutputContainer->Delete() ;
  63     delete fOutputContainer ;
  64   }
  65 }
  66
  67
  68 void AliAnalysisTaskOmegaPi0PiPi::UserCreateOutputObjects()
  69 {
  70   //Allocates histograms and puts it to the output container.
  71
  72   fOutputContainer = new TList();
  73
  74   fhM2piSel = new TH1F("hM2pi_sel","V0 inv. mass, Dxy<2cm",100,0.,1.);
  75   fOutputContainer->Add(fhM2piSel);
  76
  77   fhDxy = new TH1F("hDxy","Distance from V0 to primary vertex in XY-plane", 500, 0.,500.);
  78   fOutputContainer->Add(fhDxy);
  79
  80   fhMggSel = new TH1F("hMgg_sel","Inv. mass of two clusters, pT>1GeV",100,0.,0.3);
  81   fOutputContainer->Add(fhMggSel);
  82
  83   fhImpXY = new TH1F("hImpXY","Impact parameters in XY-plane",1000,0.,1.);
  84   fOutputContainer->Add(fhImpXY);
  85
  86
  87   //M(3pi) vs pT(3pi), ptrack < 2GeV
  88   fhM3pi0to2 = new TH2F("hM3pi_0_2","pTrack: 0-2GeV",100,0.,10., 200,0.4,1.);
  89   fOutputContainer->Add(fhM3pi0to2);
  90
  91   //M(3pi) vs pT(3pi), 2GeV < ptrack < 4GeV
  92   fhM3pi2to4 = new TH2F("hM3pi_2_4","pTrack: 2-4GeV",100,0.,10., 200,0.4,1.);
  93   fOutputContainer->Add(fhM3pi2to4);
  94
  95   //M(3pi) vs pT(3pi), 4GeV < ptrack < 6GeV
  96   fhM3pi4to6 = new TH2F("hM3pi_4_6","pTrack: 4-6GeV",100,0.,10., 200,0.4,1.);
  97   fOutputContainer->Add(fhM3pi4to6);
  98
  99   //M(3pi) vs pT(3pi), 6GeV < ptrack < 8GeV
 100   fhM3pi6to8 = new TH2F("hM3pi_6_8","pTrack: 6-8GeV",100,0.,10., 200,0.4,1.);
 101   fOutputContainer->Add(fhM3pi6to8);
 102
 103   //M(3pi) vs pT(3pi), 8GeV < ptrack < 10GeV
 104   fhM3pi8to10 = new TH2F("hM3pi_8_10","pTrack: 8-10GeV",100,0.,10., 200,0.4,1.);
 105   fOutputContainer->Add(fhM3pi8to10);
 106
 107   //M(3pi) vs pT(3pi), 10GeV < ptrack < 12GeV
 108   fhM3pi10to12 = new TH2F("hM3pi_10_12","pTrack: 10-12GeV",100,0.,10., 200,0.4,1.);
 109   fOutputContainer->Add(fhM3pi10to12);
 110
 111   //M(3pi) vs pT(3pi), ptrack > 12GeV
 112   fhM3pi12 = new TH2F("hM3pi_12","pTrack>12GeV",100,0.,10., 200,0.4,1.);
 113   fOutputContainer->Add(fhM3pi12);
 114
 115   //M(3pi) vs pT(3pi), ptrack < 8GeV
 116   fhM3pi0to8 = new TH2F("hM3pi_0_8","pTrack: 0-8GeV",100,0.,10., 200,0.4,1.);
 117   fOutputContainer->Add(fhM3pi0to8);
 118
 119 }
 120
 121
 122 void AliAnalysisTaskOmegaPi0PiPi::UserExec(Option_t* /* option */)
 123 {
 124   //Main function that does all the job.
 125
 126   AliVEvent* event = InputEvent();
 127   AliESDEvent* esd = (AliESDEvent*)event ;
 128
 129   Double_t v[3] ; //vertex ;
 130   esd->GetVertex()->GetXYZ(v) ;
 131   TVector3 vtx(v);
 132
 133   AliDebug(2,Form("Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z()));
 134
 135   TRefArray * caloClustersArr  = new TRefArray();
 136   esd->GetPHOSClusters(caloClustersArr);
 137
 138   const Int_t kNumberOfTracks = esd->GetNumberOfTracks();
 139   const Int_t kNumberOfV0s = esd->GetNumberOfV0s();
 140
 141   AliDebug(2,Form("Tracks: %d. V0s: %d\n",kNumberOfTracks,kNumberOfV0s));
 142
 143   Float_t xyImp,zImp,imp1,imp2;
 144
 145   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 146   AliDebug(2,Form("CaloClusters: %d\n", kNumberOfPhosClusters));
 147   if(kNumberOfPhosClusters<2){
 148     delete caloClustersArr;
 149     return;
 150   }
 151   TLorentzVector pc1; //4-momentum of PHOS cluster 1
 152   TLorentzVector pc2; //4-momentum of PHOS cluster 2
 153   TLorentzVector pc12;
 154
 155   TLorentzVector ppos; //4-momentum of positive track
 156   TLorentzVector pneg; //4-momentum of negative track
 157   TLorentzVector ptr12;
 158
 159   TLorentzVector pomega;
 160   Double_t etrack;
 161   Double_t p1,p2; // 3-momentum of tracks
 162
 163   Int_t relid[4] ;
 164   UShort_t *absIds;
 165
 166   Int_t minCells=3; //Ncells>2
 167   const Double_t kMpi = 0.1396;
 168   AliPHOSGeoUtils* geom = new AliPHOSGeoUtils("IHEP","");
 169
 170   for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
 171     AliESDCaloCluster *c1 = (AliESDCaloCluster *) caloClustersArr->At(iClu);
 172     if(c1->GetNCells()<minCells) continue;
 173
 174     absIds = c1->GetCellsAbsId();
 175     geom->AbsToRelNumbering(absIds[0], relid) ;
 176     TString phosMod1;
 177     phosMod1 += relid[0];
 178
 179     if(!fModules.Contains(phosMod1.Data())) continue;
 180     c1->GetMomentum(pc1,v);
 181
 182     for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
 183       AliESDCaloCluster *c2 = (AliESDCaloCluster *) caloClustersArr->At(jClu);
 184       if(c2->IsEqual(c1)) continue;
 185       if(c2->GetNCells()<minCells) continue;
 186
 187       absIds = c2->GetCellsAbsId();
 188       geom->AbsToRelNumbering(absIds[0], relid) ;
 189       TString phosMod2;
 190       phosMod2 += relid[0];
 191
 192       if(!fModules.Contains(phosMod2.Data())) continue;
 193       c2->GetMomentum(pc2,v);
 194
 195       pc12 = pc1+pc2;
 196       AliDebug(2,Form("pc12.M(): %.3f\n",pc12.M()));
 197
 198       if(pc12.M()<0.100 || pc12.M()>0.160) continue; //not a pi0 candidate!
 199       if(pc12.Pt()<1.5) continue; //pT(pi0) > 1.5GeV
 200
 201       fhMggSel->Fill(pc12.M());
 202
 203       for(Int_t iTrack=0; iTrack<kNumberOfTracks; iTrack++) {
 204         AliESDtrack* tr1 = esd->GetTrack(iTrack);
 205         p1 = tr1->P();
 206         tr1->GetImpactParameters(xyImp,zImp);
 207         imp1 = TMath::Abs(xyImp);
 208         fhImpXY->Fill(imp1);
 209         Short_t charge = tr1->Charge();
 210         if(imp1>0.004) continue; // not from the primary vertex!
 211
 212         etrack = TMath::Sqrt(p1*p1 + kMpi*kMpi);
 213         ppos.SetPxPyPzE(tr1->Px(),tr1->Py(),tr1->Pz(),etrack);
 214
 215         for(Int_t jTrack=iTrack; jTrack<kNumberOfTracks; jTrack++) {
 216           AliESDtrack* tr2 = esd->GetTrack(jTrack);
 217           p2 = tr2->P();
 218           if(tr2->Charge()==charge) continue; //same sign track (WRONG!!)
 219           tr2->GetImpactParameters(xyImp,zImp);
 220           imp2=TMath::Abs(xyImp);
 221           if(imp2>0.004) continue; // not from the primary vertex!
 222
 223           etrack = TMath::Sqrt(p2*p2 + kMpi*kMpi);
 224           pneg.SetPxPyPzE(tr2->Px(),tr2->Py(),tr2->Pz(),etrack);
 225
 226           ptr12 = ppos + pneg;
 227
 228 //        printf("ptr12.M()=%.3f, xyImp1=%f, xyImp2=%f, ch1=%d, ch2=%d,
 229 //                  p1=%.3f, p2=%.3f, m1=%.3f, m2=%.3f\n",
 230 //               ptr12.M(),imp1,imp2,tr1->Charge(),tr2->Charge(),
 231 //               tr1->P(),tr2->P(),tr1->M(),tr2->M());
 232
 233           pomega = pc12 + ptr12;
 234           AliDebug(2,Form("pomega.M(): %f\n",pomega.M()));
 235
 236           if( p1<2. && p2<2. )
 237             fhM3pi0to2->Fill(pomega.Pt(),pomega.M());
 238
 239           if( (p1>=2. && p1<4.) && (p2>=2. && p2<4.) )
 240             fhM3pi2to4->Fill(pomega.Pt(),pomega.M());
 241
 242           if( (p1>=4. && p1<6.) && (p2>=4. && p2<6.) )
 243             fhM3pi4to6->Fill(pomega.Pt(),pomega.M());
 244
 245           if( (p1>=6. && p1<8.) && (p2>=6. && p2<8.) )
 246             fhM3pi6to8->Fill(pomega.Pt(),pomega.M());
 247
 248           if( (p1>=8. && p1<10.) && (p2>=8. && p2<10.) )
 249             fhM3pi8to10->Fill(pomega.Pt(),pomega.M());
 250
 251           if( (p1>=10. && p1<12.) && (p2>=10. && p2<12.) )
 252             fhM3pi10to12->Fill(pomega.Pt(),pomega.M());
 253
 254           if( (p1>=12.) && (p2>=12.) )
 255             fhM3pi12->Fill(pomega.Pt(),pomega.M());
 256
 257           if( p1<8. && p2<8. )
 258             fhM3pi0to8->Fill(pomega.Pt(),pomega.M());
 259
 260         }
 261       }
 262
 263       for(Int_t iVert=0; iVert<kNumberOfV0s; iVert++) {
 264         AliESDv0* v0 = esd->GetV0(iVert);
 265         AliESDtrack* ptrack = esd->GetTrack(v0->GetPindex()); //positive track
 266         AliESDtrack* ntrack = esd->GetTrack(v0->GetNindex()); //negative track
 267
 268         etrack = TMath::Sqrt(ptrack->P()*ptrack->P() + kMpi*kMpi);
 269         ppos.SetPxPyPzE(ptrack->Px(),ptrack->Py(),ptrack->Pz(),etrack);
 270
 271         etrack = TMath::Sqrt(ntrack->P()*ntrack->P() + kMpi*kMpi);
 272         pneg.SetPxPyPzE(ntrack->Px(),ntrack->Py(),ntrack->Pz(),etrack);
 273
 274         ptr12 = ppos + pneg;
 275
 276         Double_t dx = vtx.X() - v0->Xv();
 277         Double_t dy = vtx.Y() - v0->Yv();
 278         Double_t dxy = TMath::Sqrt(dx*dx + dy*dy);
 279         fhDxy->Fill(dxy);
 280
 281         AliDebug(2,Form("V0: dxy=%.3f\n",dxy));
 282
 283         if(dxy<2.) fhM2piSel->Fill(ptr12.M());
 284       }
 285
 286     }
 287   }
 288
 289   delete caloClustersArr;
 290   PostData(1,fOutputContainer);
 291 }