PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraBoth.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2009, 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 //         AliAnalysisTaskSpectraBoth class\r
  18 //-----------------------------------------------------------------\r
  19 \r
  20 #include "TChain.h"\r
  21 #include "TTree.h"\r
  22 #include "TLegend.h"\r
  23 #include "TH1F.h"\r
  24 #include "TH2F.h"\r
  25 #include "TCanvas.h"\r
  26 #include "AliAnalysisTask.h"\r
  27 #include "AliAnalysisManager.h"\r
  28 #include "AliAODTrack.h"\r
  29 #include "AliAODMCParticle.h"\r
  30 #include "AliVParticle.h"\r
  31 #include "AliAODEvent.h"\r
  32 #include "AliAODInputHandler.h"\r
  33 #include "AliAnalysisTaskSpectraBoth.h"\r
  34 #include "AliAnalysisTaskESDfilter.h"\r
  35 #include "AliAnalysisDataContainer.h"\r
  36 #include "AliSpectraBothHistoManager.h"\r
  37 #include "AliSpectraBothTrackCuts.h"\r
  38 #include "AliSpectraBothEventCuts.h"\r
  39 #include "AliCentrality.h"\r
  40 #include "TProof.h"\r
  41 #include "AliPID.h"\r
  42 #include "AliVEvent.h"\r
  43 #include "AliESDEvent.h"\r
  44 #include "AliPIDResponse.h"\r
  45 #include "AliStack.h"\r
  46 #include "AliSpectraBothPID.h"\r
  47 #include <TMCProcess.h>\r
  48 \r
  49 #include <iostream>\r
  50 \r
  51 \r
  52 \r
  53 \r
  54 using namespace AliSpectraNameSpaceBoth;\r
  55 using namespace std;\r
  56 \r
  57 ClassImp(AliAnalysisTaskSpectraBoth)\r
  58 \r
  59 //________________________________________________________________________\r
  60 AliAnalysisTaskSpectraBoth::AliAnalysisTaskSpectraBoth(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0),  fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
  61 {\r
  62   // Default constructor\r
  63   \r
  64   DefineInput(0, TChain::Class());\r
  65   DefineOutput(1, AliSpectraBothHistoManager::Class());\r
  66   DefineOutput(2, AliSpectraBothEventCuts::Class());\r
  67   DefineOutput(3, AliSpectraBothTrackCuts::Class());\r
  68   DefineOutput(4, AliSpectraBothPID::Class());\r
  69   fNRebin=0;\r
  70   \r
  71 }\r
  72 //________________________________________________________________________\r
  73 //________________________________________________________________________\r
  74 void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()\r
  75 {\r
  76   // create output objects\r
  77   fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin);\r
  78 \r
  79   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
  80   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
  81   if (!fPID)       AliFatal("PID object should be set in the steering macro");\r
  82   fTrackCuts->SetAliESDtrackCuts(fCuts);\r
  83   PostData(1, fHistMan  );\r
  84   PostData(2, fEventCuts);\r
  85   PostData(3, fTrackCuts);\r
  86   PostData(4, fPID      );\r
  87 \r
  88 }\r
  89 //________________________________________________________________________\r
  90 void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)\r
  91 {\r
  92   // main event loop\r
  93         Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;\r
  94         fAOD = dynamic_cast<AliVEvent*>(InputEvent());\r
  95    \r
  96   //    AliESDEvent* esdevent=0x0;\r
  97  //     AliAODEvent* aodevent=0x0;\r
  98    \r
  99         TString nameoftrack(fAOD->ClassName());  \r
 100         if(!nameoftrack.CompareTo("AliESDEvent"))\r
 101         {\r
 102                 ifAODEvent=AliSpectraBothTrackCuts::kESDobject;\r
 103                 //esdevent=dynamic_cast<AliESDEvent*>(fAOD);\r
 104         }\r
 105         else if(!nameoftrack.CompareTo("AliAODEvent"))\r
 106         {\r
 107                 ifAODEvent=AliSpectraBothTrackCuts::kAODobject;\r
 108                 //aodevent=dynamic_cast<AliAODEvent*>(fAOD);\r
 109         }\r
 110         else\r
 111                 AliFatal("Not processing AODs or ESDS") ;\r
 112         if(fdotheMCLoopAfterEventCuts)\r
 113                 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
 114  \r
 115   \r
 116         TClonesArray *arrayMC = 0;\r
 117         Int_t npar=0;\r
 118         AliStack* stack=0x0;\r
 119   \r
 120         if (fIsMC)\r
 121         {\r
 122           \r
 123                   if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
 124                   {     \r
 125                           arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 126                           if (!arrayMC) \r
 127                           {\r
 128                                         AliFatal("Error: MC particles branch not found!\n");\r
 129                           }\r
 130                           Int_t nMC = arrayMC->GetEntries();\r
 131                           for (Int_t iMC = 0; iMC < nMC; iMC++)\r
 132                           {\r
 133                                   AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
 134                                   if(!partMC->Charge()) continue;//Skip neutrals\r
 135                                   //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance\r
 136                                   //Printf("%f     %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());\r
 137                                   if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
 138                                                 fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
 139                                   \r
 140                                   //rapidity cut\r
 141                                   if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) continue;        \r
 142                                   if(partMC->IsPhysicalPrimary())\r
 143                                          npar++;    \r
 144                                   // check for true PID + and fill P_t histos \r
 145                                   Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
 146                                   Int_t id = fPID->GetParticleSpecie(partMC);\r
 147                                   if(id != kSpUndefined) \r
 148                                   {\r
 149                                         fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
 150                                   }\r
 151                           }\r
 152                   }\r
 153                   if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
 154                   {     \r
 155                         AliMCEvent* mcEvent  = (AliMCEvent*) MCEvent();\r
 156                         //Printf("MC particles: %d", mcEvent->GetNumberOfTracks());             \r
 157                         if (!mcEvent) \r
 158                         {\r
 159                                 AliFatal("Error: MC particles branch not found!\n");\r
 160                         }\r
 161                         stack = mcEvent->Stack();\r
 162                         Int_t nMC = stack->GetNtrack();\r
 163                         for (Int_t iMC = 0; iMC < nMC; iMC++)\r
 164                         {\r
 165                                  \r
 166                                 TParticle *partMC = stack->Particle(iMC);\r
 167                                 \r
 168                                 if(!partMC)     \r
 169                                         continue;       \r
 170         \r
 171                                 if(!partMC->GetPDG(0))\r
 172                                         continue;\r
 173                                 if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) \r
 174                                         continue;//Skip neutrals\r
 175                                 if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
 176                                         fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
 177                                 if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) \r
 178                                         continue;\r
 179                                 if(stack->IsPhysicalPrimary(iMC))\r
 180                                          npar++;    \r
 181                                   // check for true PID + and fill P_t histos \r
 182                                 Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;\r
 183                                 Int_t id = fPID->GetParticleSpecie(partMC);\r
 184                                 if(id != kSpUndefined) \r
 185                                 {\r
 186                                         fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
 187                                 }\r
 188                           }\r
 189                   }\r
 190         }\r
 191   \r
 192         if(!fdotheMCLoopAfterEventCuts)\r
 193                 if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
 194 \r
 195         //main loop on tracks\r
 196         Int_t ntracks=0;\r
 197         //cout<<fAOD->GetNumberOfTracks()<<endl;\r
 198         for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) \r
 199         {\r
 200                 AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));\r
 201                 AliAODTrack* aodtrack=0;\r
 202                 AliESDtrack* esdtrack=0;\r
 203                 Float_t dca=-999.;\r
 204                 if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
 205                 {\r
 206                         esdtrack=dynamic_cast<AliESDtrack*>(track);\r
 207                         Float_t dcaz=0.0;\r
 208                         esdtrack->GetImpactParameters(dca,dcaz);\r
 209                 }\r
 210                 else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
 211                 {\r
 212                         aodtrack=dynamic_cast<AliAODTrack*>(track);     \r
 213                         dca=aodtrack->DCA();\r
 214                 }\r
 215                 else\r
 216                         continue;\r
 217                 if (!fTrackCuts->IsSelected(track,kTRUE)) \r
 218                         continue;\r
 219                 ntracks++;\r
 220                 fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
 221     \r
 222                 //calculate DCA for AOD track\r
 223                 if(dca==-999.)\r
 224                 {// track->DCA() does not work in old AOD production\r
 225                         Double_t d[2], covd[3];\r
 226                         AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
 227                         Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
 228                         delete track_clone;\r
 229                         if(!isDCA)\r
 230                                 d[0]=-999.;\r
 231                         dca=d[0];\r
 232                 }\r
 233                 fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca);  // PT histo\r
 234                 // get identity and charge\r
 235                 Bool_t rec[3]={false,false,false};\r
 236                 Int_t idRec  = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);\r
 237                 for(int irec=kSpPion;irec<kNSpecies;irec++)\r
 238                 {\r
 239    \r
 240                         if(fUseMinSigma)\r
 241                         {\r
 242                                 if(irec>kSpPion)\r
 243                                         break;\r
 244                         }\r
 245                         else\r
 246                         {       \r
 247                                 if(!rec[irec]) \r
 248                                         idRec = kSpUndefined;\r
 249                                 else    \r
 250                                         idRec=irec;\r
 251                         }               \r
 252    \r
 253                         Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
 254                 \r
 255                         // Fill histograms, only if inside y and nsigma acceptance\r
 256                         if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))\r
 257                                 fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);\r
 258                         //can't put a continue because we still have to fill allcharged primaries, done later\r
 259                 \r
 260                         /* MC Part */\r
 261                         if (arrayMC||stack) \r
 262                         {\r
 263                                 Bool_t isPrimary           = kFALSE;\r
 264                                 Bool_t isSecondaryMaterial = kFALSE; \r
 265                                 Bool_t isSecondaryWeak     = kFALSE; \r
 266                                 Int_t idGen     =kSpUndefined;\r
 267                                 Int_t pdgcode=0;\r
 268                                 if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
 269                                 {\r
 270                                         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
 271                                         if (!partMC) \r
 272                                         { \r
 273                                                 AliError("Cannot get MC particle");\r
 274                                                 continue; \r
 275                                         }\r
 276                                         // Check if it is primary, secondary from material or secondary from weak decay\r
 277                                         isPrimary           = partMC->IsPhysicalPrimary();\r
 278                                         isSecondaryWeak     = partMC->IsSecondaryFromWeakDecay();\r
 279                                         isSecondaryMaterial      = partMC->IsSecondaryFromMaterial();\r
 280                                         //cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
 281 \r
 282                                         if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial) \r
 283                                         {\r
 284                                                 AliError("old tagging");\r
 285                                                 Int_t mfl=-999,codemoth=-999;\r
 286                                                 Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
 287                                                 if(indexMoth>=0)\r
 288                                                 {//is not fake\r
 289                                                         AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
 290                                                         codemoth = TMath::Abs(moth->GetPdgCode());\r
 291                                                         mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
 292                                                 }\r
 293                                                 //Int_t uniqueID = partMC->GetUniqueID();\r
 294                                                 //cout<<"uniqueID: "<<partMC->GetUniqueID()<<"       "<<kPDecay<<endl;\r
 295                                                 //cout<<"status: "<<partMC->GetStatus()<<"       "<<kPDecay<<endl;\r
 296                                                 // if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");\r
 297                                                 if(mfl==3) \r
 298                                                         isSecondaryWeak     = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
 299                                                 else       \r
 300                                                         isSecondaryMaterial = kTRUE;\r
 301                                         }\r
 302                                         //cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
 303 \r
 304 \r
 305                                         idGen     = fPID->GetParticleSpecie(partMC);\r
 306                                         pdgcode=partMC->GetPdgCode(); \r
 307                                 }\r
 308                                 else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
 309                                 {\r
 310                                         TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));\r
 311                                         if (!partMC) \r
 312                                         { \r
 313                                                 AliError("Cannot get MC particle");\r
 314                                                 continue; \r
 315                                         }\r
 316                                         isPrimary           = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));\r
 317                                         isSecondaryWeak     = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));\r
 318                                         isSecondaryMaterial      = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));\r
 319                                         //cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;\r
 320                                         \r
 321                                         idGen     = fPID->GetParticleSpecie(partMC);\r
 322                                         pdgcode=partMC->GetPdgCode(); \r
 323                                 }\r
 324                                 else\r
 325                                         return;\r
 326                         \r
 327                         //        cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;\r
 328                         //        cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;\r
 329                   \r
 330                                 if (isPrimary&&irec==kSpPion)\r
 331                                         fHistMan->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),dca);  // PT histo of primaries\r
 332                   \r
 333                         //nsigma cut (reconstructed nsigma)\r
 334                                 if(idRec == kSpUndefined) \r
 335                                         continue;\r
 336                   \r
 337                         // rapidity cut (reconstructed pt and identity)\r
 338                                  if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;\r
 339                   \r
 340                   // Get true ID\r
 341                   \r
 342                   //if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) continue;          // FIXME: do we need a rapidity cut on the generated?\r
 343                   // Fill histograms for primaries\r
 344                   \r
 345                                  if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue,  idGen, charge)->Fill(track->Pt(),dca); \r
 346                   \r
 347                                 if (isPrimary) \r
 348                                 {\r
 349                                         fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); \r
 350                                         if(idGen != kSpUndefined) \r
 351                                         {\r
 352                                                 fHistMan->GetHistogram2D(kHistPtRecPrimary,      idGen, charge)->Fill(track->Pt(),dca);\r
 353                                                 if (idRec == idGen) \r
 354                                                         fHistMan->GetHistogram2D(kHistPtRecTruePrimary,  idGen, charge)->Fill(track->Pt(),dca); \r
 355                                         }\r
 356                                 }\r
 357                                  //25th Apr - Muons are added to Pions -- FIXME\r
 358                                 if ( pdgcode == 13 && idRec == kSpPion) \r
 359                                 { \r
 360                                         fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); \r
 361                                         if(isPrimary)\r
 362                                                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); \r
 363                                 }\r
 364                                 if ( pdgcode == -13 && idRec == kSpPion) \r
 365                                 { \r
 366                                         fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); \r
 367                                         if (isPrimary) \r
 368                                         {\r
 369                                                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); \r
 370                                         }\r
 371                                 }\r
 372                   \r
 373                                  ///..... END FIXME\r
 374                   \r
 375                                 // Fill secondaries\r
 376                                 if(isSecondaryWeak    )  \r
 377                                         fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);\r
 378                                 if(isSecondaryMaterial)  \r
 379                                         fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);\r
 380                   \r
 381                         }//end if(arrayMC)\r
 382                 }\r
 383         \r
 384         \r
 385         } // end loop on tracks\r
 386  // cout<< ntracks<<endl;\r
 387   fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);\r
 388   PostData(1, fHistMan  );\r
 389   PostData(2, fEventCuts);\r
 390   PostData(3, fTrackCuts);\r
 391   PostData(4, fPID      );\r
 392 }\r
 393 \r
 394 //_________________________________________________________________\r
 395 void   AliAnalysisTaskSpectraBoth::Terminate(Option_t *)\r
 396 {\r
 397   // Terminate\r
 398 }\r