CORRFW/test/muons/AliCFMuonResTask1.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, 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 // Example of task (running locally, on AliEn and CAF),\r
  18 // which provides standard way of calculating acceptance and efficiency\r
  19 // between different steps of the procedure.\r
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies\r
  21 // can be calculated\r
  22 //-----------------------------------------------------------------------\r
  23 // Author : R. Vernet, Consorzio Cometa - Catania (it)\r
  24 //-----------------------------------------------------------------------\r
  25 // Modification done by X. Lopez - LPC Clermont (fr)\r
  26 //-----------------------------------------------------------------------\r
  27 \r
  28 \r
  29 #ifndef ALICFMUONRESTASK1_CXX\r
  30 #define ALICFMUONRESTASK1_CXX\r
  31 \r
  32 #include "AliCFMuonResTask1.h"\r
  33 #include "AliHeader.h"\r
  34 #include "AliESDHeader.h"\r
  35 #include "AliStack.h"\r
  36 #include "TParticle.h"\r
  37 #include "TH1I.h"\r
  38 #include "AliMCEvent.h"\r
  39 #include "AliAnalysisManager.h"\r
  40 #include "AliESDEvent.h"\r
  41 #include "AliAODEvent.h"\r
  42 #include "AliCFManager.h"\r
  43 #include "AliCFCutBase.h"\r
  44 #include "AliCFContainer.h"\r
  45 #include "TChain.h"\r
  46 #include "AliESDtrack.h"\r
  47 #include "AliLog.h"\r
  48 #include "AliESDMuonTrack.h"\r
  49 #include "AliESDtrack.h"\r
  50 #include "AliESDInputHandler.h"\r
  51 #include "TCanvas.h"\r
  52 \r
  53 ClassImp(AliCFMuonResTask1)\r
  54 \r
  55 //__________________________________________________________________________\r
  56 AliCFMuonResTask1::AliCFMuonResTask1() :\r
  57   fReadAODData(0),\r
  58   fCFManager(0x0),\r
  59   fQAHistList(0x0),\r
  60   fHistEventsProcessed(0x0),\r
  61   fNevt(0)\r
  62 {\r
  63   //\r
  64   //Default ctor\r
  65   //\r
  66 }\r
  67 //___________________________________________________________________________\r
  68 AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) :\r
  69   AliAnalysisTaskSE(name),\r
  70   fReadAODData(0),\r
  71   fCFManager(0x0),\r
  72   fQAHistList(0x0),\r
  73   fHistEventsProcessed(0x0),\r
  74   fNevt(0)\r
  75 {\r
  76   //\r
  77   // Constructor. Initialization of Inputs and Outputs\r
  78   //\r
  79   Info("AliCFMuonResTask1","Calling Constructor");\r
  80 \r
  81   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;\r
  82 \r
  83   DefineOutput(1,TH1I::Class());\r
  84   DefineOutput(2,AliCFContainer::Class());\r
  85 \r
  86 }\r
  87 \r
  88 //___________________________________________________________________________\r
  89 AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) \r
  90 {\r
  91   //\r
  92   // Assignment operator\r
  93   //\r
  94   if (this!=&c) {\r
  95     AliAnalysisTaskSE::operator=(c) ;\r
  96     fReadAODData = c.fReadAODData ;\r
  97     fCFManager  = c.fCFManager;\r
  98     fQAHistList = c.fQAHistList ;\r
  99     fHistEventsProcessed = c.fHistEventsProcessed;\r
 100     fNevt = c.fNevt ;\r
 101   }\r
 102   return *this;\r
 103 }\r
 104 \r
 105 //___________________________________________________________________________\r
 106 AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) :\r
 107   AliAnalysisTaskSE(c),\r
 108   fReadAODData(c.fReadAODData),\r
 109   fCFManager(c.fCFManager),\r
 110   fQAHistList(c.fQAHistList),\r
 111   fHistEventsProcessed(c.fHistEventsProcessed),\r
 112   fNevt(c.fNevt) \r
 113 {\r
 114   //\r
 115   // Copy Constructor\r
 116   //\r
 117 }\r
 118 \r
 119 //___________________________________________________________________________\r
 120 AliCFMuonResTask1::~AliCFMuonResTask1() {\r
 121   //\r
 122   //destructor\r
 123   //\r
 124   Info("~AliCFMuonResTask1","Calling Destructor");\r
 125   if (fCFManager)           delete fCFManager ;\r
 126   if (fHistEventsProcessed) delete fHistEventsProcessed ;\r
 127   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}\r
 128 }\r
 129 \r
 130 //_________________________________________________\r
 131 void AliCFMuonResTask1::UserExec(Option_t *)\r
 132 {\r
 133   //\r
 134   // Main loop function\r
 135   // \r
 136   \r
 137   Info("UserExec","") ;\r
 138   if (!fMCEvent) {\r
 139     Error("UserExec","NO MC EVENT FOUND!");\r
 140     return;\r
 141   }\r
 142 \r
 143   fNevt++; \r
 144   Double_t containerInput[9] ;\r
 145  \r
 146 ////////\r
 147 //// MC\r
 148 //////// \r
 149 \r
 150   fCFManager->SetEventInfo(fMCEvent);\r
 151   AliStack* stack = fMCEvent->Stack();\r
 152 \r
 153   // loop on the MC event\r
 154   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { \r
 155     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);\r
 156  \r
 157     TParticle *part = mcPart->Particle(); \r
 158     TParticle *part0 = mcPart->Particle();\r
 159     TParticle *part1 = mcPart->Particle();\r
 160  \r
 161     // Selection of the resonance\r
 162     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;\r
 163 \r
 164     // Mother kinematics\r
 165     Float_t e = part->Energy();\r
 166     Float_t pz = part->Pz();           \r
 167     Float_t py = part->Py();\r
 168     Float_t px = part->Px();\r
 169     Float_t rapmc = Rap(e,pz);\r
 170     Float_t mass = part->GetCalcMass();\r
 171 \r
 172     // Decays kinematics\r
 173 \r
 174     Int_t p0 = part->GetDaughter(0);\r
 175     part0 = stack->Particle(p0); \r
 176    // selection of the rapidity for first muon\r
 177     AliMCParticle *mcpart0 = new AliMCParticle(part0);\r
 178     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart0)) continue;\r
 179     Int_t pdg0 = part0->GetPdgCode();\r
 180 \r
 181     Int_t p1 = part->GetDaughter(1);\r
 182     part1 = stack->Particle(p1);\r
 183    // selection of the rapidity for second muon\r
 184     AliMCParticle *mcpart1 = new AliMCParticle(part1);\r
 185     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart1)) continue;\r
 186     Int_t pdg1 = part1->GetPdgCode();\r
 187  \r
 188    // 0 mu- 1 mu+\r
 189     Float_t e0, pz0, py0, px0, phi0, rapmc0;\r
 190     Float_t e1, pz1, py1, px1, phi1, rapmc1;\r
 191 \r
 192     // ordering sign: first = mu-\r
 193     if(pdg0==13){\r
 194         e0 = part0->Energy();\r
 195         pz0 = part0->Pz();\r
 196         py0 = part0->Py();\r
 197         px0 = part0->Px();\r
 198         phi0 = part0->Phi(); // Warning in TParticle Phi = pi + ATan2(Py,Px) = [0,2pi] \r
 199         phi0 = Phideg(phi0);    \r
 200         rapmc0=Rap(e0,pz0);\r
 201 \r
 202         e1 = part1->Energy();\r
 203         pz1 = part1->Pz();\r
 204         py1 = part1->Py();\r
 205         px1 = part1->Px();\r
 206         phi1 = part1->Phi();\r
 207         phi1 = Phideg(phi1);\r
 208         rapmc1=Rap(e1,pz1);\r
 209     }\r
 210     else{\r
 211         if(pdg0==-13){\r
 212             e1 = part0->Energy();\r
 213             pz1 = part0->Pz();\r
 214             py1 = part0->Py();\r
 215             px1 = part0->Px();\r
 216             phi1 = part0->Phi();\r
 217             phi1 = Phideg(phi1);    \r
 218             rapmc1=Rap(e1,pz1);\r
 219             \r
 220             e0 = part1->Energy();\r
 221             pz0 = part1->Pz();\r
 222             py0 = part1->Py();\r
 223             px0 = part1->Px();\r
 224             phi0 = part1->Phi();\r
 225             phi0 = Phideg(phi0);\r
 226             rapmc0=Rap(e0,pz0); \r
 227         }\r
 228     }\r
 229     \r
 230     if(pdg0==13 || pdg1==13) { \r
 231 \r
 232         Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+\r
 233                                    (pz0+pz1)*(pz0+pz1));\r
 234         Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));\r
 235         Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);\r
 236         Float_t rapmc_check=Rap((e0+e1),(pz0+pz1));\r
 237 \r
 238         containerInput[0] = fNevt ;   \r
 239         containerInput[1] = rapmc0 ;   \r
 240         containerInput[2] = phi0 ;   \r
 241         containerInput[3] = rapmc1 ;   \r
 242         containerInput[4] = phi1 ;   \r
 243         containerInput[5] = imassmc ;   \r
 244         containerInput[6] = rapmc ;   \r
 245         containerInput[7] = ptmc;\r
 246         containerInput[8] = pmc;        \r
 247 \r
 248         // fill the container at the first step\r
 249         fCFManager->GetParticleContainer()->Fill(containerInput,0);\r
 250     }\r
 251   }    \r
 252 \r
 253 ////////\r
 254 //// ESD\r
 255 ////////\r
 256   \r
 257   AliESDEvent *fESD; \r
 258   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>\r
 259       (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
 260   fESD = esdH->GetEvent();\r
 261 \r
 262   Int_t mult1 = fESD->GetNumberOfMuonTracks() ;\r
 263 \r
 264     for (Int_t j = 0; j<mult1; j++) { \r
 265         AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));\r
 266         Float_t zr1 = mu1->Charge();\r
 267 // Select mu-\r
 268         if(zr1<0){\r
 269             Float_t pxr1 = mu1->Px();\r
 270             Float_t pyr1 = mu1->Py();\r
 271             Float_t pzr1 = mu1->Pz();\r
 272             Float_t phir1 = mu1->Phi(); // Warning in AliESDMuonTrack Phi = ATan2(Py,Px) = [-pi,pi]\r
 273             phir1 = phir1 * 57.296;\r
 274             Float_t er1 = mu1->E();\r
 275             Float_t rap1=Rap(er1,pzr1);\r
 276             // selection of the rapidity for first muon\r
 277             if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mu1)) continue;\r
 278 \r
 279             for (Int_t jj = 0; jj<mult1; jj++) {\r
 280                 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));\r
 281                 Float_t zr2 = mu2->Charge();\r
 282 // Select mu+\r
 283                 if(zr2>0 && jj !=j){\r
 284                     Float_t pxr2 = mu2->Px();\r
 285                     Float_t pyr2 = mu2->Py();\r
 286                     Float_t pzr2 = mu2->Pz();\r
 287                     Float_t phir2 = mu2->Phi();\r
 288                     phir2 = phir2 * 57.296;\r
 289                     Float_t er2 = mu2->E();\r
 290                     Float_t rap2=Rap(er2,pzr2);\r
 291                     // selection of the rapidity for second muon\r
 292                     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mu2)) continue;\r
 293 \r
 294                     Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+\r
 295                                                 (pzr1+pzr2)*(pzr1+pzr2));\r
 296                     Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));\r
 297                     Float_t raprec=Rap((er1+er2),(pzr1+pzr2));\r
 298                     Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);\r
 299 \r
 300                     containerInput[0] = fNevt ;   \r
 301                     containerInput[1] = rap1 ;   \r
 302                     containerInput[2] = phir1 ;   \r
 303                     containerInput[3] = rap2 ;\r
 304                     containerInput[4] = phir2 ;   \r
 305                     containerInput[5] = imassrec ;   \r
 306                     containerInput[6] = raprec ;   \r
 307                     containerInput[7] = ptrec;\r
 308                     containerInput[8] = prec;\r
 309                     \r
 310                     // fill the container at the second step\r
 311                     fCFManager->GetParticleContainer()->Fill(containerInput,1);\r
 312 \r
 313                 }  // mu+ Selection\r
 314             }      // second mu Loop\r
 315         }          // mu- Selection\r
 316     }        \r
 317 \r
 318 //  ----------\r
 319   fHistEventsProcessed->Fill(0);\r
 320   PostData(1,fHistEventsProcessed) ;\r
 321   PostData(2,fCFManager->GetParticleContainer()) ;\r
 322 }\r
 323 //________________________________________________________________________\r
 324 const Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,\r
 325                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) \r
 326 {\r
 327 // invariant mass calculation\r
 328     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+\r
 329                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));\r
 330     return imassrec;\r
 331 }\r
 332 //________________________________________________________________________\r
 333 const Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) \r
 334 {\r
 335 // calculate rapidity\r
 336     Float_t rap;\r
 337     if(e!=pz){\r
 338         rap = 0.5*TMath::Log((e+pz)/(e-pz));\r
 339         return rap;\r
 340     }\r
 341     else{\r
 342         rap = -200;\r
 343         return rap;\r
 344     }\r
 345 }\r
 346 //________________________________________________________________________\r
 347 const Float_t AliCFMuonResTask1::Phideg(Float_t phi) \r
 348 {\r
 349 // calculate Phi from TParticle in range [-180,180] \r
 350     Float_t phideg;\r
 351     \r
 352         phideg = phi-TMath::Pi();\r
 353         phideg = phideg*57.296;\r
 354         return phideg;\r
 355 }\r
 356 //________________________________________________________________________\r
 357 void AliCFMuonResTask1::Terminate(Option_t *) \r
 358 {\r
 359   // draw result of the Invariant mass MC and ESD\r
 360 \r
 361     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));   \r
 362 \r
 363     TH1D *kmass = cont->ShowProjection(5,0);\r
 364     TH1D *rmass = cont->ShowProjection(5,1);\r
 365 \r
 366     TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);\r
 367     c1->Divide(1,2);\r
 368     c1->cd(1);\r
 369     kmass->Draw("HIST");\r
 370     c1->cd(2);\r
 371     rmass->Draw("HIST");\r
 372 }\r
 373 //________________________________________________________________________\r
 374 \r
 375 #endif\r