PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoEventReaderAOD - the reader class for the Alice AOD                //
   4 // Reads in AOD information and converts it into internal AliFemtoEvent       //
   5 // Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl                          //
   6 //          Adam Kisiel kisiel@mps.ohio-state.edu                             //
   7 //                                                                            //
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliFemtoEventReaderAOD.h"
  11
  12 #include "TFile.h"
  13 #include "TTree.h"
  14 #include "AliAODEvent.h"
  15 #include "AliAODTrack.h"
  16 #include "AliAODVertex.h"
  17
  18 #include "AliFmPhysicalHelixD.h"
  19 #include "AliFmThreeVectorF.h"
  20
  21 #include "SystemOfUnits.h"
  22
  23 #include "AliFemtoEvent.h"
  24 #include "AliFemtoModelHiddenInfo.h"
  25
  26 ClassImp(AliFemtoEventReaderAOD)
  27
  28 #if !(ST_NO_NAMESPACES)
  29   using namespace units;
  30 #endif
  31
  32 using namespace std;
  33 //____________________________
  34 //constructor with 0 parameters , look at default settings
  35 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
  36   fNumberofEvent(0),
  37   fCurEvent(0),
  38   fEvent(0x0),
  39   fAllTrue(160),
  40   fAllFalse(160),
  41   fFilterBit(0),
  42   fPWG2AODTracks(0x0),
  43   fInputFile(" "),
  44   fFileName(" "),
  45   fTree(0x0),
  46   fAodFile(0x0)
  47 {
  48   // default constructor
  49   fAllTrue.ResetAllBits(kTRUE);
  50   fAllFalse.ResetAllBits(kFALSE);
  51 }
  52
  53 AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
  54   AliFemtoEventReader(),
  55   fNumberofEvent(0),
  56   fCurEvent(0),
  57   fEvent(0x0),
  58   fAllTrue(160),
  59   fAllFalse(160),
  60   fFilterBit(0),
  61   fPWG2AODTracks(0x0),
  62   fInputFile(" "),
  63   fFileName(" "),
  64   fTree(0x0),
  65   fAodFile(0x0)
  66 {
  67   // copy constructor
  68   fInputFile = aReader.fInputFile;
  69   fFileName  = aReader.fFileName;
  70   fNumberofEvent = aReader.fNumberofEvent;
  71   fCurEvent = aReader.fCurEvent;
  72   fEvent = new AliAODEvent();
  73   fAodFile = new TFile(aReader.fAodFile->GetName());
  74   fAllTrue.ResetAllBits(kTRUE);
  75   fAllFalse.ResetAllBits(kFALSE);
  76   fFilterBit = aReader.fFilterBit;
  77   fPWG2AODTracks = aReader.fPWG2AODTracks;
  78 }
  79 //__________________
  80 //Destructor
  81 AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
  82 {
  83   // destructor
  84   delete fTree;
  85   delete fEvent;
  86   delete fAodFile;
  87   if (fPWG2AODTracks) {
  88     fPWG2AODTracks->Delete();
  89     delete fPWG2AODTracks;
  90   }
  91 }
  92
  93 //__________________
  94 AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
  95 {
  96   // assignment operator
  97   if (this == &aReader)
  98     return *this;
  99
 100   fInputFile = aReader.fInputFile;
 101   fFileName  = aReader.fFileName;
 102   fNumberofEvent = aReader.fNumberofEvent;
 103   fCurEvent = aReader.fCurEvent;
 104   if (fTree) delete fTree;
 105   if (fEvent) delete fEvent;
 106   fEvent = new AliAODEvent();
 107   if (fAodFile) delete fAodFile;
 108   fAodFile = new TFile(aReader.fAodFile->GetName());
 109   fAllTrue.ResetAllBits(kTRUE);
 110   fAllFalse.ResetAllBits(kFALSE);
 111   fFilterBit = aReader.fFilterBit;
 112   fPWG2AODTracks = aReader.fPWG2AODTracks;
 113
 114   return *this;
 115 }
 116 //__________________
 117 AliFemtoString AliFemtoEventReaderAOD::Report()
 118 {
 119   // create reader report
 120   AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
 121   return temp;
 122 }
 123
 124 //__________________
 125 void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
 126 {
 127   //setting the name of file where names of AOD file are written
 128   //it takes only this files which have good trees
 129   char buffer[256];
 130   fInputFile=string(inputFile);
 131   cout<<"Input File set on "<<fInputFile<<endl;
 132   ifstream infile(inputFile);
 133
 134   fTree = new TChain("aodTree");
 135
 136   if(infile.good()==true)
 137     {
 138       //checking if all give files have good tree inside
 139       while (infile.eof()==false)
 140         {
 141           infile.getline(buffer,256);
 142           TFile *aodFile=TFile::Open(buffer,"READ");
 143           if (aodFile!=0x0)
 144             {
 145               TTree* tree = (TTree*) aodFile->Get("aodTree");
 146               if (tree!=0x0)
 147                 {
 148                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
 149                   fTree->AddFile(buffer);
 150                   delete tree;
 151                 }
 152               aodFile->Close();
 153             }
 154           delete aodFile;
 155         }
 156     }
 157 }
 158
 159 AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
 160 {
 161   // read in a next hbt event from the chain
 162   // convert it to AliFemtoEvent and return
 163   // for further analysis
 164   AliFemtoEvent *hbtEvent = 0;
 165
 166   if (fCurEvent==fNumberofEvent)//open next file
 167     {
 168       if(fNumberofEvent==0)
 169         {
 170           fEvent=new AliAODEvent();
 171           fEvent->ReadFromTree(fTree);
 172
 173           // Check for the existence of the additional information
 174           fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
 175
 176           if (fPWG2AODTracks) {
 177             cout << "Found additional PWG2 specific information in the AOD!" << endl;
 178             cout << "Reading only tracks with the additional information" << endl;
 179           }
 180
 181           fNumberofEvent=fTree->GetEntries();
 182           cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
 183           fCurEvent=0;
 184         }
 185       else //no more data to read
 186         {
 187           cout<<"no more files "<<hbtEvent<<endl;
 188           fReaderStatus=1;
 189           return hbtEvent;
 190         }
 191     }
 192
 193   cout<<"starting to read event "<<fCurEvent<<endl;
 194   fTree->GetEvent(fCurEvent);//getting next event
 195   cout << "Read event " << fEvent << " from file " << fTree << endl;
 196
 197   hbtEvent = new AliFemtoEvent;
 198
 199   CopyAODtoFemtoEvent(hbtEvent);
 200
 201   fCurEvent++;
 202
 203   return hbtEvent;
 204 }
 205
 206 void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
 207 {
 208   // A function that reads in the AOD event
 209   // and transfers the neccessary information into
 210   // the internal AliFemtoEvent
 211
 212   // setting global event characteristics
 213   tEvent->SetRunNumber(fEvent->GetRunNumber());
 214   tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
 215   tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 216   tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 217   tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 218   tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 219   tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
 220   tEvent->SetZDCParticipants(-1);
 221   tEvent->SetTriggerMask(fEvent->GetTriggerMask());
 222   tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
 223
 224   // Primary Vertex position
 225   double fV1[3];
 226   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 227
 228   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
 229   tEvent->SetPrimVertPos(vertex);
 230
 231   //starting to reading tracks
 232   int nofTracks=0;  //number of reconstructed tracks in event
 233
 234   // Check to see whether the additional info exists
 235   if (fPWG2AODTracks)
 236     nofTracks=fPWG2AODTracks->GetEntries();
 237   else
 238     nofTracks=fEvent->GetNumberOfTracks();
 239
 240   int realnofTracks=0;   // number of track which we use in a analysis
 241   cout << "Event has " << nofTracks << " tracks " << endl;
 242
 243   for (int i=0;i<nofTracks;i++)
 244     {
 245       AliFemtoTrack* trackCopy = new AliFemtoTrack();
 246
 247       if (fPWG2AODTracks) {
 248         // Read tracks from the additional pwg2 specific AOD part
 249         // if they exist
 250         // Note that in that case all the AOD tracks without the
 251         // additional information will be ignored !
 252         AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
 253
 254         // Getting the AOD track through the ref of the additional info
 255         AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
 256         if (!aodtrack->TestFilterBit(fFilterBit))
 257           continue;
 258
 259         CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
 260
 261         double pxyz[3];
 262         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 263         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 264         // Check the sanity of the tracks - reject zero momentum tracks
 265         if (ktP.mag() == 0) {
 266           delete trackCopy;
 267           continue;
 268         }
 269       }
 270       else {
 271         // No additional information exists
 272         // Read in the normal AliAODTracks
 273         const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
 274
 275         if (!aodtrack->TestFilterBit(fFilterBit))
 276           continue;
 277
 278         CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
 279
 280         double pxyz[3];
 281         aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
 282         const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 283         // Check the sanity of the tracks - reject zero momentum tracks
 284         if (ktP.mag() == 0) {
 285           delete trackCopy;
 286           continue;
 287         }
 288       }
 289
 290
 291       tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
 292       realnofTracks++;//real number of tracks
 293     }
 294
 295   tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
 296
 297   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
 298 }
 299
 300 void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack,
 301                                                  AliFemtoTrack *tFemtoTrack,
 302                                                  AliPWG2AODTrack *tPWG2AODTrack)
 303 {
 304   // Copy the track information from the AOD into the internal AliFemtoTrack
 305   // If it exists, use the additional information from the PWG2 AOD
 306
 307   // Primary Vertex position
 308   double fV1[3];
 309   fEvent->GetPrimaryVertex()->GetPosition(fV1);
 310
 311   tFemtoTrack->SetCharge(tAodTrack->Charge());
 312
 313   //in aliroot we have AliPID
 314   //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
 315   //we use only 5 first
 316
 317   // AOD pid has 10 components
 318   double aodpid[10];
 319   tAodTrack->GetPID(aodpid);
 320   tFemtoTrack->SetPidProbElectron(aodpid[0]);
 321   tFemtoTrack->SetPidProbMuon(aodpid[1]);
 322   tFemtoTrack->SetPidProbPion(aodpid[2]);
 323   tFemtoTrack->SetPidProbKaon(aodpid[3]);
 324   tFemtoTrack->SetPidProbProton(aodpid[4]);
 325
 326   double pxyz[3];
 327   tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
 328   AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
 329   tFemtoTrack->SetP(v);//setting momentum
 330   tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 331   const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
 332   //setting track helix
 333   const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
 334   AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
 335   tFemtoTrack->SetHelix(helix);
 336
 337   // Flags
 338   tFemtoTrack->SetTrackId(tAodTrack->GetID());
 339   tFemtoTrack->SetFlags(1);
 340   tFemtoTrack->SetLabel(tAodTrack->GetLabel());
 341
 342   // Track quality information
 343   float covmat[6];
 344   tAodTrack->GetCovMatrix(covmat);
 345   tFemtoTrack->SetImpactD(covmat[0]);
 346   tFemtoTrack->SetImpactZ(covmat[2]);
 347   tFemtoTrack->SetCdd(covmat[3]);
 348   tFemtoTrack->SetCdz(covmat[4]);
 349   tFemtoTrack->SetCzz(covmat[5]);
 350   // This information is only available in the ESD
 351   // We put in fake values or reasonable estimates
 352   tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
 353   tFemtoTrack->SetITSncls(1);
 354   tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
 355   tFemtoTrack->SetTPCncls(1);
 356   tFemtoTrack->SetTPCnclsF(-1);
 357   tFemtoTrack->SetTPCsignalN(-1);
 358   tFemtoTrack->SetTPCsignalS(-1);
 359
 360   if (tPWG2AODTrack) {
 361     // Copy the PWG2 specific information if it exists
 362     tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
 363     tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
 364
 365     double xtpc[3] = {0,0,0};
 366     tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
 367     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 368     tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
 369     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 370   }
 371   else {
 372     // If not use dummy values
 373     tFemtoTrack->SetTPCClusterMap(fAllTrue);
 374     tFemtoTrack->SetTPCSharedMap(fAllFalse);
 375
 376     double xtpc[3] = {0,0,0};
 377     tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
 378     tFemtoTrack->SetNominalTPCExitPoint(xtpc);
 379   }
 380
 381   int indexes[3];
 382   for (int ik=0; ik<3; ik++) {
 383     indexes[ik] = 0;
 384   }
 385   tFemtoTrack->SetKinkIndexes(indexes);
 386 }
 387
 388 void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
 389 {
 390   fFilterBit = (1 << (ibit));
 391 }
 392
 393
 394
 395
 396
 397