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
  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
  57   Float_t xPos, yPos, zPos;
  58   Int_t xpad = -1, ypad = -1;
  59   Float_t edep;
  60   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
  61
  62   for (Int_t ievt = 0; ievt < nevt; ievt++)
  63     {
  64
  65       printf("Event Number = %d\n",ievt);
  66       Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
  67       printf("Number of Particles = %d\n",nparticles);
  68       fRunLoader->GetEvent(ievt);
  69       // ------------------------------------------------------- //
  70       // Pointer to specific detector hits.
  71       // Get pointers to Alice detectors and Hits containers
  72
  73       TTree* treeH = fPMDLoader->TreeH();
  74
  75       Int_t ntracks    = (Int_t) treeH->GetEntries();
  76       printf("Number of Tracks in the TreeH = %d\n", ntracks);
  77
  78
  79       TClonesArray* hits = 0;
  80       if (fPMD) hits = fPMD->Hits();
  81
  82       // Start loop on tracks in the hits containers
  83
  84       for (Int_t track=0; track<ntracks;track++)
  85         {
  86           gAlice->ResetHits();
  87           treeH->GetEvent(track);
  88           if (fPMD)
  89             {
  90               npmd = hits->GetEntriesFast();
  91               for (int ipmd = 0; ipmd < npmd; ipmd++)
  92                 {
  93                   fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
  94                   trackno = fPMDHit->GetTrack();
  95
  96                   //fprintf(fpw,"trackno = %d\n",trackno);
  97
  98                   //  get kinematics of the particles
  99
 100                   TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
 101                   trackpid  = mparticle->GetPdgCode();
 102
 103                   Int_t igatr = -999;
 104                   Int_t ichtr = -999;
 105                   Int_t igapid = -999;
 106                   Int_t imo;
 107                   Int_t igen = 0;
 108                   Int_t idmo = -999;
 109
 110                   Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 111                   if (mparticle->GetFirstMother() == -1)
 112                     {
 113                       tracknoOld  = trackno;
 114                       trackpidOld = trackpid;
 115                       statusOld   = -1;
 116
 117                       vx = mparticle->Vx();
 118                       vy = mparticle->Vy();
 119                       vz = mparticle->Vz();
 120
 121                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, -1, trackpid, vx, vy, vz);
 122
 123                     }
 124                   Int_t igstatus = 0;
 125                   while((imo = mparticle->GetFirstMother()) >= 0)
 126                     {
 127                       igen++;
 128
 129                       mparticle =  gAlice->GetMCApp()->Particle(imo);
 130                       idmo = mparticle->GetPdgCode();
 131
 132                       vx = mparticle->Vx();
 133                       vy = mparticle->Vy();
 134                       vz = mparticle->Vz();
 135
 136                       //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 137                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 138
 139                       if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
 140                         {
 141                           igatr = imo;
 142                           igapid = idmo;
 143                           igstatus = 1;
 144                         }
 145                       if(igstatus == 0)
 146                         {
 147                           if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 148                             {
 149                               igatr = imo;
 150                               igapid = idmo;
 151                             }
 152                         }
 153                       ichtr = imo;
 154                     } // end of while loop
 155
 156                   if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 157                     {
 158                       mtrackno = igatr;
 159                       mtrackpid = igapid;
 160                     }
 161                   else
 162                     {
 163                       mtrackno  = ichtr;
 164                       mtrackpid = idmo;
 165                     }
 166                   if (statusOld == -1)
 167                     {
 168                       mtrackno  = tracknoOld;
 169                       mtrackpid = trackpidOld;
 170                     }
 171
 172                   xPos = fPMDHit->X();
 173                   yPos = fPMDHit->Y();
 174                   zPos = fPMDHit->Z();
 175
 176                   edep       = fPMDHit->GetEnergy();
 177                   Int_t vol1 = fPMDHit->GetVolume(1); // Column
 178                   Int_t vol2 = fPMDHit->GetVolume(2); // Row
 179                   Int_t vol3 = fPMDHit->GetVolume(3); // UnitModule
 180                   Int_t vol6 = fPMDHit->GetVolume(6); // SuperModule
 181
 182                   // -----------------------------------------//
 183                   // For Super Module 1 & 2                   //
 184                   //  nrow = 96, ncol = 48                    //
 185                   // For Super Module 3 & 4                   //
 186                   //  nrow = 48, ncol = 96                    //
 187                   // -----------------------------------------//
 188
 189                   smnumber = (vol6-1)*6 + vol3;
 190
 191                   if (vol6 == 1 || vol6 == 2)
 192                     {
 193                       xpad = vol2;
 194                       ypad = vol1;
 195                     }
 196                   else if (vol6 == 3 || vol6 == 4)
 197                     {
 198                       xpad = vol1;
 199                       ypad = vol2;
 200                     }
 201
 202                   //printf("Zposition = %f Edeposition = %f",zPos,edep);
 203                   Float_t zposition = TMath::Abs(zPos);
 204
 205                 }
 206             }
 207         } // Track Loop ended
 208
 209     }
 210   fRunLoader->UnloadgAlice();
 211   fRunLoader->UnloadHeader();
 212   fRunLoader->UnloadKinematics();
 213
 214   fPMDLoader->UnloadHits();
 215
 216   timer.Stop();
 217   timer.Print();
 218
 219 }