PMD/AliPMDHitsRead.C

   1 //
   2 // This macro reads the Hits Tree
   3 //
   4 void AliPMDHitsRead(Int_t nevt = 1)
   5 {
   6
   7   TStopwatch timer;
   8   timer.Start();
   9   TH2F *h2 = new TH2F("h2"," Y vs. X",200,-100.,100.,200,-100.,100.);
  10   //  FILE *fpw = fopen("alipmdhits.dat","w");
  11
  12   AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root");
  13
  14   if (!fRunLoader)
  15    {
  16      printf("Can not open session for file ");
  17    }
  18
  19   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
  20   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
  21   if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
  22
  23   gAlice = fRunLoader->GetAliRun();
  24
  25   if (gAlice)
  26     {
  27       printf("Alirun object found\n");
  28     }
  29   else
  30     {
  31       printf("Could not found Alirun object\n");
  32     }
  33
  34   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
  35   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
  36   if (fPMDLoader == 0x0)
  37     {
  38       printf("Can not find PMDLoader\n");
  39     }
  40
  41
  42   fPMDLoader->LoadHits("READ");
  43
  44   // This reads the PMD Hits tree and assigns the right track number
  45   // to a cell and stores in the summable digits tree
  46   //
  47
  48   const Int_t kPi0 = 111;
  49   const Int_t kGamma = 22;
  50   Int_t   npmd;
  51   Int_t   trackno;
  52   Int_t   smnumber;
  53   Int_t   trackpid;
  54   Int_t   mtrackno;
  55   Int_t   mtrackpid;
  56   Int_t   xpad = -1, ypad = -1;
  57   Float_t edep;
  58   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
  59   Float_t xPos, yPos, zPos;
  60   Float_t xx, yy;
  61
  62   AliPMDUtility cc;
  63
  64   for (Int_t ievt = 0; ievt < nevt; ievt++)
  65     {
  66
  67       printf("Event Number = %d\n",ievt);
  68       Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
  69       printf("Number of Particles = %d\n",nparticles);
  70       fRunLoader->GetEvent(ievt);
  71       // ------------------------------------------------------- //
  72       // Pointer to specific detector hits.
  73       // Get pointers to Alice detectors and Hits containers
  74
  75       TTree* treeH = fPMDLoader->TreeH();
  76
  77       Int_t ntracks    = (Int_t) treeH->GetEntries();
  78       printf("Number of Tracks in the TreeH = %d\n", ntracks);
  79
  80
  81       TClonesArray* hits = 0;
  82       if (fPMD) hits = fPMD->Hits();
  83
  84       // Start loop on tracks in the hits containers
  85
  86       for (Int_t track=0; track<ntracks;track++)
  87         {
  88           gAlice->ResetHits();
  89           treeH->GetEvent(track);
  90           if (fPMD)
  91             {
  92               npmd = hits->GetEntriesFast();
  93               for (int ipmd = 0; ipmd < npmd; ipmd++)
  94                 {
  95                   fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
  96                   trackno = fPMDHit->GetTrack();
  97
  98                   //fprintf(fpw,"trackno = %d\n",trackno);
  99
 100                   //  get kinematics of the particles
 101
 102                   TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
 103                   trackpid  = mparticle->GetPdgCode();
 104
 105                   Int_t igatr = -999;
 106                   Int_t ichtr = -999;
 107                   Int_t igapid = -999;
 108                   Int_t imo;
 109                   Int_t igen = 0;
 110                   Int_t idmo = -999;
 111
 112                   Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 113                   if (mparticle->GetFirstMother() == -1)
 114                     {
 115                       tracknoOld  = trackno;
 116                       trackpidOld = trackpid;
 117                       statusOld   = -1;
 118
 119                       vx = mparticle->Vx();
 120                       vy = mparticle->Vy();
 121                       vz = mparticle->Vz();
 122
 123                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, -1, trackpid, vx, vy, vz);
 124
 125                     }
 126                   Int_t igstatus = 0;
 127                   while((imo = mparticle->GetFirstMother()) >= 0)
 128                     {
 129                       igen++;
 130
 131                       mparticle =  gAlice->GetMCApp()->Particle(imo);
 132                       idmo = mparticle->GetPdgCode();
 133
 134                       vx = mparticle->Vx();
 135                       vy = mparticle->Vy();
 136                       vz = mparticle->Vz();
 137
 138                       //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 139                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 140
 141                       if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
 142                         {
 143                           igatr = imo;
 144                           igapid = idmo;
 145                           igstatus = 1;
 146                         }
 147                       if(igstatus == 0)
 148                         {
 149                           if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 150                             {
 151                               igatr = imo;
 152                               igapid = idmo;
 153                             }
 154                         }
 155                       ichtr = imo;
 156                     } // end of while loop
 157
 158                   if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 159                     {
 160                       mtrackno = igatr;
 161                       mtrackpid = igapid;
 162                     }
 163                   else
 164                     {
 165                       mtrackno  = ichtr;
 166                       mtrackpid = idmo;
 167                     }
 168                   if (statusOld == -1)
 169                     {
 170                       mtrackno  = tracknoOld;
 171                       mtrackpid = trackpidOld;
 172                     }
 173
 174                   xPos = fPMDHit->X();
 175                   yPos = fPMDHit->Y();
 176                   zPos = fPMDHit->Z();
 177
 178                   edep       = fPMDHit->GetEnergy();
 179                   Int_t vol1 = fPMDHit->GetVolume(1); // Column
 180                   Int_t vol2 = fPMDHit->GetVolume(2); // Row
 181                   Int_t vol3 = fPMDHit->GetVolume(7); // UnitModule
 182                   Int_t vol6 = fPMDHit->GetVolume(8); // SuperModule
 183
 184                   // -----------------------------------------//
 185                   // For Super Module 1 & 2                   //
 186                   //  nrow = 96, ncol = 48                    //
 187                   // For Super Module 3 & 4                   //
 188                   //  nrow = 48, ncol = 96                    //
 189                   // -----------------------------------------//
 190
 191                   smnumber = (vol6-1)*6 + (vol3-1);
 192
 193                   xpad = vol1 - 1;
 194                   ypad = vol2 - 1;
 195
 196                   if(zPos > 361.5)
 197                     {
 198                       cc.RectGeomCellPos(smnumber,xpad,ypad,xx,yy);
 199                       h2->Fill(xx,yy);
 200                     }
 201
 202                 }
 203             }
 204         } // Track Loop ended
 205
 206     }
 207
 208   h2->Draw();
 209
 210   fRunLoader->UnloadgAlice();
 211   fRunLoader->UnloadHeader();
 212   fRunLoader->UnloadKinematics();
 213   fPMDLoader->UnloadHits();
 214
 215   timer.Stop();
 216   timer.Print();
 217
 218 }